Non-linear patterns in age-related DNA methylation may reflect CD4⁺ T cell differentiation

ABSTRACT

DNA methylation (DNAm) is an important epigenetic process involved in the regulation of gene expression. While many studies have identified thousands of loci associated with age, few have differentiated between linear and non-linear DNAm trends with age. Non-linear trends could indicate early- or late-life gene regulatory processes. Using data from the Illumina 450K array on 336 human peripheral blood samples, we identified 21 CpG sites that associated with age (P<1.03E-7) and exhibited changing rates of DNAm change with age (P<1.94E-6). For 2 of these CpG sites (cg07955995 and cg22285878), DNAm increased with age at an increasing rate, indicating that differential DNAm was greatest among elderly individuals. We observed significant replication for both CpG sites (P<5.0E-8) in a second set of peripheral blood samples. In 8 of 9 additional data sets comprising samples of monocytes, T cell subtypes, and brain tissue, we observed a pattern directionally consistent with DNAm increasing with age at an increasing rate, which was nominally significant in the 3 largest data sets (4.3E-15<P<0.039). cg07955995 and cg22285878 reside in the promoter region of KLF14, which encodes a protein involved in immune cell differentiation via the repression of FOXP3. These findings may suggest a possible role for cg07955995 and cg22285878 in immunosenescence.

KEYWORDS:

• Aging
• DNA methylation
• FOXP3, immunosenescence
• KLF14

Introduction

Gene expression is regulated via highly coordinated epigenetic changes, of which the most studied is DNA methylation (DNAm). DNAm is the binding of a methyl group (−CH₃) to DNA, which, in mammals, occurs most commonly at a cytosine nucleotide that resides 5′ to a guanine nucleotide, referred to as a CpG site.¹ Clusters of CpG-rich regions, known as CpG islands (CGIs), are often found in gene promoters and are typically hypomethylated.² The position of the CpG site within the transcription unit has a substantial influence on how DNAm impacts expression of the downstream gene,³ and evidence suggests that DNAm in gene promoters can result in gene silencing.⁴

While patterns of DNAm across the genome vary according to tissue type and environmental exposure, many studies have shown that age explains a substantial portion of the variation in human DNAm.^5-13 DNAm has been observed to decrease genome-wide with age,¹⁴ although at a regional level the relationship between DNAm and age is much more nuanced. For example, CGIs are associated with an increasing rate of DNAm with age, whereas non-CGIs are associated with a decreasing rate of DNAm with age.⁷ Furthermore, postnatal DNA is hypomethylated and undergoes a rapid increase in DNAm in early life before stabilizing in adulthood, followed by a gradual decrease later in life.¹⁵

While linear trends can provide a useful summary of age-related DNAm patterns, deviations from linearity can provide insight into the relationship between DNAm and senescence. Heretofore, studies investigating non-linear patterns in age-related DNAm have uncovered patterns of rapid DNAm changes in early life, followed by stabilization in later life. For example, differences in the rates of DNAm have been observed between pediatric and adult subjects,¹³ with many CpG sites showing a decreasing rate of change with age. Such patterns may reflect early-life developmental processes. Perhaps more pertinent to the study of senescence and age-related diseases is an investigation of the reverse pattern: stable DNAm levels in early life followed by rapid methylation/demethylation in late life. To our knowledge, this work is the first to investigate the reverse pattern. Here, we analyze genome-wide DNAm data in a set of peripheral blood samples with the goal of identifying loci for which DNAm increases or decreases at an increasing rate with age. We follow up this analysis in 10 additional data sets comprising distinct tissue types including peripheral blood, dorsolateral prefrontal cortex (dlPFC) tissue, neurons isolated from PFC tissue, glial cells isolated from PFC tissue, and purified blood cell types isolated from peripheral blood, including monocytes, CD4⁺ T cells, and CD8⁺ T cells.

Results

Using DNAm data obtained from peripheral blood samples from 336 individuals aged 15–78 from the Grady Trauma Project (GTP, Table 1), we performed separate linear regressions for each of 483,399 CpG sites to model DNAm as a function of age and other covariates, including estimated cell type proportions (Methods). A total of 25,723 CpG sites were significantly associated with age (P<1.03E-7). Subsequently, we tested for non-linear trends in these 25,723 CpG sites by fitting quadrilinear models that included a term for age². For the subsequent analysis, 21 CpG sites had a Bonferroni-significant age² term (P<1.94E-6=0.05/25,723). Among these 21 CpG sites (Table S1 and Figure S1), we observed 3 general patterns:

Table 1. Reference information regarding each of the 11 data sets analyzed.

Gene series	Study	Abbrev.	Tissue type	Publication
GSE72680	Grady Trauma Project	GTP	Peripheral blood	Zannas, Arloth, Carrillo-Roa, Iurato, Röh, Ressler, Nemeroff, Smith, Bradley, Heim, Menke, Lange, Brückl, Ising, Wray, Erhardt, Binder and Mehta^53
GSE60132	TOPS Family Study of Epigenetics	TOPS	Peripheral blood	Ali, Cerjak, Kent, James, Blangero, Carless and Zhang^54
GSE56581	MESA Epigenomics and Transcriptomics Study (human T cells)	MESA-T	Purified CD4^+ T cells	Reynolds, Taylor, Ding, Lohman, Johnson, Siscovick, Burke, Post, Shea, Jacobs, Stunnenberg, Kritchevsky, Hoeschele, McCall, Herrington, Tracy and Liu^56
GSE56046	MESA Epigenomics and Transcriptomics Study (human monocytes)	MESA-M	Purified monocytes	Reynolds, Taylor, Ding, Lohman, Johnson, Siscovick, Burke, Post, Shea, Jacobs, Stunnenberg, Kritchevsky, Hoeschele, McCall, Herrington, Tracy and Liu^56
GSE59065	Estonian Genome Center Investigation of Age-related epigenetics and immune system function in PBL, CD4+ and CD8+ T cells	EGC-PBL	Peripheral blood leukocytes (PBL), CD4+ and CD8+ T cells	Tserel, Kolde, Limbach, Tretyakov, Kasela, Kisand, Saare, Vilo, Metspalu, Milani and Peterson^40
		EGC-CD4
		EGC-CD8
phs000741.v1.p1	Genetics Of Lipid Lowering Drugs and Diet Network	GOLDN	Purified CD4^+ T cells	Irvin, Zhi, Joehanes, Mendelson, Aslibekyan, Claas, Thibeault, Patel, Day, Jones, Liang, Chen, Yao, Tiwari, Ordovas, Levy, Absher and Arnett^63
GSE74193	Schizophrenia-related DNAm and gene expression	CTX	Dorsolateral prefrontal cortex (dlPFC)	Jaffe, Gao, Deep-Soboslay, Tao, Hyde, Weinberger and Kleinman^57
GSE41826	National Institute of Child Health and Human Development (NICHD) Brain Bank of Developmental Disorders – Cell Heterogeneity	NICHD-G	Prefrontal cortex (PFC) Glial cells	Guintivano, Aryee and Kaminsky^58
		NICHD-N	Prefrontal cortex neurons

