Biological stability of DNA methylation measurements over varying intervals of time and in the presence of acute stress

Figures & data

Figure 1. Description of test-retest scenarios.

Note: Blue lines and dots indicate the NoStress session and red lines and dots indicate the Stress session. CrossSessionT1 was the only scenario in which probes across sessions were used, which is represented by a purple open dot. Each open dot represents a time point that was used to calculate the probe stabilities under the corresponding scenarios. All scenarios included two time points to calculate probe stabilities, apart from NoStressT1-2-3-4, in which four time points were used.

Figure 2. Associations of DNAm probe stability with sample characteristics.

Note: Heat-scatter plots are set on a three-colour gradient, with red being the densest, green being moderately dense, and blue being the least dense. White indicates no points in the region. A) Differences in probe stability associated with the method of calculating probe ICC values. B) Differences in probe stability associated with either adjusting or not adjusting the ICC values based on the proportions of immune cell subtypes. C) Four panels showing differences in probe stability associated with varying sample sizes (x-axis ICC with N = 14 samples; y-axis: ICC with N = 31 samples). D) Four panels show differences in probe stability associated with varying numbers of repeated measurements.

Figure 3. Stability of DNAm probes in test-retest scenarios.

Note: All probes with ICC values less than 0.01 were excluded from analyses in this figure. a) Distribution of probe stability for each test-retest scenario. The vertical dashed line indicates the median probe stability (excluding probes with ICC < 0.01). b) Cluster dendrogram showing the similarity of test-retest scenarios (based on probe ICC values).

Table 1. Test-retest scenario descriptive statistics.

Scenario	Mean	SD	Median	Skew	Kurtosis
NoStressT1–2	0.47	0.22	0.47	0.08	−0.70
NoStressT1–3	0.46	0.22	0.45	0.12	−0.66
NoStressT1–4	0.35	0.22	0.31	0.70	−0.09
NoStressT3–4	0.49	0.25	0.50	−0.01	−0.92
NoStressT1-2-3-4	0.44	0.21	0.44	0.18	−0.42
StressT1–2	0.60	0.25	0.64	−0.47	−0.74
StressT1–3	0.60	0.24	0.64	−0.50	−0.64
StressT1–4	0.61	0.24	0.65	−0.55	−0.56
StressT3–4	0.59	0.25	0.62	−0.40	−0.76
CrossSessionT1	0.33	0.22	0.28	0.85	0.13

All descriptive statistics were calculated excluding probes with stabilities less than 0.01.

Table 2. Test-retest scenario probe stability distributions.

Scenario	ICC<0.01	0.01<ICC<0.50	0.50<ICC<0.75	0.75<ICC<0.90	ICC>0.90
NoStressT1–2	7.09%	51.50%	30.00%	9.00%	2.41%
NoStressT1–3	6.27%	54.44%	29.84%	7.61%	1.84%
NoStressT1–4	13.31%	66.79%	15.24%	3.34%	1.32%
NoStressT3–4	10.25%	45.41%	28.96%	11.25%	4.13%
NoStressT1-2-3-4	3.05%	59.23%	30.64%	5.18%	1.90%
StressT1–2	9.70%	30.11%	30.16%	21.51%	8.51%
StressT1–3	9.93%	28.66%	31.17%	21.51%	8.72%
StressT1–4	11.28%	27.20%	31.66%	21.82%	8.03%
StressT3–4	10.11%	31.43%	30.64%	19.93%	7.89%
CrossSessionT1	13.51%	68.05%	13.78%	3.22%	1.44%

Percentage of probes in each scenario stratified by ICC value. The total number of probes in each scenario was 843,328.

Figure 4. DNAm probe set enrichment analyses.

