Epigenome-wide meta-analysis of PTSD symptom severity in three military cohorts implicates DNA methylation changes in genes involved in immune system and oxidative stress.

Epigenetic factors modify the effects of environmental factors on biological outcomes. Identification of epigenetic changes that associate with PTSD is therefore a crucial step in deciphering mechanisms of risk and resilience. In this study, our goal is to identify epigenetic signatures associated with PTSD symptom severity (PTSS) and changes in PTSS over time, using whole blood DNA methylation (DNAm) data (MethylationEPIC BeadChip) of military personnel prior to and following combat deployment. A total of 429 subjects (858 samples across 2 time points) from three male military cohorts were included in the analyses. We conducted two different meta-analyses to answer two different scientific questions: one to identify a DNAm profile of PTSS using a random effects model including both time points for each subject, and the other to identify a DNAm profile of change in PTSS conditioned on pre-deployment DNAm. Four CpGs near four genes (F2R, CNPY2, BAIAP2L1, and TBXAS1) and 88 differentially methylated regions (DMRs) were associated with PTSS. Change in PTSS after deployment was associated with 15 DMRs, of those 2 DMRs near OTUD5 and ELF4 were also associated with PTSS. Notably, three PTSS-associated CpGs near F2R, BAIAP2L1 and TBXAS1 also showed nominal evidence of association with change in PTSS. This study, which identifies PTSD-associated changes in genes involved in oxidative stress and immune system, provides novel evidence that epigenetic differences are associated with PTSS.

Introduction

Posttraumatic stress disorder (PTSD) can develop in some people following trauma and results in severe symptoms including intrusive thoughts, avoidance of trauma-related stimuli, negative cognitive and mood changes, and hyperarousal that disturb mental and physical wellbeing[1]. Although a vast majority of the population experiences trauma to at some point in their life[2], PTSD prevalence is only 6.8% among US population[3]. Because only a fraction of people who experience trauma go on to develop PTSD, it is important to understand the factors that increase risk for the disorder or contribute to its symptom severity. DNA methylation (DNAm), an epigenetic modification, is one such factor involved in adaptation to traumatic stress[4-6].

Epigenome-wide association studies (EWASs) of PTSD have discovered differentially methylated CpGs in genes related to neuronal and immune pathways[7-14]. The majority of these studies have a cross-sectional design; DNAm being examined at a single time-point. In addition to previously identified, cross-sectional associations, understanding whether DNAm changes as individuals develop PTSD or experience changes in PTSD symptom severity (PTSS) is crucial. Recently, two longitudinal studies reported DNAm changes associated with PTSD development in individuals exposed to combat trauma[11, 12]. Both studies used modest sample sizes of 93 and 266 subjects, respectively, and the HumanMethylation450 BeadChip to identify CpGs and differentially methylated regions (DMRs) associated with PTSD development. Rutten et al. observed lower DNAm levels in PTSS at genomic regions in ZFP57, RNF39 and HIST1H2APS2[11]. Snijder et al. reported contributions from the immune system through the HLA locus, HEXDC and MAD1L1 in development of PTSD, using different subjects of the three military cohorts that participated in this study[12].

Building on the prior work of Snijder et al.[12], this study features a larger sample size, a denser and more comprehensive array, and additional statistical models to gain more insight into the epigenetics of PTSD. We first performed a meta-analysis in 858 samples (429 subjects with pre- and post- deployment samples) to identify CpGs and DMRs that associate with PTSS. Then, we conducted a second meta-analysis in 429 subjects to identify associations between DNAm and change in PTSS pre- to post-deployment. Finally, we evaluated CpGs identified in previously published Psychiatric Genomics Consortium PTSD Workgroup (PGC-PTSD) EWAS[10-14] in targeted longitudinal analyses. We focus on PTSS in order to overcome case-control selection bias, as some participants had elevated PTSD symptom scores before deployment, and to gain statistical power through the use of continuous variables[15].

Methods

Cohorts

This study includes 429 subjects from three military cohorts that are presented in Table 1: Marine Resilience Study (MRS), Army Study to Assess Risk and Resilience in Servicemembers (Army STARRS), and Prospective Research in Stress-related Military Operations (PRISMO). Details of each cohort are in the Supplement. PTSS was measured by each individual study pre- and post- deployment. All participants gave informed consent, and all studies were approved by respective institutional review boards.