Pattern 1

Two of the 21 CpG sites (cg07955995 and cg22285878 in KLF14) exhibited a low level of DNAm (≤ 0.04) and a low rate of change (slope or ΔDNAm proportion per year = −8.7E-4 and −5.7E-4) at the minimum age of the GTP age range (age_min= 15.9 years) that positively accelerated as age increased (slope at age 70 = 3.1E-3 and 1.7E-3; first row of Tables 2 and 3).

Table 2. Statistics corresponding to the age² term of the regression model fitted separately to each of the 11 data sets for the CpG site, cg07955995 (Equation 1(Equation 1) ${\beta }_i={\mathit{\alpha }}_0+{\mathit{\gamma }}_{\mathit{a}\mathit{g}\mathit{e}}\left(AG{E}_i\right)+{\mathbf{\gamma }}_{\mathit{a}\mathit{g}{\mathit{e}}^2}\left(AG{E}_i^2\right)+{\displaystyle \mathbf{\sum }}_{j = 1}^k{\mathit{\gamma }}_{\mathit{j}}co{v}_{ij}+{\mu }_i+{ϵ}_i$(Equation 1) ). The regression model includes covariates. Columns 1–3 include abbreviations of the data set, cell/tissue type, and the sample size. Columns 4–6 include the estimated rate of change in DNAm proportion per year at 50, 60, and 70 years of age. Columns 7–10 include the coefficient, standard error, T-statistic, and P-value on the age² term.

Dataset	Cell/tissue type	n	${\gamma }_{age = 50}$	${\gamma }_{age = 60}$	${\gamma }_{age = 70}$	${\gamma }_{ag{e}^2}$	$S{E}_{ag{e}^2}$	$T_{ag{e}^2}$	$P_{ag{e}^2}$
GTP	Whole blood	336	1.6E-03	2.4E-03	3.1E-03	3.7E-05	4.1E-06	9.1	1.6E-17
TOPS	Whole blood	192	2.0E-03	2.7E-03	3.4E-03	3.4E-05	3.9E-06	8.6	4.3E-15
EGC-PBL	Whole blood	97	1.3E-03	1.1E-03	8.1E-04	−1.3E-05	2.3E-05	−0.57	0.57
MESA-M	Monocyte	1202	1.1E-03	1.6E-03	2.0E-03	2.2E-05	1.0E-05	2.2	0.030
EGC-CD8	CD8^+ T cell	100	1.2E-03	1.8E-03	2.4E-03	2.8E-05	2.6E-05	1.1	0.29
MESA-T	CD4^+ T cell	214	7.7E-04	1.2E-03	1.6E-03	2.0E-05	1.6E-05	1.2	0.22
EGC-CD4	CD4^+ T cell	99	9.5E-04	1.0E-03	1.1E-03	3.6E-06	1.7E-05	0.21	0.83
GOLDN	CD4^+ T cell	991	4.3E-03	4.9E-03	5.5E-03	4.0E-06	7.7E-07	5.2	2.0E-07
CTX	dlPFC	346	1.7E-03	1.8E-03	1.9E-03	5.3E-06	2.6E-06	2.1	0.038
NICHD-G	Glial cell (FC)	50	1.3E-03	1.3E-03	1.2E-03	−2.9E-06	8.7E-06	−0.33	0.74
NICHD-N	Neuron (FC)	50	2.0E-03	2.3E-03	2.5E-03	1.3E-05	1.0E-05	1.2	0.23

The regression model includes covariates. Columns 1–3 include abbreviations of the data set, cell/tissue type, and the sample size. Columns 4–6 include the estimated rate of change in DNAm proportion per year at 50, 60, and 70 years of age. Columns 7–10 include the coefficient, standard error, T-statistic, and P-value on the age² term.

Table 3. Statistics corresponding to the age² term of the regression model fitted separately to each of the 11 data sets for the CpG site, cg22285878 (Equation 1(Equation 1) ${\beta }_i={\mathit{\alpha }}_0+{\mathit{\gamma }}_{\mathit{a}\mathit{g}\mathit{e}}\left(AG{E}_i\right)+{\mathbf{\gamma }}_{\mathit{a}\mathit{g}{\mathit{e}}^2}\left(AG{E}_i^2\right)+{\displaystyle \mathbf{\sum }}_{j = 1}^k{\mathit{\gamma }}_{\mathit{j}}co{v}_{ij}+{\mu }_i+{ϵ}_i$(Equation 1) ).

Dataset	Cell/tissue type	n	${\gamma }_{age = 50}$	${\gamma }_{age = 60}$	${\gamma }_{age = 70}$	${\gamma }_{ag{e}^2}$	$S{E}_{ag{e}^2}$	$T_{ag{e}^2}$	$P_{ag{e}^2}$
GTP	Whole blood	336	8.8E-04	1.3E-03	1.7E-03	2.1E-05	3.3E-06	6.4	5.3E-10
TOPS	Whole blood	192	9.8E-04	1.3E-03	1.6E-03	1.5E-05	2.7E-06	5.7	5.0E-08
EGC-PBL	Whole blood	97	6.0E-04	7.7E-04	9.5E-04	8.6E-06	1.4E-05	0.62	0.54
MESA-M	Monocytes	1202	5.4E-04	8.7E-04	1.2E-03	1.6E-05	6.3E-06	2.6	0.010
EGC-CD8	CD8^+ T cells	100	6.3E-04	1.4E-03	2.1E-03	3.7E-05	1.7E-05	2.2	0.033
MESA-T	CD4^+ T cells	214	6.2E-04	7.7E-04	9.2E-04	7.5E-06	1.2E-05	0.64	0.52
EGC-CD4	CD4^+ T cells	99	7.8E-04	8.2E-04	8.5E-04	1.7E-06	2.1E-05	0.082	0.93
GOLDN	CD4^+ T cells	991	1.2E-03	1.3E-03	1.3E-03	3.9E-06	8.5E-07	4.6	4.2E-06
CTX	dlPFC	346	1.3E-03	1.5E-03	1.6E-03	6.9E-06	2.2E-06	3.1	0.0018
NICHD-G	Glial cell (FC)	50	1.5E-03	1.5E-03	1.4E-03	−1.9E-06	8.3E-06	−0.22	0.82
NICHD-N	Neuron (FC)	50	1.1E-03	1.2E-03	1.4E-03	7.9E-06	8.3E-06	0.95	0.35

