Epigenetic Aging in Major Depressive Disorder

Study participants were from the Netherlands Study of Depression and Anxiety (NESDA), an ongoing longitudinal multicenter cohort study designed to investigate the long-term course and consequences of depressive and anxiety disorders (26). NESDA’s 2,981 participants (ages 18–65) include patients with a current or lifetime diagnosis of depression and/or an anxiety disorder and control subjects with no lifetime depressive disorder or anxiety disorders. Participants were recruited from the general population, general practices, and mental health organizations in order to reflect various settings and the entire range of psychopathology. Presence of major depression was ascertained with the DSM-IV-based Composite International Diagnostic Interview, version 2.1 (CIDI) (27), administered by trained research staff. Exclusion criteria were a clinically overt primary diagnosis of other psychiatric conditions (e.g., psychotic disorders, obsessive-compulsive disorder, bipolar disorder, and severe substance use disorders) and not being fluent in Dutch. The study was approved by the ethical committees of all participating centers, and participants provided written informed consent.

The total sample of 1,130 participants was divided into a control group (no lifetime psychiatric disorders and low depressive symptoms, as indicated by a score <14 on the Inventory of Depressive Symptomatology [28]) (N=319) and a group with current major depressive disorder (a score ≥14 on the Inventory of Depressive Symptomatology within the past 6 months) (N=811). Individuals who did not meet criteria for either of the two groups were excluded. The sample selection was further based on good-quality genome-wide association study genotype information available from a previous investigation (29).

Data on sex, education level (in years), and body mass index (BMI) were collected during interviews. Alcohol use was recorded as mean number of drinks per week. Smoking behavior was measured by plasma cotinine levels, an adequate marker for calculating recent tobacco exposure (30). Physical activity was assessed using the International Physical Activity Questionnaire and indicated by total metabolic equivalent (MET) minutes per week. Health status was assessed as number of chronic diseases for which participants received medical treatment.

In all participants, depression severity was measured with the self-report version of the 30-item Inventory of Depressive Symptomatology (28). Childhood trauma was assessed using the childhood trauma interview from the Netherlands Mental Health Survey and Incidence Study, with personal history questions that included a structured inventory of trauma exposure during childhood. Finally, use of antidepressants was assessed through container inspection and categorized using the World Health Organization’s Anatomical Therapeutic Chemical Classification System: tricyclic antidepressants, selective serotonin reuptake inhibitors, and other antidepressants.

In participants with major depressive disorder, the duration of depression was measured by the Life Chart Interview, utilizing a calendar method to assess the percentage of time in which symptoms were present during the past 4 years (31). Current comorbid anxiety (panic disorder, generalized anxiety disorder, agoraphobia, social phobia) and alcohol disorder, as well as age at onset of depression were assessed with the CIDI. A more detailed description of all phenotypes is provided in the data supplement that accompanies the online edition of this article.

To assay the methylation status of the approximately 28 million common CpG sites in the human genome, we used an optimized protocol for MBD-seq (25). With this approach, genomic DNA is fragmented and the methylated fragments are then bound to the MBD2 protein, which has a high affinity for methylated DNA. The nonmethylated fraction is washed away, and only the methylation-enriched fraction is sequenced (for more details, see the online data supplement). This optimized protocol assesses about 94% of the CpGs in the methylome (25). The sequenced reads were aligned to the reference genome (build hg19/GRCh37) with Bowtie 2 (32) using local and gapped alignment. Aligned reads were further processed using the RaMWAS Bioconductor package (33) to perform quality control and calculate methylation scores for each CpG.

While existing algorithms (9, 10) have gone through demonstrations of utility and reliability, estimating DNAm age with those prediction models would have been suboptimal for this study. These algorithms were derived using methylation data from a different platform in study populations with different characteristics (e.g., age distribution). Commonly used methods for assaying DNA methylation depend on the Illumina arrays, platforms that generate variables representing percentage methylated (ranging from 0 to 1). We used MBD-seq, generating methylation data that is semiquantitative (scores may range from 0–20) (24). Because the weights assigned to individual CpGs when making age predictions directly depend on the platform and study population, they will not optimally capture the effects of CpGs on age in the present study. Therefore, we “recalibrated” the DNAm age estimate in a way that is optimal for this study. It is important to emphasize that we aimed to obtain the best possible DNAm age estimates for MBD-seq data in our sample. We have not developed a new clock to be generalized to data from platforms other than MBD-seq. Furthermore, because we already collected MBD-seq data, we did not attempt to reduce the predictor set to the smallest number of CpGs, as pruning sites may reduce the precision of the DNAm age estimates.