Table 1:

Demographic and clinical characteristics of MRS, Army STARRS, and PRISMO

	Cases	Controls	p value	Total
Number
MRS	64	63	-	127
Army STARRS	92	92	-	184
PRISMO	43	75	-	118
Age, mean (SD)
MRS	22.16 (2.35)	21.97 (2.12)	0.64	22.06 (2.24)
Army STARRS	24.41 (4.85)	24.52 (4.86)	0.88	24.47 (4.84)
PRISMO	27.84 (9.69)	27.07 (8.74)	0.67	27.35 (9.06)
PTSD symptoms pre-deployment, mean (SD)
MRS, PCL-17	24.88 (8.69)	20.00 (4.54)	0.0001	22.46 (7.34)
Army STARRS, PCL-6	6.92 (1.34)	6.48 (0.95)	0.01	6.7 (1.18)
PRISMO, SRIP	29.00 (4.14)	26.72 (4.09)	0.005	27.55 (4.24)
PTSD symptoms post-deployment, mean (SD)
MRS, PCL-17	49.23 (11.17)	22.03 (6.06)	< 2.2e-16	35.74 (16.34)
Army STARRS, PCL-17	43.83 (16.04)	20.74 (3.77)	< 2.2e-16	31.18 (16.01)
PRISMO, SRIP	42.14 (4.67)	27.11 (4.76)	< 2.2e-16	32.56 (8.66)
Self-reported Race/Ethnicity, N (%)
MRS			0.98
European	43 (67.2)	45 (71.4)	-	88 (69.3)
African American	3 (4.8)	2 (3.2)	-	5 (3.9)
Latino	9 (14.0)	7 (11.1)	-	16 (12.6)
East Asian	1 (1.5)	1 (1.6)	-	2 (1.6)
Other	8 (12.5)	8 (12.7)	-	16 (12.6)
Army STARRS			0.90
European	61 (66.3)	63 (68.5)	-	124 (67.4)
African American	10 (10.9)	11 (12.0)	-	21 (11.4)
Other	21 (22.8)	18 (19.5)	-	39 (21.2)
PRISMO			1.00
European	43 (100)	75 (100)	-	118 (100)
Smoking Score pre-deployment, mean (SD)
MRS	−7.69 (14.14)	−8.20 (13.89)	0.84	−7.94 (13.96)
Army STARRS	−6.30 (18.98)	−10.26 (16.43)	0.13	−8.28 (17.81)
PRISMO	2.42 (21.08)	2.96 (24.42)	0.90	2.77 (23.16)
Smoking Score post-deployment, mean (SD)
MRS	−7.21 (16.40)	−9.33 (14.43)	0.44	−8.26 (15.43)
Army STARRS	−4.78 (18.33)	−9.60 (16.64)	0.063	−7.19 (17.62)
PRISMO	4.17 (21.68)	3.53 (25.35)	0.88	3.76 (23.98)
Early life trauma, mean (SD)
MRS, CTQ[a]	41.55 (13.08)	36.85 (11.27)	0.047	39.22 (12.39)
Army Starrs, NSS	6.76 (2.80)	6.21 (2.05)	0.13	6.48 (2.46)
PRISMO, ETI[b]	5.09 (3.29)	3.37 (3.23)	0.007	4.01 (3.35)

Each study used different scales for PTSD, the corresponding scales are included in the row names: CTQ, Childhood Trauma Questionnaire; ETI, Early Trauma Inventory; NSS, The Army STARRS New Soldier Survey; PCL-17, 17-item PTSD Checklist; PCL-6, 6-item PTSD Checklist; SRIP, Self-Report Inventory for PTSD. SD: standard deviation. Smoking Score was estimated using DNAm data. The p-values were computed using t-test (for continuous variables) and Fisher’s exact test (for categorical variables) for comparison of PTSD case and control groups.

Missing data for 18 subjects.

Missing data for 2 subjects.

Quality Control (QC) Procedures