The regression model includes covariates. Columns 1–3 include abbreviations of the data set, cell/tissue type, and the sample size. Columns 4–6 include the estimated rate of change in DNAm proportion per year at 50, 60, and 70 years of age. Columns 7–10 include the coefficient, standard error, T-statistic, and P-value on the age² term.

Pattern 2

One CpG site (cg14293223 in GFI1B) exhibited high DNAm at age_min and a small slope that became increasingly negative with age (Table S1, Figure S1). It appeared that the oldest individual, who also had the lowest level of DNAm, was driving the quadrilinear effect between DNAm and age: when the oldest individual was removed from the model, both the age term and the age² term were no longer significant.

Pattern 3

The remaining 18 CpG sites that were modeled exhibited a positive relationship between DNAm and age at age_min, whose slope decreased with age, becoming negative among older individuals (Table S1, Figure S1).

Patterns 1 and 2 are of particular interest because they indicate stable DNAm in early life, and a more rapidly changing rate of DNAm late in life. To follow up on the 2 CpG sites (cg07955995 and cg22285878) exhibiting Pattern 1 and the one CpG site (cg14293223) exhibiting Pattern 2 in the GTP, we attempted replication in a second set of samples derived from peripheral whole blood from the Take Off Pounds Sensibly (TOPS) Family Study of Epigenetics, whose subjects were members of the TOPS Club. We observed strong replication at cg07955995 (P = 4.3E-15) and cg22285878 (P = 5.0E-8), with extremely similar effect sizes in TOPS and GTP (${\gamma }_{ag{e}^2}$ in Tables 2–3). We did not find a significant quadrilinear relationship between DNAm and age at the CpG site exhibiting Pattern 2 (cg14111928; P = 0.16).

To investigate whether the quadrilinear patterns observed for cg07955995 and cg22285878 were specific to whole blood, we fit a similar model in 9 additional data sets that were derived from whole blood, subtypes of whole blood cells (monocytes, CD8⁺ T cells, and CD4⁺ T cells), and brain [tissue samples from dorsolateral prefrontal cortex (dlPFC) and samples of glial cells and neurons, both isolated from the prefrontal cortex tissue] (Table 1). All data sets exhibited increasing DNAm with age at cg07955995 and cg22285878 (Table S2 and S3). In the quadrilinear model, the increasing rate of increase in DNAm with age that was observed in the GTP and TOPS analyses was observed in 7 of 9 data sets for cg07955995 and 8 of 9 data sets for cg22285878 (Tables 2 and 3, Fig. 1). The age² terms were significant in MESA-M (samples of isolated monocytes), GOLDN (CD4⁺ T-cells), and CTX (dlPFC) for both CpG sites (4.3E-15<P<0.039; Tables 2 and 3). The coefficient on age² had a positive sign with comparable effect sizes for nearly all data sets, with the exceptions of EGC-PBL at cg07955995 and NICHD-G at both cg07955995 and cg22285878. These exceptions could potentially be attributable to a lack of older individuals in NICHD-G, and a lack of individuals of ages 23–72 years in the EGC study, given the difficulty of identifying non-linear patterns using just 2 extreme endpoints. However, the pattern is consistent overall and the comparable effect sizes suggest that the main difference between data sets with a significant vs. non-significant age² term is statistical power due to differences in both age range and sample size.

Figure 1. Scatterplots of DNAm (y-axis) vs. age (x-axis) with fitted quadratic regression lines for each of the 11 data sets, holding covariates constant at their mean values.

These results show a consistent pattern in age-related DNAm at these 2 CpG sites across white blood cell subtypes. In particular, MESA-M (samples of isolated monocytes) and GOLDN (samples of isolated CD4⁺ T cells) exhibited evidence of an age-quadratic relationship with DNAm. However, it is well-known that peripheral blood is a heterogeneous tissue, and cell type proportions may vary across samples.¹⁶ Analyses of peripheral blood were adjusted for cell type proportion (see Methods), but to further investigate the possible role of cell type in the observed relationships, we examined correlations between age and estimated cell type proportions in the GTP. Of the 6 cell types considered, only CD8⁺ T cells showed a small but significant correlation with age (r=−0.188, P=0.005; Figure S2). Notably, the 2 oldest individuals (who are also the 2 individuals with the highest DNAm) had similar proportions of CD8⁺ T cells compared with the remaining individuals (Figure S2). For MESA-M, we investigated whether age associated with the cell type contamination proportions estimated by Reynolds et al. (Figure S3).⁵⁶ We found small but significant correlations with age with proportions of B cells (r = 0.0584, P=0.0428), natural killer cells (r = 0.1034, P = 0.0003), and T cells (r=−0.0703, P=0.0148). We also plotted cell type contamination proportion against DNAm at both cg07955995 and cg22285878 (Figures S4 and S5). At cg07955995, we observed a small but significant correlation between DNAm and estimated B cell contamination (r=−0.0906, P = 0.0017). At cg22285878, we observed small but significant correlations with estimated proportions of B cells (r=−0.0689, P = 0.0169) and natural killer cells (r = 0.0645, P=0.0254).

Because an increased variation in DNAm among older subjects appeared to be driving the quadrilinear relationship between DNAm and age in GTP, TOPS, MESA-M, and GOLDN, we investigated whether the variance in DNAm differed significantly between younger and older individuals in each data set. The EGC data set consisted solely of 2 age groups: a young age group (≤ 34 years) and an old age group (≥ 73 years). For comparability, we split the GTP, TOPS, CTX, and GOLDN data sets into young (≤ 34 years) and old ≥ 73 years) age groups. MESA-M and MESA-T have age ranges of 45–79 and 44–83, so we could not split these data sets in the same age groups; similarly, NICHD-G and NICHD-N only included 2 individuals ≥ 73 years. For these data sets, we instead used the median values (58 for MESA-M, 60 for MESA-T, and 23 for both NICHD-G and NICHD-N) to split the data sets into younger (≤median) and older (>median) age groups. We then computed the variance of DNAm at cg07955995 and cg22285878 separately for the 2 age groups. Fig. 2A and 2B plot the deviation of these variances from the $y = x$ line for all data sets. When comparing all data sets GTP and TOPS had the highest variances among the old age groups at cg07955995. GTP also had the highest variance among the old age group at cg22285878, followed by CTX. When comparing EGC-CD4, EGC-CD8, and EGC-PBL (derived from the same set of subjects) EGC-CD8 exhibited the highest variance in the old age group at cg07955995 while EGC-CD4 exhibited the lowest variance. At cg22285878, EGC-CD4 exhibited the highest variance of the 3 in the old age group whereas PBL exhibited the lowest variance. Across all data sets, variance was significantly greater in the old age groups compared with the young age groups (Fig. 2A and 2B; Table 4), suggesting that the quadrilinear pattern is in part driven by increased variability in DNAm across older subjects.