Note: All probes with ICC values less than 0.01 were excluded from analyses in this figure. The dashed horizontal lines indicate the median probe stability in each scenario (excluding probes with ICC < 0.01) and the y-axes indicate probe ICC values. A) Common DNAm clock probe set and immune cell estimator probe set distributions plotted against the median probe stability of each test-retest scenario. Horvath, Hannum, PhenoAge, PC Clock, and PACE are probes used in epigenetic clocks, while IDOL-6 and IDOL-Extended are probes used in previous and current immune cell estimators. B) Biologically relevant probe set distributions plotted against the median probe stability of each test-retest scenario. The core circadian clock genes included CLOCK, NPAS2, BMAL1, BMAL2, PER1, PER2, PER3, CRY1, and CRY2. C) Comparisons of probe stability and probe correlations with gene expression values for all test-retest scenarios. The x-axis indicates the correlation between the probe and gene expression values, and the y-axis indicates the ICC value of the probe. Blue dots are probes that were not significantly associated with gene expression values and red dots are probes that were significantly associated with gene expression values. All included probes were located in TSSs.

Figure 5. Effects of acute psychosocial stress on highly stable DNAm probes.

Note: a) Manhattan plot showing T1 – T2 stress effect significance (y-axis) and chromosomal location of the probe (x-axis). The dotted line indicates a significance threshold of 0.05. b) Manhattan plot showing T1 – T3 stress effect significance (y-axis) and chromosomal location of the probe (x-axis). The dotted line indicates a significance threshold of 0.05 and the solid line indicates a significance threshold of 3.87E-05. C) Manhattan plot showing T1 – T4 stress effect significance (y-axis) and chromosomal location of the probe (x-axis). The dotted line indicates a significance threshold of 0.05 and the solid line indicates a significance threshold of 3.87E-05. Labelled GSR probes cg19913448 and cg16183701 included.

Table 3. Effects of acute stress on highly stable DNAm probes.

Probe	Gene	β-Intercept	β Time 4 Estimate	M-Intercept	M Time 4 Estimate	P-Value
cg19913448	GSR	0.026	−3.21E–03	−5.276	−0.191	1.30E–10
cg16183701	GSR	0.017	−1.66E–03	−5.908	−0.156	1.06E–07
cg13722517	YARS2	0.030	−3.58E–03	−5.093	−0.211	1.63E–07
cg26031613	KLC1	0.045	−3.77E–03	−4.972	−0.210	1.24E–06
cg25114630	CHSY1	0.044	−3.44E–03	−5.236	−0.147	3.09E–06
cg00635556	EI24	0.021	−2.08E–03	−5.611	−0.155	3.13E–06
cg13353337	LZTR1	0.126	−5.34E–03	−3.903	−0.190	8.14E–06
cg25613667	WDR55	0.164	−7.77E–03	−3.376	−0.157	1.13E–05
cg00287122	RAB19	0.046	−5.96E–03	−4.444	−0.224	1.43E–05
cg24400183	ZNF563	0.014	−2.61E–03	−6.128	−0.254	1.73E–05
cg03007623	CDKN2D	0.058	−2.29E–03	−5.206	−0.185	1.92E–05

All Bonferroni significant TSS probes and their associated genes obtained when comparing T4 and T1 of the stress-session. Intercepts and T4 estimates are given as both β-values and M-values. P-values are given with respect to M-values. Significance threshold: P < 3.87E–05.

Figure 6. ELA vs. Non-ELA stabilities of DNAm probes in the presence of acute psychosocial stress.

Note: Heat-scatter plots are set on a three-colour gradient, with red being the densest, green being moderately dense, and blue being the least dense. White indicates no points in the region. Probe stability in the presence of acute psychosocial stress for non-ELA individuals (x-axis) and ELA individuals (y-axis) in four different scenarios.

Data availability statement

The code used to perform all the analyses is available at: https://github.com/abnerapsley1/DNAmStabilityMeasurements. Data used in this publication have been deposited in NCBI’s Gene Expression Omnibus (GEO) and are accessible through the GEO Series accession number GSE227815 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE227815).