Whole blood DNAm was measured using the Illumina MethylationEPIC BeadChip. The same QC pipeline was applied separately to each of the cohorts. We used the R package CpGassoc to filter out samples with probe detection call rates <90% and an average intensity value of either <50% of the experiment-wide sample mean or <2000 arbitrary units (AU)[16]. We set low quality probes (detection p-values >0.01) as missing. We filtered out probes that were missing for >10% of samples within studies. We removed cross hybridizing probes[17]. A total of 820 498 probes passed QC in all cohorts and were included in our analyses. We performed single-sample Noob (ssNoob) normalization using R package minfi[18]. To remove chip and positional batch effects, we applied ComBat, protecting age and PTSD status[19]. We used logit transformed beta values (M-values) in our analyses[20].

For each sample, cellular heterogeneity (i.e. the proportion of CD8+T, CD4+T, natural killer (NK), B cells, monocytes and neutrophils) was predicted using the Robust Partial Correlation (RPC) method implemented in Epidish[21] using the reference data reported by Salas and colleagues[22]. Ancestry principal components (PCs) were generated from DNAm, following the method described by Barfield et al.[23], as previously implemented[24]. The components that correlate most with self-reported race/ethnicity (PCs 2–3) were used to adjust for ancestry (Supplementary Figure 1)[23, 24]. DNAm data was used to estimate smoking information as previously described[25]. Computation of ancestry PCs and smoking scores are described in the Supplement.

Statistical analysis

Since different measures of PTSS were used across studies and timepoints, heterogeneity was minimized by rescaling PTSS using a min-max normalization method to scale the range in [0, 1]. To identify CpGs associated with PTSS (Meta-Analysis 1), we used a linear mixed model with DNAm values at both time points as the dependent variable, PTSS at both time points as a main effect, and a random intercept for subject. Age, CD8+T, CD4+T, NK, B cell, and monocyte cell proportions, and ancestry PCs derived from DNAm data were included as covariates. To identify CpGs associated with change in PTSS (Meta-Analysis 2), we conducted a longitudinal analysis using a linear regression model, where post-deployment DNAm was modeled as a function of change in PTSS while adjusting for pre-deployment DNAm, PCs for ancestry, and changes in age (i.e. time passed between pre- and post-deployment data collection), CD8T+, CD4T+, NK, B cell, and monocyte cell proportions.

Meta-analysis of cohorts was performed using weighted sum of z-scores method, as Cochran’s Q test did not show substantial heterogeneity[26-28]. To control for multiple testing, we used the epigenome-wide significance threshold proposed for the MethylationEPIC BeadChip (p < 9.0E-08)[29]. Post-hoc sensitivity analysis explored the possible confounding effects of smoking (Meta-Analysis 1) and changes in smoking status (Meta-Analysis 2) by including DNAm derived smoking scores as a covariate. A post-hoc sensitivity analysis also examined the impact of early life trauma on PTSS-related DNAm changes by including early life trauma burden as a covariate in the models. To assess the variability of the PTSS-associated CpGs over time, a third post-hoc analysis in which PTSS was excluded from the model was performed. Post-hoc power analyses were performed as described in the Supplement.

In addition to our two primary meta-analyses, we also performed DMR analyses to identify i) DMRs associated with PTSS (i.e. using the same framework as Meta-Analysis 1), and ii) DMRs associated with change in PTSS by conditioning post-deployment DNAm on baseline DNAm (i.e. using the same framework as Meta-Analysis 2). We used DMRcate to calculate the significance of regions (at least 2 probes within 1kb of each other) based on EWAS summary statistics[30]. This included at least one strongly associated CpG site (p <0.0001). DMRs with a Stouffer transformed false-discovery rate (FDR) of 5% across the region were considered significant.

Finally, to provide additional insight into earlier findings, we evaluated the PTSD-associated CpGs from previous PGC-PTSD EWAS[10-14] via targeted analyses using the framework for Meta-Analysis 1 described above. Bonferroni correction was used to account for multiple comparisons.

Blood-Brain Correlations

Correlation between blood and brain DNAm of associated CpGs was examined using the IMAGE-CpG database[31]. Specifically, this database maintains Spearman correlation coefficients (rho) and associated p-values for CpGs from 27 individuals with paired blood and live brain samples.