Figure 2. Plot of the variance of methylation across subjects at cg07955995 (left) and cg22285878 (right) for the older age group (73 years or older) against the younger age group for each data set. MESA-M, MESA-T, NICHD-G, and NICHD-N appear with different symbols because the age groups were calculated differently than the other data sets. For these 4 data sets, the median ages (in years) were used (58 for MESA-M, 60 for MESA-T, and 23 for both NICHD-G and NICHD-N) to create young groups (≤median) and old groups (>median) because MESA-M and MESA-T had no individuals less than 34 years of age and NICHD-G and NICHD-N only had 2 individuals older than 73 years.

Table 4. Variances of the young age group, variances of the old age group, and P-values corresponding to the F-statistic calculated as the ratio of variances of the young and old age groups for each of the 11 data sets.

	$cg07955995$			$cg22285878$
Dataset	$Va{r}_{young}$	$Va{r}_{old}$	$P_{Var}$	$Va{r}_{young}$	$Va{r}_{old}$	$P_{Var}$
GTP	8.3E-05	2.7E-03	1.4E-11	6.6E-05	2.0E-03	5.6E-11
TOPS	6.5E-05	3.0E-03	2.8E-26	6.6E-05	6.0E-04	2.0E-08
EGC-PBL	1.8E-04	1.5E-03	3.8E-05	1.5E-04	4.3E-04	9.7E-03
MESA-M	4.2E-04	1.6E-03	6.4E-56	1.8E-04	7.3E-04	1.9E-62
EGC-CD8	3.9E-04	1.9E-03	1.0E-11	3.6E-04	8.0E-04	1.4E-04
MESA-T	1.9E-04	4.2E-04	2.1E-03	1.4E-04	2.1E-04	1.0E-09
EGC-CD4	2.9E-04	6.6E-04	4.4E-08	1.9E-04	1.2E-03	3.0E-03
GOLDN	7.5E-05	3.8E-04	2.1E-03	7.8E-05	2.7E-04	1.1E-03
CTX	4.3E-04	1.2E-03	6.7E-19	4.5E-04	1.4E-03	9.8E-12
NICHD-G	3.4E-04	9.2E-05	0.76	3.0E-04	1.9E-04	0.54
NICHD-N	4.4E-04	1.2E-07	0.9997	4.3E-04	1.9E-04	0.65

Columns 2–4 correspond to cg07955995 and columns 5–7 correspond to cg22285878.

Although many of the studies described above had limited numbers of older individuals, this pattern was also observed in an additional blood-based data set (peripheral blood mononuclear cells) that included 122 nonagenarians and 21 younger controls aged 19–30 from the Vitality 90+ Study (Gene Expression Omnibus accession number GSE58888);¹⁷ this data set was not analyzed with the others because individual ages were unavailable. Figure S6 shows that for both CpG sites, the nonagenarians in the study had greater mean DNAm levels (P<6E-13) and variance of DNAm (P=8.7e-6 for cg07955995 and 0.0014 for cg22285878) than the younger controls, supporting the results observed in the other data sets.

To test the robustness of the quadrilinear model and ensure we did not miss other interesting non-linear patterns, we also performed a secondary analysis in which we fit spline regressions for the 25,723 age-associated CpG sites in GTP (see Methods for details). Results were similar to those of the quadrilinear model: 20 CpG sites were significant for the spline-based model, whereas 21 CpG sites were significant for the quadrilinear model. Among these, 15 were significant for both models (Figure S7).

As technical validation of the array-based results for cg07955995 and cg22285878, we used the Sequenom EpiTyper technology to interrogate a section of the promoter region that included cg07955995 and cg22285878 for the 2 oldest and 2 of the youngest GTP individuals in the GTP data set (subject ages = 21, 21, 75, and 78). Standard curve analysis based on commercially available methylation standards (Figure S8) suggested that the assay underestimates DNAm at these 2 sites. For cg07955995, DNAm levels from the targeted assay of young vs. old subjects were consistent with the original array-based results, though the difference between older and younger subjects was somewhat compressed (Figure S9). For cg22285878, the EpiTyper DNAm levels show little difference between the older and younger subjects; this may be due to a combination of the smaller differences observed in GTP and other data sets for cg22285878 and the underestimation of DNAm indicated in Figure S8. Notably, increased DNAm for the older subjects was observed at many of the 22 CpG sites in the region, displaying more striking differences at some of the 20 flanking CpG sites than at cg07955995 and cg22285878 (Figure S10).

The 2 CpG sites (cg07955995 and cg22285878) exhibiting Pattern 1 are located 14 base pairs apart in the promoter region of the gene KLF14, and show high pairwise correlation in all data sets (0.70 < r < 0.87; Figure S11). KLF14 is an imprinted, maternally expressed gene with a hypomethylated CpG island.¹⁸ Previous work has found associations between DNAm at KLF14 and age in pancreatic islets, adipose tissue, and whole blood.^19-22 In addition, KLF14 has been suggested to be involved in the regulation of inflammation.²³ To investigate a possible relationship between DNAm and inflammation, we compared interleukin 6 (IL-6) and C-reactive protein (CRP) across GTP subjects (Figure S12 and S13). The 2 oldest individuals, who also exhibited the highest β-values at both cg07955995 and cg22285878, did not have plasma levels of inflammatory markers that differed from the rest of the GTP subjects. There was also no association between levels of either inflammatory marker and DNAm levels at cg07955995 or cg22285878 (0.4972 ≤ P ≤ 0.9780). We also investigated CRP and IL-6 in GOLDN, along with 3 additional inflammatory markers (TNFα, MCP1, and sIL2Rα). No significant associations were found between DNAm and these 5 markers at either CpG site.

Discussion

Many studies have investigated the relationship between DNAm and human age. While significantly different rates of DNAm have been observed between pediatric and adult populations,¹³ we are unaware of studies investigating changing rates of DNAm with age in the same cohort. Such investigations are important to unravel the complex relationship between DNAm and senescence. In particular, CpG sites demonstrating an increasing rate of DNAm change with age may be involved in processes relevant to aging. To understand the importance of these results, we discuss the tissue and cell type specific patterns we observe, and we outline a possible pathway whereby suppression of KLF14 via DNAm at cg07955995 and cg22285878 may play a role in immunosenescence. In light of previous findings linking KLF14 to metabolic outcomes, we discuss whether our results are consistent with aging-related onset of metabolic disorders, immunosenescence, or both. Finally, we discuss limitations of our study and implications for future research.