Our approach for estimating DNAm age is similar to the one taken by Horvath (10). Specifically, we used elastic nets, a variable selection method that is particularly useful when the number of predictors is much larger than the number of observations (34). Parameter alpha was set to zero (i.e., ridge regression, retaining all sites in the model) where chronological age was used as outcome and methylation sites as predictors. To estimate predictive power and obtain DNAm age estimates for each subject, k-fold cross-validation was used, with k=10. Thus, the sample was randomly partitioned into 10 equally sized subsamples. Of the 10 subsamples, nine were used as training data and the remaining subsample as validation data. This ensured that in samples with the same properties and platform, our results would “replicate” and provide unbiased estimates of DNAm age. In the RaMWAS implementation, a cycle of methylome-wide association studies (MWAS), marker selection, and estimation via ridge regression was repeated in each training data set, with the resulting model applied to the test data to obtain unbiased estimates of DNAm age for each of the k=10 iterations.

Several analyses were performed to validate the model. First, the model used to estimate DNAm age contained 80,000 CpGs (see Table S1 in the data supplement). Tenfold cross-validation showed that chronological age could be predicted very well, with a correlation of 0.95 (p<0.001) (Figure 1). Second, when analyzing assessed phenotypes in NESDA with DNAm age, we confirmed some similar determinants of DNAm age found in previous studies validating our outcome measure (see Table S3 in the data supplement). Male sex (35) and higher BMI (13, 36, 37) were associated with higher epigenetic aging. Third, to validate calculation of DNAm age, we used both ridge regression and the lasso method (used by Horvath). The additional elastic net model, with parameter alpha set to 0.5, resulted in a comparable correlation of 0.93 between chronological age and predicted DNAm age in our data set (compared with 0.95 with ridge regression and alpha=0), indicating that parameter set point did not have a large impact on our outcome measure. Finally, to ensure that no systematic bias was introduced by training the model on both case subjects and control subjects, we also trained the prediction model in control subjects only. This resulted in a slightly lower correlation between DNAm age estimates and chronological age (r=0.93), since the control subjects represented only a third of our total sample (i.e., lower statistical power). However, the correlation between the DNAm age estimates obtained in the full sample and those obtained using only control subjects was high (r=0.98, p<0.001), indicating that psychiatric status did not have an impact on the estimation of DNAm age.

We pooled data from five brain sample collections from four different brain banks (the Victorian Brain Bank Network, the Harvard Brain Bank, the Netherlands Brain Bank, and the Stanley Medical Research Institute), including a total of 141 brain samples from the Brodmann’s areas 10 and 25 regions. Presence of major depressive disorder (N=74) was determined by at least one psychiatrist by using information obtained from a family member who was well acquainted with the deceased. Control subjects (N=67) had no history of psychiatric disorders. Postmortem interval and pH were recorded in the brain collections, with the exception of the Harvard Brain Bank.

To further test the reliability and validity of our methods, we used the same approach to predict DNAm age in these samples with MBD-seq methylation data generated from the SOLiD 5500 W platform (Life Technologies, Carlsbad, Calif.). The model used to predict DNAm age contained 100,000 CpGs (see Table S2 in the data supplement) and obtained a correlation of 0.69 between predicted DNAm age and chronological age. The lower correlation with age is likely because the methylation data were generated with an older platform with lower quality (e.g., lower alignment of reads) from a more heterogeneous data set. More details about the samples and methods are provided in the data supplement.

To investigate case-control differences in epigenetic aging, we conducted linear regression models with epigenetic aging as the outcome and all covariates as predictors. To correct for the relative abundance of cell types that may be differentially associated with major depression, an additional model included cell type proportions as covariates (38). Other linear regression models were used to examine the relationship between epigenetic aging and Inventory of Depressive Symptomatology score across groups and clinical characteristics within patients with major depression. All analyses were corrected for all sociodemographic, lifestyle, and health covariates and used two-tailed tests with a significance threshold of 0.05.

Within postmortem brain samples, we constructed a linear mixed model in R using the nlme package to account for the heterogeneity of epigenetic aging across brain collections. Thus, brain collection was entered as random effect and sex as fixed effect. The p value was derived by a likelihood ratio test, a hypothesis-driven one-sided test, with p<0.05 considered significant.