Genetic influence of CpGs associated with PTSS

To evaluate the effect of nearby polymorphisms on DNAm levels of CpGs associated with PTSS, we leveraged cis-methylation quantitative trait locus (cis-meQTL) data from BIOS QTL browser[32]. Cis-meQTL, here was described as the correlation between a CpG and a SNP within 250 kb with a CpG-level FDR threshold of 5% (p ≤ 1.38E-04)[32]. To evaluate whether the meQTLs from BIOS QTL browser had similar effects in our study, we tested the associations between post-deployment methylation levels of CpGs and their respective meQTL SNPs in all three cohorts, using linear regression models that adjust for the cohorts. For CpGs with identified meQTLs, we performed post-hoc sensitivity analyses by adding their respective meQTL SNP as a covariate to Meta-Analysis 2.

Pathway enrichment analysis

We conducted exploratory gene ontology (GO) pathway enirchment analyses using MissMethyl[33] and including variably methylated probes (VMPs) from any of the 3 cohorts that nominally associate with PTSS (Meta-Analysis 1) or change in PTSS (Meta-Analysis 2). FDR of 5% was used to define significant pathways. To identify VMPs, we first calculated longitudinal DNAm differences for all CpG sites (Δβ = βpost − βpre); then we computed the median value of absolute DNAm differences (|median(Δβ)|), defining those with >1% difference from the median (|median(Δβ)|≥0.01) as variable.

Code availability

The scripts generated to perform the Meta-Analysis 1, Meta-Analysis 2, and DMR analysis are available in https://github.com/PGC-PTSD-EWAS/PGC-PTSD-Longitudinal-Analysis.

Results

Demographics of the cohorts

Demographic and clinical information of participants from all studies (total N subjects = 429) are summarized in Table 1. All participants were male and were primarily of European ancestry (N = 330, 77%). Age and smoking did not differ between PTSD cases and trauma-exposed controls. Eventual cases had higher pre-deployment PTSS, compared to controls in all three cohorts, potentially due to higher rates of early life trauma (Table 1).

Meta-Analysis 1: Evaluating CpGs associated with PTSS

We identified four significant CpGs (Table 2, Figure 1A, Supplementary Table 1). These sites were located near the coagulation factor II thrombin receptor (F2R), canopy FGF signaling regulator 2 (CNPY2), Brain-specific angiogenesis inhibitor 1-associated protein 2-like protein 1 (BAIAP2L1), and thromboxane A synthase 1 (TBXAS1) genes. For all sites, lower DNAm levels associated with higher PTSS (Supplementary Figures 2–4). These observed associations were not due to ageing (Supplementary Table 1). All four sites remained significant with the same direction of association in our sensitivity analysis adjusted for smoking score (Supplementary Table 1). Only the CpG in CNPY2 did not exceed the genome-wide significance threshold when we adjusted for early life trauma (Supplementary Table 1), indicating that the PTSS-associated DNAm changes were largely uninfluenced by early life trauma burden.

Table 2:

Genome-wide significant CpG sites (p < 9.0E-08) associated with PTSS

CpG	Location	Gene	Features	Z	p	Z 2	p 2
cg11627632	chr5:76011698	F2R	TSS200	−5.95	2.75E-09	−2.38	0.02*
cg12961546	chr12:56709730	CNPY2	5’UTR; 1stExon	−5.62	1.95E-08	−1.78	0.07
cg00277769	chr7:97922759	BAIAP2L1	3’UTR	−5.41	6.39E-08	−2.55	0.01*
cg03604364	chr7:139705703	TBXAS1	Body	−5.41	6.47E-08	−2.66	0.01*

Genome-wide significant results of the EWAS meta-analysis 1 of three cohorts (N samples = 858). Association analyses of each cohort are based on a random intercept model with a random effect of subject. Z2 and p2 represents the statistics of meta-analysis 2. The sites that were nominally significant in meta-analysis 2 were indicated by an asterisk (*) in column p2.

Figure 1.

PTSS associates with DNAm across the genome.