While an increased rate of change in DNAm with age suggests relevance to aging-specific processes, it remains unclear what biologic phenomenon might be responsible for such age-related DNAm. Many studies have found not only tissue-specific but also cell type-specific DNAm patterns.^24,25 Furthermore, such tissue- and cell type-specific DNAm patterns could be involved in cell lineage differentiation.¹⁶ We investigated age-related DNAm patterns among data sets of different tissue and cell types to explore the possibility that the changing rate of DNAm that we observed was specific to a tissue, a cell type, or group of cell types. We found that CpG sites cg07955995 and cg22285878 had DNAm that increased at an increasing rate with age across nearly all of the 11 data sets we investigated. While our investigation is not exhaustive, the patterns so far suggest that our findings generalize across several sub-types of whole blood. It remains to be seen whether this pattern will generalize to other tissues. In PFC tissue (CTX, NICHD-N), the 2 CpG sites exhibited a pattern consistent with the analyses in blood (an increasing rate of increase and larger variance in older subjects). While this pattern was significant in dlPFC tissue (CTX), it was insignificant in PFC neurons (NICHD-N), and was not observed in PFC glial cells (NICHD-G). Future studies should investigate these patterns in additional tissues as additional data sets become available.

The analyses of peripheral whole-blood data (GTP and TOPS) revealed a significant quadrilinear relationship between β-values and age that was consistent with an increasing rate of DNAm change with age. In subsequent analyses, we investigated whether this finding replicated in blood cell subtypes (Tables 2 and 3). Isolated monocytes (MESA-M) also showed a significant quadrilinear relationship between β and age at both cg07955995 and cg22285878. In isolated CD8⁺ T cells (EGC-CD8), the relationship was not significant, but the coefficient for the age² term was comparable to those observed in monocytes, suggesting that the age range (no individuals between the ages of 23 and 72 years) and small sample sizes (n = 100) hindered our ability to detect a quadrilinear association. However, the older subjects exhibited significantly greater variance than the younger subjects, consistent with what was observed for other cell subtypes. In isolated CD4⁺ T cells (EGC-CD4, MESA-T, and GOLDN), we observed quadrilinear associations the same direction with one exception (EGC-CD4 for cg22285878) but the associations were significant in GOLDN only. Like EGC-CD8, EGC-CD4 had a small sample size (n = 99) and was limited to individuals at the endpoints of the age distribution; MESA-T also had a small sample size (n = 214). Effect sizes were comparable across the 3 CD4⁺ T cell data sets: the coefficients on the quadratic terms were actually larger in MESA-T than in GOLDN for both CpG sites, and the coefficients for EGC-CD4 and GOLDN were similar in magnitude for cg07955995.

Given the patterns observed in T cell subtypes, it is notable that cg07955995 and cg22285878 reside in the promoter region of a gene called KLF14, which is involved in immune cell differentiation, particularly the conversion of CD4⁺CD25⁻ naïve T cells to CD4⁺CD25⁺ regulatory T cells. KLF14 belongs to a family of 17 proteins known as Kruppel-like factor (KLF) proteins, which are involved in immune cell differentiation. KLF proteins exert their effect via the binding to gene promoters and enhancers.²³ A recent study found that KLF14 represses FOXP3 via epigenetic regulation at the Treg-specific demethylated region (TSDR) of FOXP3 in vitro and in vivo using a mouse model.²³ Previous work has shown that FOXP3 induces and is required for the conversion of CD4⁺CD25⁻ naïve T cells to CD4⁺CD25⁺ regulatory T cells via the joint stimulation of TCR and TGF-β.^26-28 A study of CD25-deficient mice found that the overexpression of proinflammatory cytokines, including IL-2 and IFN-γ, in response to bacterial superantigen stimulation, particularly staphylococcal enterotoxin B, was curbed by injection of CD4⁺CD25⁺, suggesting that CD4⁺CD25⁺ plays a role in regulating antigen-induced inflammatory responses.²⁹ Moreover, evidence suggests a relationship between immunosenescence (the general deterioration of the immune system over the course of the lifespan) and CD4⁺CD25⁺ regulatory T cell levels and function.^30,31 In consideration of the aforementioned studies and our findings, we outline one possible biologic pathway whereby late-life DNAm at cg07955995 and cg22285878 in CD4⁺ T cells could be a response to immunosenescence (Figure S14): (1) Age-related DNAm in the promoter region of KLF14 downregulates expression of KLF14; (2) In the absence of normal levels of KLF14, which otherwise inhibits FOXP3, FOXP3 expression is upregulated; (3) Higher levels of FOXP3 induces a greater number of naïve CD4⁺25⁻ T cells to convert to CD4⁺25⁺ T cells (sometimes referred to as Treg cells); and (4) CD4⁺25⁺ T cells regulate the heightened inflammation observed in immunosenescent individuals. To investigate the conjectured pathway outlined above, we assessed correlation between DNAm at cg07955995 and cg22285878 and levels of inflammatory markers in GTP and GOLDN individuals. However, there was no association between inflammatory markers and DNAm at cg07955995 and cg22285878 in either GTP (Figures S12 and S13) or GOLDN, and the 2 oldest individuals in the GTP data set, who also exhibited the highest β-values at both cg07955995 and cg22285878, were well within the normal range for plasma levels of IL-6 and CRP. Nevertheless, Pontoux et al. investigated levels of IL-2 and IFN-γ in response to injection of CD4⁺CD25⁺.²⁹ Thus, it is possible that we were unable to find any indication of immunosenescence because this particular pathway affects inflammatory markers that we did not investigate. Future studies should collect additional inflammatory markers, including IL-2 and IFN-γ, on a larger subset of elderly individuals to further investigate this putative pathway and its relationship to immunosenescence.

While we suggest that DNAm at cg22285878 and cg07955995 could regulate KLF14 expression and possibly contribute to an immunosenescent phenotype in older adults, we did not investigate expression of KLF14 because it did not show sufficient expression to pass standard quality control procedures in the available blood-based expression data sets (GTP, MESA-M, MESA-T, and a small subset of GOLDN). Because KLF14 did not show expression in these data sets, it is unlikely KLF14 is expressed over long periods of time. Data from the Common Fund (CF) Genotype-Tissue Expression Project (GTEx) indicates that KLF14 expression in whole blood exhibits a very low mean gene expression level with a large upward skew relative to other tissue types (Figure S15), which is what one would expect if a gene is transiently expressed in a tissue.³² Many genes have been found to undergo transient stints of expression.^33-35 In fact, one study found transient expression of cell surface proteins in monocytes in response to inflammatory cytokines.³⁶ If KLF14 regulates T cell conversion via intermittent expression followed by long periods of inactivity, then it is unlikely that we would be able to detect KLF14 gene expression in our samples. If this is the case, it is possible that KLF14 is activated in the blood to decrease excessive levels of CD4⁺CD25⁺ regulatory T cells. If heightened DNAm at cg22285878 and cg07955995 in immunosenescent individuals regulates KLF14 expression by blocking or partially blocking brief and intermittent periods of KLF14 expression in some CD4⁺ T cells, then we would expect Treg cells to accumulate in the blood of immunosenescent individuals. Interestingly, previous research has found that Treg cells accumulate with age.^30,31 Another possible explanation for lack of detectable levels of KLF14 gene expression is a technical failure of the Illumina probe, which has been suggested in a previous study.⁴⁴