To perform enrichment tests of top MWAS findings in brain and blood, we used the shiftR R package with 1 million permutations for each test and used three thresholds (0.5%, 1%, and 5%) to define “top findings.” To account for this “multiple testing,” shiftR uses the same thresholds in the permutations where the test statistic distribution under the null hypothesis is generated from the most significant thresholds or combination of thresholds. A more in-depth description is provided in the data supplement. To gain insight into the overlapping biological pathways affecting epigenetic aging in blood and brain, we used ConsensusPathDB (39) to test whether genes harboring epigenetic aging–associated CpGs were enriched for level 5 Gene Ontology (GO) terms. Methylation sites with p<1×10⁻⁵ were selected and had to be within gene boundaries. At least four genes had to be present in the GO term to be considered. Finally, we also evaluated the overlap between chronological age–associated and epigenetic aging–associated CpG sites to examine whether similar biological processes were involved in chronological and biological aging.

The mean age of the NESDA sample was 41.5 years (SD=13.0, range=18–64), and 64.5% of the sample were female (Table 1). The depression and control groups did not differ significantly in age, but the depression group had a higher proportion of females (p=0.02) and had fewer years of education on average (p<0.001). As anticipated, patients in the depression group reported higher levels of depression severity and use of antidepressants (all p values, <0.001). The mean childhood trauma score was also higher in the depressed group (p<0.001).

Epigenetic aging showed (by design) a mean of zero (SD=3.58 years), ranging from −13.26 to 15.00. Patients in the depression group had significantly higher epigenetic aging compared with the control group (b=0.64, t=2.65, p=0.008; effect size, Cohen’s d=0.18), indicating that the patient was estimated to be 0.64 years (or 7.68 months) older on average than the control group after full adjustment for covariates (Table 2). Additional analyses correcting for cell type proportions did not change results and produced a Cohen’s d of 0.14 (see the data supplement). Consistent with a dose-response effect, a fully adjusted linear regression showed that greater epigenetic aging was significantly associated with higher Inventory of Depressive Symptomatology score in the overall sample (β=0.10, p=0.001). As expected from the high correlation between the DNAm age estimates generated by both models (r=0.98), the above-mentioned results remained unchanged when the same analyses were performed with the DNAm age estimates from the control-subjects-only model (see the data supplement).

Within the depression groups, we found epigenetic aging to be positively associated with childhood trauma score (β=0.09, p=0.02, see Table 3). The association between epigenetic aging and Inventory of Depressive Symptomatology score in the overall sample did not remain significant when analyzed only within the depression group (β=0.05, p=0.21), likely because of lower variation in symptom severity. No further significant associations with clinical characteristics were found.

Further analyses revealed that patients in the depression group who had childhood trauma showed the highest epigenetic aging compared with control subjects without childhood trauma (p=0.001, Cohen’s d=0.29), highlighting the fact that this major depression and childhood trauma subgroup is associated with the highest epigenetic aging (see Figure S1 in the data supplement). Notably, greater severity of symptoms of chronic major depression was correlated with childhood trauma (r=0.39, p<0.001), making it difficult to discern which of these two factors drives increased epigenetic aging. Linear regression indeed showed that both childhood trauma (β=0.08, p=0.01) and depression severity score (β=0.07, p=0.03) were significant predictors of epigenetic aging when analyzed in the same model.

The mean age of the postmortem brain samples was 55.2 years (SD=19.3, range=20–100), and 45.4% of the sample were female. Groups were matched on age and sex. The mean postmortem interval was 35.1 hours (SD=21.1), and the mean pH was 6.51 (SD=0.25) across samples. Epigenetic aging was uncorrelated with pH or postmortem interval. Table 4 summarizes the descriptive characteristics by brain collection. Only the control (N=67) and depression (N=74) samples from the same brain collection were included in analyses (see the data supplement for more detail).

Our replication findings in independent brain samples supported our findings in NESDA and again showed that the epigenetic aging was higher in the depression group than in the control group (b=1.11, χ²=3.41, p=0.03). The beta indicates that the depression group was estimated to be on average 1.11 years older than the control group. The phenotype information available from the postmortem samples was limited, and therefore we were unable to attempt any replication of the exploratory clinical associations observed, such as childhood trauma, in NESDA.

When evaluating the overlap between both epigenetic aging indicators, we found that after correcting for multiple testing, the top 1% findings from the epigenetic aging MWAS in blood were significantly enriched for CpGs in the top 0.5% of the epigenetic aging MWAS from brain (odds ratio=1.19, p<0.001). To examine possible processes underlying epigenetic aging in both tissues, we performed pathway analyses on the 1,084 overlapping CpGs associated with epigenetic aging, leading to 330 genes (90.7%) that were present in at least one GO category. Subsequently, this resulted in 53 significantly enriched GO terms (see Table S4 in the data supplement). The top GO terms included neurogenesis (p=9.79×10⁻⁹), neuron differentiation (p=5.34×10⁻⁸), and regulation of neuron death (p=4.67×10⁻⁵), indicating that several depression-relevant pathways were enriched in the cross-tissue epigenetic aging indicators.