A) Manhattan plot for Meta-Analysis 1 across 3 cohorts (N samples = 858). Association analyses of each cohort are based on a random intercept model with a random effect of subject. B) Manhattan plot for Meta-Analysis 2 across 3 cohorts (N subjects = 429). Association analyses of each cohort are conducted by conditioning post-deployment methylation on baseline DNAm. The x-axis is the chromosomal location of each site across the genome. The y-axis is the −log10 of the unadjusted p-value for the association with PTSD symptom severity. The red line indicates genome-wide EWAS statistical significance at p < 9.0E-8.

Of the four CpGs, blood DNAm levels of cg00277769 in BAIAP2L1 and cg03604364 in TBXAS1 were correlated with brain DNAm levels (Supplementary Table 2).

In addition, we identified 88 DMRs that were significantly associated with PTSS (Stouffer p <0.05; Supplementary Table 3). All DMRs except two were still significant with the same direction of association in the sensitivity analysis adjusted for smoking score (Supplementary Table 3).

Meta-Analysis 2: Evaluating CpGs associated with change in PTSS

We identified 47 660 CpGs where post-deployment methylation was nominally (p <0.05) associated with change in PTSS; however, none exceeded the genome-wide significance threshold (Figure 1B). Of the four significant CpGs from Meta-Analysis 1, three associated with change in PTSS after deployment at the p <0.05 level, all with the same direction of association (cg11627632 in F2R, cg00277769 in BAIAP2L1, cg03604364 in TBXAS1; Table 2).

In addition, we identified 15 DMRs whose post-deployment DNAm significantly differed from pre-deployment DNAm based on their change in PTSS (Stouffer p <0.05; Supplementary Table 4). All DMRs were still significant with the same direction of effect in the sensitivity analysis that accounted for smoking (Supplementary Table 4). Of those 15 DMRs that were associated with change in PTSS, two DMRs located on genes OTUD5 and ELF4 were also associated with PTSS, which was not more than would be expected by chance (Fisher’s exact test p = 0.10; Table 3, Supplementary Figure 5).

Table 3:

Differentially methylated regions (DMRs) associated with both PTSS and change in PTSS

	DMR analysis with Meta-Analysis 1 Framework				DMR analysis with Meta-Analysis 2 Framework
Overlapping Genes	Position	N CpGs	p-value	Direction[a]	Position	N CpGs	p-value	Direction[b]
OTUD5	chrX:48814580–48814955	8	2.46E-05	+	chrX:48814205–48815125	13	9.53E-09	+
ELF4	chrX:129244725–129244742	3	1.63E-02	+	chrX:129244725–129245245	9	8.50E-07	+

Chr: chromosome, Start/End: start and end position of DMR (hg19), N CpG: the number of CpGs measured within the DMR, Overlapping Genes: genes that span the DMR. p-value: Stouffer p value of the DMR. Direction: direction of the relation between change in PTSD symptom severity and DNAm levels.

increased methylation with increase in PTSS.

increased methylation with increased change in PTSS.

Genetic effects of CpGs associated with PTSS

Out of four CpGs associated with PTSD symptom severity, three (cg11627632 in F2R, cg00277769 in BAIAP2L1, cg03604364 in TBXAS1) correlated with at least one nearby SNP within 250 kb of the CpG, according to BIOS QTL browser (Supplementary Table 5). However, in our study, only cg03604364 in TBXAS1 was associated with its meQTL SNP (rs3779130) from BIOS QTL browser (Supplementary Table 5). Controlling for genotypes in the Meta-Analysis 1 to evaluate the effect of SNPs from BIOS QTL browser did not substantially affect the observed results, and all three CpGs maintained genome-wide significance (Note that BAIAP2L1 does not reach genome-wide significance for one of the two SNPs; Supplementary Table 6).

Evaluation of CpGs from previously published PGC-PTSD EWAS

We compared our results to previous PGC-PTSD EWAS results[10-14] to evaluate published genome-wide significant CpGs (Table 4). Out of 31 CpGs from five studies, we observed nominal evidence of association for three CpGs: cg05575921 in AHRR (Meta-Analysis 1), cg26703534 in AHRR (Meta-Analysis 1), and cg19534438 on G0S2 (in both Meta-Analysis 1 and 2). The direction of association was the same as that reported in the original studies for cg05575921 and cg26703534, but opposite for cg19534438 (Table 4).