While aging affects all aspects of the immune system, T cells are the most severely affected.³⁷ While many factors may contribute to an immunosenescent phenotype, one such factor, namely, cytomegalovirus (CMV), infects an estimated 50–80% of North Americans by age 40 and is most commonly symptomatic among immunodeficient individuals, such as the elderly.³⁸ CMV has been observed to induce CD8⁺ T cell proliferation,³⁹ which has been hypothesized to play a role in age-related DNAm change.⁴⁰ Moreover, a substantial portion of peripheral blood CD4⁺ T cell and CD8⁺ T cell responses has been observed in CMV-seropositive individuals comprising approximately 10% of T cell memory compartments.⁴¹ Our investigation found a significant quadrilinear effect between DNAm and age at cg07955995 and cg22285878 in peripheral whole-blood (GTP and TOPS), isolated CD4⁺ T cells, and, to a lesser (but still significant) extent, in isolated monocytes (MESA-M). These findings suggest that CD4⁺ T cells could be partially driving the quadrilinear effect, which is consistent with the putative pathway outlined above. In GTP (at cg07955995 and cg22285878) and TOPS (at cg07955995), the oldest individuals had the highest DNAm levels, and in several of the remaining data sets (MESA-M, EGC-PBL, EGC-CD8, EGC-CD4; Fig. 1), we observed the highest DNAm levels among a subset of the older individuals, which could indicate a possible immunosenescent phenotype in these individuals.

In addition to the possible role for KLF14 in immune cell differentiation and inflammation, previous studies have linked it to metabolic outcomes. Genetic variants in the promoter region of KLF14 associate robustly with risk for type II diabetes (T2D) and high-density lipoprotein cholesterol.^42,43 It has also been reported that DNAm of cg22285878 in pancreatic islets associates both with aging and with measures of insulin secretion independent of age.¹⁹ GTEx data indicate that KLF14 is expressed in the adrenal gland as well as in both visceral and subcutaneous adipose tissue (Figure S15), and it has been reported that KLF14 acts as a master regulator in adipose tissue, regulating a set of genes correlated with metabolic traits.⁴⁴ A growing body of literature suggests that inflammation is involved in insulin resistance and chronic metabolic disorders, such as T2D,^45-49 so it is possible that the roles of KLF14 in inflammation and in chronic metabolic disorders represent 2 sides of the same coin. The high levels of DNAm at cg07955995 and cg22285878 that we observed among elderly individuals may reflect immunosenescence, onset of a chronic metabolic disorder, or both.

To further our understanding of the role of age-related DNAm at cg07955995 and cg22285878 and its possible relationship with metabolic disorders and immunosenescence, future research should focus on cell type specific data with large sample size and a wide age range. Moreover, our study suggests a possible role of T cell subtypes, specifically, Treg cells (CD4⁺CD25⁺ T cells), which are difficult to isolated from other CD4⁺CD25⁺ T cells. As we overcome challenges involving the isolation of cell subtypes, particularly T cell subtypes, future studies should investigate DNAm of T cell subtypes to further investigate the biologic pathway conjectured here. In addition, research should pay particular attention to DNAm data on T2D, obese, and immunosenescent individuals, including those who are CMV-seropositive.

Several limitations of this study are worth noting. First, all 11 data sets are cross-sectional, so it is not possible to infer that an increase in age yields an increase in DNAm, merely that there is an association. Ideally, longitudinal data could better assess the relationship in question. Second, the primary analysis was conducted on a data set with a small sample size (n = 336) and few older individuals, as were some of the follow-up analyses (NICHD-N and NICHD-G). These data sets may have lacked the statistical power to detect all CpG sites that exhibit a non-linear trend. A related concern is that the small number of older individuals in these samples could potentially exert undue influence as outliers, which could lead to false positive results. However, the consistency of our initial result in 9 of 11 data sets derived from diverse tissues and the corroboration of observed patterns in a larger data set of nonagenarians suggest that the patterns observed for these 2 CpG sites reflect replicable and generalizable biologic differences and are robust rather than outlier-driven. Third, it is possible that unobserved environmental factors, such as pollutants, accumulate in the body over time and help drive the age-related relationship at these 2 CpG sites. Finally, we were unable to detect KLF14 expression in the available data sets, which limited the possible follow-up analyses. While it is possible that KLF14 is generally not expressed in blood or its subtypes, it is also possible that it is expressed only transiently, since it is unlikely that transient expression of KLF14 would be detected in a small cross-sectional sample. Future studies in other tissues, such as adipose tissue, will be necessary to assess the possible role of these CpG sites in KLF14 expression, though it is notable that a large genome-wide study of methylation and expression in adipose tissue did not report significant associations (FDR<0.01) involving cg07955995, cg22285878, or KLF14.⁵⁰

Overall, our study is the first to utilize methylome-wide association studies across data sets derived from diverse cell and tissue types to investigate a quadratic relationship between DNAm and age. We found 2 CpG sites that exhibit stable DNAm early in life followed by a rapid increase in DNAm in late life in peripheral whole blood, monocytes, and isolated CD4⁺ T cells. These CpG sites reside in the promoter region of KLF14, which has recently been shown to be involved in CD4⁺ T cell differentiation via the suppression of FOXP3. These findings highlight the importance of DNAm as a means to further our understanding of aging, immunology, and biologic pathways.

Methods

Samples

A total of 11 data sets were analyzed (Table 1). The primary discovery analysis was performed in data from the Grady Trauma project, and the remaining data sets were analyzed to replicate or follow up findings of the primary analysis.

Grady Trauma Project (GTP)

The primary analysis was performed using a subset of 336 African-American individuals ranging in age from 16 to 78 from data collected as part of the GTP, a study investigating the effects of genetic and environmental factors on individuals' response to stressful life events. Participants were recruited from waiting rooms at Grady Memorial Hospital in Atlanta, GA between 2005 and 2008. The Institutional Review Boards of Emory University School of Medicine and Grady Memorial Hospital approved all procedures of the Grady Trauma Project.^51-53

Take Off Pounds Sensibly family study (TOPS)