Table 4:

Evaluation of CpGs identified in previously published PGC-PTSD EWAS

CPG	Reference	Tissue	Gene	Meta-Analysis 1		Meta-Analysis 2
				z	p	z	p
cg19534438	Logue et al., 2020	Blood	G0S2	−3.10[a]	0.002	−2.37[a]	0.018
cg04130728	Logue et al., 2020	Brain	CHST11	−1.63	0.103	−1.95	0.051
cg14911689	Rutten et al., 2018	Blood	NINJ2	−0.98	0.33	0.78	0.44
cg24406898	Rutten et al., 2018	Blood	COL1A2	−0.75	0.45	−0.76	0.45
cg01516881	Rutten et al., 2018	Blood	DUSP22	−0.66	0.51	0.12	0.90
cg11763394	Rutten et al., 2018	Blood	PAX8	−0.60	0.55	0.61	0.54
cg11235426	Rutten et al., 2018	Blood	DUSP22	−0.55	0.58	0.27	0.79
cg06417478	Rutten et al., 2018	Blood	HOOK2	−0.49	0.62	−1.07	0.29
cg04657146	Rutten et al., 2018	Blood	HOOK2	−0.33	0.74	1.37	0.17
cg26654770	Rutten et al., 2018	Blood	NINJ2	−0.31	0.75	−0.36	0.72
cg18110333	Rutten et al., 2018	Blood	DUSP22	−0.30	0.77	0.12	0.90
cg21548813	Rutten et al., 2018	Blood	DUSP22	−0.21	0.83	0.81	0.42
cg05785424	Rutten et al., 2018	Blood	intergenic	−0.13	0.89	−0.19	0.85
cg07249765	Rutten et al., 2018	Blood	SDK1	−0.10	0.92	−0.55	0.58
cg10075506	Rutten et al., 2018	Blood	MYT1L	−0.09	0.92	0.07	0.94
cg03517284	Rutten et al., 2018	Blood	HIST1H2APS2	0.06	0.95	0.47	0.64
cg11738485	Rutten et al., 2018	Blood	HOOK2	−0.04	0.97	0.35	0.73
cg03395511	Rutten et al., 2018	Blood	DUSP22	−0.03	0.98	1.12	0.26
cg04657146	Rutten et al., 2018	Blood	intergenic	0.00	1.00	0.25	0.80
cg05575921	Smith et al., 2020	Blood	AHRR	−3.36	0.001 [b,c]	−1.30	0.194
cg26703534	Smith et al., 2020	Blood	AHRR	−2.51	0.012 c	−0.50	0.614
cg25648203	Smith et al., 2020	Blood	AHRR	1.56	0.120	0.55	0.579
cg21161138	Smith et al., 2020	Blood	AHRR	−1.04	0.298	−1.20	0.231
cg05901543	Snijders et al., 2020	Blood	CDH15	−1.39	0.165	−0.84	0.399
cg18917957	Snijders et al., 2020	Blood	CTRC	1.05	0.292	0.56	0.574
cg12169700	Snijders et al., 2020	Blood	MAD1L1	0.70	0.484	0.19	0.852
cg05656210	Snijders et al., 2020	Blood	intergenic	0.49	0.627	0.11	0.915
cg16956686	Snijders et al., 2020	Blood	SDK1	0.38	0.706	0.06	0.953
cg20756026	Snijders et al., 2020	Blood	HEXDC	0.23	0.816	1.21	0.225
cg19577098	Uddin et al., 2018	Blood	HGS	1.62	0.106	1.10	0.272
cg23637605	Uddin et al., 2018	Blood	NRG1	0.04	0.969	−0.60	0.550

CpGs with nominal significance of association in our study are shown in bold.

Opposite direction of effect from the original study.

Significant after Bonferonni correction for 31 CpGs (0.05/31 = 0.0016).

Significant when 82 duplicate samples (41 subjects) from PRISMO cohort were removed.

It is important to note that 41 subjects of the PRISMO cohort who participated this study were also included in some of the previous studies[11, 12, 14]. To perform an independent analysis for the AHRR CpGs, cg05575921 and cg26703534, we repeated our targeted meta-analysis by removing these 41 subjects. The two CpGs were still significant and showed a same direction of effect (cg05575921, p = 0.002, z = −3.15; cg26703534, p = 0.004, z = −2.88). Of note, cg05575921 remained significant after multiple test correction for 31 genome-wide significant CpGs identified in previous PGC-PTSD EWASs.

We identified 157 809 VMPs, of which 16 974 associated with PTSS and 8569 with change in PTSS. GO enrichment analysis revealed 173 biological processes enriched in CpGs that associate with PTSS and 9 with change in PTSS (FDR < 0.05). Many processes relate to immune function (Supplementary Table 7), and 7 processes, including leukocyte migration, are enriched in both PTSS and change in PTSS analyses.

Discussion

It is unclear why some develop PTSD after trauma while others do not[34-36]. A likely underlying mechanism is epigenetic alteration, which links environmental circumstances and experiences to biological response. Here, we employed two study designs: i) to identify CpGs that are associated with PTSS measured at pre- and post-deployment (Meta-Analysis 1), and ii) to investigate associations of DNAm with change in PTSS pre-to-post deployment (Meta-Analysis 2).

The first meta-analysis showed that increased PTSS is associated with lower methylation levels at four CpGs located in F2R, CNPY2, BAIAP2L1 and TBXAS1. The CpGs in F2R, BAIAP2L1 and TBXAS1 were also nominally associated with change in PTSS, and therefore show only a small change in pre- to post-deployment DNAm. F2R participates thrombotic response regulation, and is involve in mediating the cross-talk between coagulation and inflammation[37]. A study that investigated gene-expression levels in peripheral blood samples reported lower F2R expression in PTSD cases[38]. Together, this information suggests that transcriptional regulation of F2R may contribute to PTSD, conceivably by modulating immune response. CNPY2 functions in the endoplasmic reticulum (ER) and plays a key role in the transitioning from the non‐stressed to the stressed state[39]. In addition, CNPY2 contributes to central nervous system development by stimulating neurite outgrowth[40]. This evidence suggests the possible role of CNPY2 in oxidative stress and central nervous system development[40], which are critical pathways in PTSD[41]. BAIAP2L1 promotes cell proliferation by stimulating the EGFR-ERK pathway[42] and regulating short actin bundles during cell movement[43]. A recent study that investigated alterations in brain transcriptomics associated with intergenerational stress transmission reported upregulation of BAIAP2L1 in neonatal and adult mice[44]. In addition, methylation levels of cg00277769 in BAIAP2L1 correlated in blood and brain tissues of human subjects as reported in IMAGE-CpG[31]. Hence, BAIAP2L1 expression might be regulated in response to stress and trauma.

Of particular interest in relation to PTSD are the findings relating to TBXAS1. The ER membrane protein TBXAS1 metabolizes Prostaglandin H2 (PGH2), which regulates the dilation of blood vessels[45] to i) Thromboxane A2 (TXA2), which is critical during inflammation[46, 47], ii) 12-Hydroxyheptadecatrienoic acid (12-HHT), which may participate in monocyte- and neutrophil-based inflammation[48-50], and iii) Malonyldialdehyde, a marker for oxidative stress[51]. This metabolic reaction may be regulated by DNAm, as a study conducted in endothelial cells reported that TBXAS1 demethylation resulted increased thromboxane B2 (TXB2), the product of TXA2 breakdown[52]. TXA2 promotes platelet aggregation by binding to thromboxane receptor (TP)[53]. In addition, TXA2–TP signaling has been suggested to amplify dopamine overflow from the striatum[54]. Altered striatal dopamine function has been linked to early life and adulthood adversity, as well as psychiatric disorders, such as schizophrenia and PTSD[55-59]. In addition, SNPs in TBXAS1 have been reported to predict gray matter volume changes of left collateral sulcus of visual cortex (hOC3vL) in schizophrenia[60]. Interestingly, the SNP in TBXAS1 that associates with cg03604364 methylation levels (rs3779130) was nominally associated with PTSS in the recent PGC-PTSD meta-analysis (z = 2.00, p = 0.045)[61]. According to the BIOS meQTL browser and PGC-PTSD meta-analysis, carriers of the rs3779130 T allele have lower methylation levels and higher PTSS. This aligns with our own findings that individuals with higher PTSS have lower cg03604364 methylation levels. Moreover, cg03604364 showed the strongest within-individual correlation of methylation measured in blood and brain among the four CpGs associated with PTSS. Together, these results suggest that alterations within TBXAS1 associate with traumatic events and could contribute to psychiatric disorders, possibly though inflammatory and/or oxidative stress pathways.