The TOPS Family Study of Epigenetics included methylation data collected from peripheral blood of 192 individuals. Individual ages ranged from 6 to 85 and each individual belonged to 1 of 7 extended families of Northern European descent. To be included in the study, each nuclear family was required to have 2 obese siblings and at least one parent or sibling who was never obese. The NCBI Gene Expression Omnibus accession number corresponding to the TOPS Family Study is GSE60132.⁵⁴

Multi-Ethnic Study of Atherosclerosis (MESA)

The Multi-Ethnic Study of Atherosclerosis (MESA) collected methylation data for samples of CD4⁺ T cells (MESA-CD4; 214 subjects) and monocytes (MESA-M; 1,202 subjects) isolated from peripheral blood. The age range was 45–79 for MESA-CD4 subjects and 44–83 for MESA-M subjects. The MESA study was conducted to collect population-based information on the prevalence and progression of subclinical cardiovascular disease. Subjects were recruited from 6 sites: Baltimore City and Baltimore County, Maryland; Chicago, Illinois; Forsyth County, North Carolina; Los Angeles County, California; New York, New York; and St. Paul, Minnesota. The NCBI Gene Expression Omnibus accession number corresponding to MESA-CD4 and MESA-M are GSE56581 and GSE56046, respectively.^55,56

Estonian genome center investigation of age-related epigenetics and immune system function in PBL, CD4⁺, and CD8⁺ T cells (EGC-PBL, EGC-CD4, EGC-CD8)

Peripheral blood leukocytes (EGC-PBL) were collected, in addition to CD4⁺ T cells (EGC-CD4) and CD8⁺ T cells (EGC-CD8), which were isolated from peripheral blood from healthy donors of the Estonian Genome Center of the University of Tartu. A total of 101 subjects were in the study, and were divided into a young age group and an old age group whose ranges were 22–34 and 73–84, respectively. After quality control, a total of 296 samples were collected from the 101 subjects: 99 CD4⁺ T cell samples, 100 CD8⁺ T cell samples, and 97 PBL samples. The NCBI Gene Expression Omnibus accession number corresponding to the Estonian Genome Center Investigation of Age-related epigenetics and immune system function in PBL, CD4⁺, and CD8⁺ T cells is GSE59065.⁴⁰

Schizophrenia-related DNAm and gene expression (CTX)

Samples of dlPFC tissue were surgically removed from post-mortem brain samples donated to the NIMH Brain Tissue Collection at the National Institute of Health in Bethesda, Maryland. Our investigation focused on 346 samples from the control group with an age range of 2–85. The NCBI Gene Expression Omnibus accession number corresponding to the Schizophrenia-related DNAm and gene expression data are GSE74193.⁵⁷

Genetics Of Lipid Lowering Drugs and Diet Network (GOLDN)

The GOLDN study recruited families that had at least 2 siblings from the National Heart, Lung, and Blood Institute Family Heart Study from Minneapolis and Salt Lake City. The study collected methylation data from CD4⁺ T cells isolated from peripheral blood samples. After quality control, 991 individuals remained in the analysis. All individuals self-identified as European American and provided written informed consent, and the study protocol was approved by Institutional Review Boards at the University of Minnesota, University of Utah, Tufts University/New England Medical Center, and University of Alabama at Birmingham.

National Institute of Child Health and Human Development (NICHD) brain bank of developmental disorders – cell heterogeneity study (NICHD-G, NICHD-N)

The NICHD study included post-mortem samples from subjects diagnosed with major depressive disorder (MDD) and matched controls. The study collected methylation data from glial cells (n = 50) isolated from prefrontal cortex tissue, in addition to methylation data collected from neurons (n = 50) isolated from prefrontal cortex tissue. The NCBI Gene Expression Omnibus accession number corresponding to NICHD is GSE41826⁵⁸.

Analysis

For each of the studies described above, DNA was collected from tissue samples and was bisulfite-treated for cytosine to thymine conversion and hybridized to the Illumina Infinium HumanMethylation450 (450K) BeadChip and processed according to the instructions of the manufacturer. For GTP and GOLDN, the methylated signal (M) and unmethylated signal (U) as well as detection P-values were obtained from GenomeStudio after processing. For all other studies, data were downloaded from the NCBI Gene Expression Omnibus (GEO).

For all data sets except GOLDN, the cpg.qc function within the R package CpGassoc⁵⁹ was used to perform quality control. Quality control included the removal of probes with missing data for >5% of samples, removal of samples with >5% missing probe sites, and the removal of samples with a detection P-value > 0.001. Methylated and unmethylated signals were then quantile normalized for each data set separately. A β-value was then computed from the methylated signal (M) and unmethylated signal (U) as follows: $\beta =\frac{M}{U+M}$. In GOLDN, quality control was performed in R to exclude samples with >1.5% missing data, β-values with a detection P-value >0.01, CpG sites where probe sequences mapped multiple loci or to a location that did not match the annotation file, and CpG sites where >10% of samples had inadequate intensity. In data sets derived from whole-blood samples (GTP and TOPS), cell type proportions were estimated using the Houseman method as implemented in the R package minfi,^60,61 which utilizes Illumina 450K methylation data from 6 isolated cell types as reference data.¹⁶

For our initial discovery analysis in GTP, the R function cpg.assoc within the CpGassoc package⁵⁹ was used to perform a separate linear regression for each CpG site where β was regressed on age, age², sex, and estimated cell type proportions, and an indicator for row on chip as fixed effects. Chip ID was included as a random effect to absorb potential batch effects. The quadratic term for age was included to allow identification of CpG sites demonstrating an increasing or decreasing rate of change with age. The model used for each of the 11 data sets appears below:(Equation 1) ${\beta }_i={\mathit{\alpha }}_0+{\mathit{\gamma }}_{\mathit{a}\mathit{g}\mathit{e}}\left(AG{E}_i\right)+{\mathbf{\gamma }}_{\mathit{a}\mathit{g}{\mathit{e}}^2}\left(AG{E}_i^2\right)+{\displaystyle \mathbf{\sum }}_{j = 1}^k{\mathit{\gamma }}_{\mathit{j}}co{v}_{ij}+{\mu }_i+{ϵ}_i$(Equation 1) ${\beta }_i$: proportion of cells methylated at a CpG site of the i^th individuals${\mathit{\alpha }}_0$: intercept${\mathit{\gamma }}_{\mathit{j}}$: slope coefficient corresponding to the j^th independent variable$k$: number of covariates in the model. Covariates differ between data sets (Table S4).${\mathit{\mu }}_i$: between-chip error of the i^th individual${ϵ}_i$: within-chip error of the i^th individual