In the second meta-analysis, no CpGs associated with the change in PTSS. Nevertheless, in our region-based analyses, we identified 15 DMRs associated with change in PTSS. For all 15 DMRs, direction of effect was positive, indicating increased DNAm with increased post-deployment PTSS, conditional on baseline DNAm. Interestingly, three DMRs were located on the X chromosome, which is largely overlooked or excluded in EWAS due to sex-dimorphic distribution[62]. Since this study only included males, we could interpret X chromosome results. The three X chromosome DMRs are located on ovarian tumor deubiquitinase 5 (OTUD5), 74 like ETS transcription factor 4 (ELF4) and RB binding protein 7, chromatin remodeling factor (RBBP7). Notably, OTUD5, ELF4 and RBBP7 expression has previously been implicated in PTSS development in women upon trauma[63]. DMRs near OTUD5 and ELF4 were also associated with PTSS in our first DMR analysis that uses the same framework as Meta-Analysis 1. The other notable DMR findings are adenylate cyclase 6 (ADCY6) and the GNAS Complex Locus (GNAS), which have previously been implicated in PTSD. ADCY6 expression was altered in the amygdala of a PTSD-like mouse model[64]. The GNAS locus is known for its complex imprinted expression pattern and produces multiple transcripts due to promoter DMRs[65]. Differential methylation of GNAS has been shown to be associated with schizophrenia in multiple studies[66-68]. In addition, GNAS was differentially expressed in blood, hemibrain and spleen of a PTSD mouse model[69].

Finally, we assessed whether 31 CpGs from previous PGC-PTSD EWAS results[10-14] were also significant in our study. Only three (cg19534438 in G0S2, cg05575921 and cg26703534 in AHRR) were nominally significant with the same direction in Meta-Analysis 1. The inconsistency of the remaining CpGs might be due to previous overestimation of the effect, or to the operationalization of the PTSD phenotype, as all previous studies used dichotomous case-control phenotype, whereas we used symptom scores. We did not attempt to evaluate DMR findings due to methodological challenges related to the DMR analysis, that is, DMR results may vary by methylation array type and the DMR analysis methodology.

Our study is not without limitations. First, the MethylationEPIC BeadChip captures only ~3% of the human methylome and, thus, we do not interrogate all CpGs that may associate with PTSS or change in PTSS[70]. Second, in Meta-Analysis 1, we used a linear mixed model to identify CpGs that are associated with PTSS. Hence, that analysis does not distinguish whether the identified epigenetic differences are a cause or consequence of the symptomatology. Third, the methylation data was generated from blood. Though this approach captures PTSD-related DNAm changes and is informative to discover biomarkers of PTSD, it may not reflect DNAm patterns in brain. However, two out of four significant CpGs are moderately correlated between blood and brain tissues, suggesting a similar PTSD-related DNAm pattern for these CpGs in brain. Finally, the cohorts participating in this study are comprised of male and predominantly European participants who were exposed to military trauma. Hence, it is not clear how these findings will translate to females, civilians, or other ancestry groups.

Despite these limitations, this work represents the largest longitudinal study of epigenetics in relation to traumatic stress to date. Our results support a role for epigenetic mechanisms in PTSD severity and implicate genes that are involved in immune system and oxidative stress, and they align with previous studies that support the role of inflammatory processes in PTSD through HPA-axis reactivity, co-morbid metabolic conditions, and behavioral changes common in individuals with PTSD[71]. Our findings also support the need to fully understand the regulation of biologically significant genes, such as CNPY2, BAIAP2L1, TBXAS1, ADCY6 and GNAS, with regards to PTSD, particularly though functional studies that can delineate directionality in PTSD development, symptom severity, or treatment.