CpG sites were considered to have a significantly increasing rate of methylation change if ${\widehat{\gamma }}_{ag{e}^2}$ had the same sign as ${\widehat{\gamma }}_{age}$ and the P-value corresponding to the t-statistic on ${\widehat{\gamma }}_{ag{e}^2}$ was smaller than the Bonferroni-adjusted level of significance (α = 1.03E-07).

In addition to the quadrilinear model described in the preceding paragraph, we also performed a spline analysis. We chose 19 knots corresponding to 19 quantiles across the age-distribution and fit a model using the same fixed effects and random effects in the above paragraph, with the exclusion of the age-quadratic term and the inclusion of an age-knot interaction term. A separate model was fit for each of the 19 knots across each of the 25,723 CpG sites found to be linearly associated with age in the primary analysis. To compare with our quadrilinear results, we plotted the t-statistics corresponding to the knot that maximized the absolute value of the t-statistic on the age-knot interaction term among the 19 knots across each of 25,723 CpG sites against the t-statistics on the age-quadratic term from the quadrilinear model (Figure S7). The minimum P-value was selected among the 19 P-values corresponding to the age-knot interaction term for each of the 25,723 CpG sites in the spline model, and tested against a Bonferroni threshold of $\alpha =\frac{0.05}{19*25,723}=$ 9.5E-8 to determine significance while adjusting for both the 19 knots considered and the 25,723 CpG sites tested.

For CpG sites demonstrating an increasing rate of change with age in GTP, follow-up analyses were performed on the remaining data sets. Only CpG sites with a significant quadratic term were included in the subsequent analyses, which were performed using the R function lme if random effects were included and lm if no random effects were included (with the exception of GOLDN, where the model was fit using SOLAR⁶² to account for familial relationships between subjects). The response variable β was regressed on age and age² using the model shown in Equation 1(Equation 1) ${\beta }_i={\mathit{\alpha }}_0+{\mathit{\gamma }}_{\mathit{a}\mathit{g}\mathit{e}}\left(AG{E}_i\right)+{\mathbf{\gamma }}_{\mathit{a}\mathit{g}{\mathit{e}}^2}\left(AG{E}_i^2\right)+{\displaystyle \mathbf{\sum }}_{j = 1}^k{\mathit{\gamma }}_{\mathit{j}}co{v}_{ij}+{\mu }_i+{ϵ}_i$(Equation 1) . Chip ID was included as a random effect for GTP, EGC-PBL, MESA-M, EGC-CD8, MESA-T, EGC-CD4, NICHD-G, NICHD-N. Details on fixed and random effect covariates for all data sets are shown in Table S4. The MESA-T and MESA-M models included estimated measures of contamination by other cell types, namely, proportions of B-cells, monocytes, natural killer cells, and neutrophils, which were provided with GEO data sets GSE56581 and GSE56046⁵⁶. Likewise, the CTX model included estimated cell type proportions (embryonic stem cells, neural progenitor cells, mature adult neurons, non-neuronal cells, and ES-derived dopamine neurons), which were provided with the GEO data set (GSE74193).⁵⁷

Based on Equation 1(Equation 1) ${\beta }_i={\mathit{\alpha }}_0+{\mathit{\gamma }}_{\mathit{a}\mathit{g}\mathit{e}}\left(AG{E}_i\right)+{\mathbf{\gamma }}_{\mathit{a}\mathit{g}{\mathit{e}}^2}\left(AG{E}_i^2\right)+{\displaystyle \mathbf{\sum }}_{j = 1}^k{\mathit{\gamma }}_{\mathit{j}}co{v}_{ij}+{\mu }_i+{ϵ}_i$(Equation 1) , we derived the instantaneous slope between DNAm and age at ages 50, 60, and 70 y using the following computations (Tables 2 and 3):${\gamma }_{age = 50}:{\gamma }_{age}+2{\gamma }_{ag{e}^2}\times 50$${\gamma }_{age = 60}:{\gamma }_{age}+2{\gamma }_{ag{e}^2}\times 60$${\gamma }_{age = 70}:{\gamma }_{age}+2{\gamma }_{ag{e}^2}\times 70$

For comparisons of variance in DNAm between older and younger subjects, GTP, TOPS, CTX, EGC-PBL, EGC-CD4, and EGC-CD8 were each partitioned into an older (≥ 73 years) and a younger group (≤ 34 years). Because the age ranges of MESA-M and MESA-T were 45–79 and 44–83 years, respectively, the younger group in these data sets comprised individuals ≤ the median age (58 for MESA-M and 60 for MESA-T) and the older age group comprised individuals > the median age. We used the median (23 for both NICHD-G and NICHD-N) to separate NICHD-N and NICHD-G into young and old age groups in a similar fashion as MESA because NICHD-G and NICHD-N both only had 2 individuals aged ≥ 73 years. The variances of DNAm at cg07955995 and cg22285878 were computed for the young age group and the old age group and an F statistic was computed as the ratio of these variances (Table 4).

As technical validation to address the possibility that artifacts were driving our key results, we performed targeted analysis of samples from the 2 oldest and 2 of the youngest individuals in the GTP cohort. We designed an assay for the promoter region of KLF14 that included 20 additional CpG sites flanking cg07955995 and cg22285878. We used commercially available standards to validate the assay across the range of possible methylation values (0–100%) and generated standard curves of observed vs. expected DNAm at key CpG sites (Figure S8). We then used Sequenom EpiTyper technology to measure DNAm at 22 CpG sites for the 2 older vs. the 2 younger subjects.

Disclosure of potential conflicts of interest

The authors declare no competing financial interests.

Author contributions

N.J. performed all data analyses except those specific to GOLDN, wrote the main manuscript text, and prepared all Figures and Tables. H.W. performed all data analyses in the GOLDN study, and contributed figure panels from these analyses. A.S. contributed to study design and data collection in the GTP and provided a critical read of the manuscript. S.N. and A.S. performed the Sequenom EpiTyper assay for the technical validation of results. D.M.A., D.K.A., and S.A. contributed to study design and data collection in the GOLDN study. S.A. supervised data analyses in the GOLDN study and provided a critical read of the manuscript. K.C. conceived of the study, supervised all data analyses, contributed to study design and data collection in the GTP, and made critical revisions to the manuscript. All authors reviewed and approved the manuscript.

Acknowledgments

We gratefully acknowledge Drs. Kerry Ressler and Elisabeth Binder for use of data from the Grady Trauma Project. The Grady Trauma Project was supported by the National Institutes of Health Grants MH071537 and MH096764 (Ressler) and MH085806 (A.K.S.). This work was also supported by a NARSAD YI award (A.K.S.) and the Howard Hughes Medical Institute (Ressler). GOLDN was supported by National Institutes of Health grant HL104135 (D.K.A.). We are also very appreciative of Zachary V. Johnson for his help formatting Fig. 1.