Epigenome-wide meta-analysis of PTSD across 10 military and civilian cohorts identifies methylation changes in AHRR

Epigenetic differences may help to distinguish between PTSD cases and trauma-exposed controls. Here, we describe the results of the largest DNA methylation meta-analysis of PTSD to date. Ten cohorts, military and civilian, contribute blood-derived DNA methylation data from 1,896 PTSD cases and trauma-exposed controls. Four CpG sites within the aryl-hydrocarbon receptor repressor (AHRR) associate with PTSD after adjustment for multiple comparisons, with lower DNA methylation in PTSD cases relative to controls. Although AHRR methylation is known to associate with smoking, the AHRR association with PTSD is most pronounced in non-smokers, suggesting the result was independent of smoking status. Evaluation of metabolomics data reveals that AHRR methylation associated with kynurenine levels, which are lower among subjects with PTSD. This study supports epigenetic differences in those with PTSD and suggests a role for decreased kynurenine as a contributor to immune dysregulation in PTSD.

P ost-traumatic stress disorder (PTSD) is characterized by reexperiencing, avoidance, and hyperarousal symptoms that can negatively impact mood and physiologic health 1 .
Although not everyone who experiences trauma goes on to develop PTSD, those who do often experience severe and disabling symptomatology 2,3 . Research suggests that both genetic and environmental factors contribute to risk for developing PTSD 4,5 .
Given the necessary, but not sufficient, role of environmental exposure (i.e., trauma) in the development of PTSD, it is critical to characterize the pathways underlying differential risk and resilience. Studies that contexualize the role of environmental influences provide additional insight into modifiable factors that may promote post-trauma resilience. For example, lack of social support at the time of trauma is associated with increased risk of developing PTSD in both military and civilian contexts 6 . Similarly, studies that investigate how presumed environmental influences might affect biological pathways could provide insights into the genes whose regulation patterns differ in those with PTSD. A growing line of research aimed at elucidation of mechanisms by which environmental factors contribute to PTSD, has focused on epigenetic modifications, which regulate gene function in response to environmental triggers. Epigenetic modifications, such as DNA methylation at cytosine-guanine dinucleotides (CpG sites), can induce changes in gene expression that are maintained through each round of cell division.
Multiple reviews have linked traumatic stress to differences in the proportion of methylated DNA intensity to non-methylated DNA intensity (ß) at specific CpG sites [7][8][9][10] . Indeed, a number of epigenome-wide association studies (EWAS) of individual cohorts have identified PTSD-associated CpG sites in genes and pathways related to neurotransmission and immune function [11][12][13][14][15] . Similarly, studies using methylation of specific CpG sites that predict chronological age with reasonably high accuracy to capture age acceleration (deviations from a line of prefect prediction), demonstrate associations with PTSD and link differences in peripheral DNA methylation to memory formation and neural integrity 16,17 . Although promising, the extant literature of EWAS studies on PTSD is limited by use of individual cohorts with small sample sizes and varying analytic pipelines that can make it challenging to synthesize findings. Consortia efforts can overcome these limitations by providing a shared analytic pipeline and increasing sample size and thus statistical power. The goal of this study is to capitalize on consortium strengths by conducting a meta-analysis of DNA methylation across 10 military and civilian cohorts participating in the Psychiatric Genomics Consortium (PGC) PTSD Epigenetics Workgroup. Our results suggest that lower aryl-hydrocarbon receptor repressor (AHRR) methylation in those with PTSD correlates with lower kynerunine levels, which may contribute to immune dysregulation in PTSD.

Results
Participating cohorts. Sample characteristics for the 10 cohorts that have contributed data are listed in Table 1 (N = 1896). All participants were exposed to trauma, and 42% had a current diagnosis of PTSD. There were no significant age differences in PTSD cases and trauma-exposed controls. However, the demographic characteristics for each cohort varied substantially (Supplementary Data 1), with cases more likely to be male (X 2 = 22.4; p < 0.05) and current smokers (X 2 = 54.2; p < 0.05) across the majority of cohorts and in the overall cohort.
PTSD-associated CpG sites from Meta-Analysis. In our primary meta-analysis, we found 4 CpG sites associated with current PTSD after correction for multiple comparisons (−6.60 < z < −5.57; p < 3.6E-08; Fig. 1; Table 2; Supplementary Data 2). AHRR contains the top 4 PTSD-associated CpGs, with lower methylation in PTSD cases relative to controls ( Supplementary Fig. 1). This was consistent when stratified by both sex (Supplementary Fig. 2) and ancestry (Supplementary Fig. 3). Data from iMETHYL, which reports eQTMs associated at false-discovery rate (FDR) < 0.05, suggests that methylation of both cg05575921 and cg25648203 are inversely associated with AHRR expression in peripheral blood mononucular cells (PBMCs).
Sensitivity analysis with smoking status. As lower methylation of AHRR CpG sites has been associated with smoking 18-20 , we controlled for smoking status in our sensitivity analyses of the 4 significant CpGs (Supplementary Fig. 4). The association between AHRR methylation and PTSD was attenuated for all four CpGs. Since there is a higher rate of smoking among PTSD cases, we evaluated potential differences of effect between PTSD and smoking status by testing the associated CpGs separately in smokers and non-smokers (Fig. 2) and by plotting the associated CpGs separtely by smoking status (Supplementary Fig. 5). The associations between AHRR CpGs and PTSD were most prominent in non-smokers. Of the multiple stratified analyses performed (Supplementary Figs. 2, 3, and 5), the effect sizes were consistent across different strata, with the smokers being the only one did not show an association with PTSD.
We further delineated the association between PTSD and smoking-associated CpG sites from a large meta-analysis on smoking conducted by the CHARGE consortium 20 . We hypothesized that unreported smoking among those that identify as non-smokers would result in the appearance of an association between smokingassociated CpG sites and PTSD. This evaluation included 21 CpGs representing the most significantly associated CHARGE smoking sites and CpGs in AHRR. If the association between AHRR and PTSD was due to unreported smoking in PTSD controls, we would expect to see a similar effect across all smoking-CpGs. However, after adjusting for the number of tests performed, AHRR CpG sites remained associated with PTSD in the non-smokers only, and no smoking CpG showed a significant association with PTSD in the controls or cases. These findings are consistent with the hypothesis that the association between PTSD and DNA methylation of AHRR is independent of smoking status ( Supplementary Fig. 6).
Tryptophan catabolism in PTSD. We next evaluated other biological factors that could contribute to aryl-hydrocarbon receptor (AHR) expression using linear regression models. AHR is activated in many immune cell types, including T cells, B cells, and NK cells, by indolamine-mediated tryptophan catabolism 21 . Inflammation can shift tryptophan metabolism away from serotonin synthesis towards kynurenine synthesis 22 . To evaluate the association of tryptophan catabolism in PTSD and its relationship to AHRR methylation, we leveraged data from the 116 subjects from MRS that had both DNA methylation and tryptophan metabolite data. In this group, we first note that kynurenine levels were lower in subjects with PTSD relative to controls (t = −2.00; p = 0.048; Supplementary Fig. 7). Consistent with this observation, lower methylation of AHRR CpGs associated with lower kynurenine (cg21161138; t = 2.09, p = 0.039) and lower kynurenic acid (cg05575921; t = 2.95, p = 0.004).
Further, metabolome data in the MRS confirmed that selfreported smoking status is consistent with empirical cotinine levels; 87% of self-reported non-smokers and light smokers had cotinine levels consistent with no-smoking, second-hand smoking or light exposure, while 91% of regular or heavy smokers had cotinine levels consistent with their endorsement. In addition, these data also provide insight into why controlling for smoking may attenuate the association between DNA methylation and PTSD. We noted an inverse relationship between cotinine and kynurenine (r = −0.267; p = 0.0004) that was consistent among both smokers and non-smokers ( Supplementary Fig. 8).

Discussion
An individual's risk of developing PTSD depends on both the nature of the trauma and the physiological response to that trauma 23 . Not all individuals exposed to trauma develop PTSD, and a better understanding of the modifiable biological factors underlying risk and resilience will inform the development of new prevention and treatment strategies. In this study we identified CpG sites associated with PTSD. We observed that, on average, PTSD cases had lower methylation at several CpG sites in the AHRR gene when compared to trauma-exposed controls. Methylation of AHRR CpGs has been strongly linked to smoking [18][19][20] . As substantially more of the PTSD cases reported smoking compared to controls, this suggested that we should control for smoking in the EWAS. A comparable approach was taken by Marzi and colleagues 24 . In their study of DNA methylation in relation to victimization stress in adolescents, they noted experiment-wide significant associations in 3 of the 4 AHRR CpGs associated with PTSD. Similar to our findings,    the association of AHRR and PTSD was no longer significant after statistically controlling for smoking. Marzi and colleagues concluded that there are no robust changes in DNA methylation in victimized young people. Our study went beyond this approach to conduct a stratified analysis that revealed the most prominent PTSD-associated difference in AHRR methylation was evident in the non-smokers 24 . To evaluate the possibility of confounding due to unreported smoking among those that identify as nonsmokers, we tested top smoking-related CpGs, including AHRR, for association with PTSD, and identified no evidence of association between smoking-related CpG sites and PTSD, supporting our conclusion that the association between AHRR CpG sites and PTSD was independent of smoking.
Though it is possible that stratifying on smoking could introduce a collider bias that distorts the stratified analysis because of a SNP or other unmeasured confounder, this scenario seems unlikely given that the associations between AHRR CpGs and PTSD are stronger in non-smokers. As such, this cross-sectional analysis is unable to determine whether methylation of AHRR CpGs is likely a cause or consequence of smoking. However, two recently published studies have provided insight into this question. The first sought to identify DNA methylation differences associated with PTSD in blood of 378 PTSD cases and 135 controls from the Translational Research Center for Traumatic Brain Injury and Stress Disorders (TRACTS) cohort 25 . This study used a methylation-based score for smoking as a covariate, and still reported a significant association between PTSD and cg05575921 (AHRR). The second examined DNA methylation differences in blood associated with polycyclic aromatic hydrocarbon (PAH) levels unrelated to smoking 26 . In 708 non-smokers from Taiwan, the authors report associations between PM 2.5 levels, an indicator of regional air quality, and lower methylation at cg05575921. Taken together, these studies suggest that methylation of this CpG can occur through processes other than smoking.
Over the last decade, multiple studies described the role of the aryl-hydrocarbon receptor in the immune system, with specific roles in T cells, B cells, monocytes and dendritic cells (DC) 21,[27][28][29] . The aryl-hydrocarbon receptor (AhR) pathway works in a regulatory capacity. Briefly, when a ligand binds to the AhR, it translocates to the nucleaus where it drives expression of its target genes, including the aryl-hydrocarbon receptor repressor (AHRR), which begins a feedback loop in which it can competitively bind the AhR. The AhR pathway can either limit or stimulate an inflammatory response. One mechanism by which this occurs is by promoting differentiation of T cells into T regulatory cells (Tregs) or T helper 17 (Th17), though this appears to be done in a ligand-specific manner. For example, endogenous ligands, such as kynurenine or dietary indoles 30 , tend to promote differentiation into Tregs, which reduce the immune response in a self-limiting manner, resulting in expression of anti-inflammatory genes, reduced inflammation, and less Treg generation. In contrast, exogenous ligands, such as dioxin or polyaromatic hydrocarbons in cigarette smoke, tend to promote differentiation into Th17, expression of pro-inflammatory cytokines, heightened inflammation, and activation the drug metabolizing enzymes 27,31,32 .
In a subset of subjects with tryptophan metabolite data, we reported that kynurenine levels were significantly lower in subjects with PTSD relative to trauma-exposed controls and that lower methylation of AHRR CpGs associated with lower kynurenine and its metabolites. This pattern is similar to what is observed following chronic exposure to nicotine and possibly other AHR-stimulating ligands 33 . This chronic exposure scenario is reflected in our cohort among both controls and PTSD cases by the finding that higher levels of cotinine were associated with decreased levels of kynurenine. Independent of its source, lower levels of kynurenine would likely result in a counter-regulatory increase in pro-inflammatory activity 34 , and may help to account for the frequent observation of heightened inflammatory activity observed in subjects with PTSD 35 . Our secondary analysis of DNA methylation and kynurenine was limited to the only cohort with cotinine and metabolomics data available on the subjects included in the methylation study, which limited sample size, and power for this analysis, and potentially generalizability to other cohorts. However, we believe that the ability to verify selfreported smoking status with a biological measure and to evaluate our kynurenine hypothesis was a considerable strength of the study. We hope these results will prompt other investigators to replicate and extend these findings.
While smokers are more common among our PTSD group, approximately half of the cases in our cohorts are non-smokers, who exhibited the most prominent associations between AHRR CpGs and PTSD. Methylation of AHRR CpGs decreases only so much even among heavy smokers, suggesting that there may be limited variability in these subjects and a limited degree to which methylation at this site may be decreased. Though the vast majority of the AHRR methylation literature is focused on characterizing its variation in the context of smoking, the types of polycyclic aromatic hydrocarbons that stimulate AHR are present from multiple exogenous sources, including burning wood or charcoal, auto emissions, industrial exhaust, and urban dust 36 . Similarly, participants in the WTC and military cohorts are likely to have substantial levels of occupational exposure that could result in AHR activation 37-39 . Though we were not able to directly assess exposure to endogenous or exogenous AhR ligands in this study, it is reasonable to hypothesize that smoking is just one of many potential environmental exposures that could promote inflammation following AHR stimulation.
Both epidemiologic and immunologic studies report immune dysregulation in those with PTSD compared to controls. For example, autoimmune and inflammatory disorders, such as rheumatoid arthritis, have been linked to PTSD [40][41][42] . Studies of the immune system generally support higher levels of pro-inflammatory cytokines in PTSD cases relative to controls 35 . A study investigating Tregs reported a lower number of Tregs in the blood of PTSD patients 43 , while another reported a difference in the composition of Treg populations that suggested higher susceptibility for autoimmunity 44 . Another study reported lower proportions of Tregs and higher proportions of Th17 cells in PTSD cases before linking Th17 counts to higher clinical symptom scores 45 . Finally, a randomized control trial of subjects undergoing narrative exposure psychotherapy noted higher levels of Tregs and lower PTSD symptoms following treatment 46 . Unfortunately, we were unable to directly measure or impute the proportion of different types of T cells in this study. Though the degree to which the immune system is involved remains speculative without functional studies, we hope these data will support more detailed immunophenotyping of subjects with PTSD. This study has a number of additional limitations that should be considered. First, we used existing data generated on the HumanMethylation450 array, which captures only a fraction of the CpG sites in the genome. Though sequencing methods would capture a larger proportion of the epigenome that may include regions important for PTSD, focusing on this array allowed us to evaluate a larger and more diverse cohort of subjects and enabled processing and analysis of these diverse samples with a common analytic pipeline. Similarly, this meta-analysis uses only crosssectional data. In order to establish whether these PTSDassociated differences are a cause or consequence of the disorder or both, longitudinal and functional studies will be required. Second, there is phenotypic heterogenetity among the cohorts contributing to this meta-analysis. Some cohorts assess individuals that were recently exposed to trauma while others include individuals with chronic PTSD. Also, few of the cohorts included in this meta-analysis have detailed physical or psychiatric information on subjects prior to trauma exposure, making it difficult to evaluate the role of lifestyle factors, such as obesity, or comorbidities, such as substance use. Thus, it is possible that some of the epigenetic differences observed in this study may have been in place prior to PTSD development, including the possibility of previous episode of PTSD in subjects that no longer meet current diagnostic criteria, or reflect other factors such as differences in genetic background between cases and controls. Third, our study examined whole blood-derived DNA. Though this approach likely captures part of the PTSD sequelae and may be informative for future biomarker studies, it is unlikely to reflect DNA methylation in brain regions most relevant for PTSD. As studies of tissues from PTSD Brain Banks are conducted, it will be important to look for parallels between these different tissues.
Taken together, the results of this study implicate the immune system in PTSD and suggest that epigenetic mechanisms may play a role in that process. A substantial fraction of those diagnosed with PTSD do not respond to pharmacologic or psychological interventions, and clinical and preclinical studies have begun to evaluate strategies to limit inflammation as a first line treatment 47,48 . Future studies should evaluate the role of ligandspecific AHR activation in the development and progression of PTSD.

Methods
Post-traumatic stress disorder cohorts and assessments. The participating cohorts, presented in Table 1 cohort from the National Center for PTSD (VA-NCPTSD). All subjects participating in these studies provided informed consent, and all studies were approved by respective institutional review boards. Current PTSD diagnosis was assessed by each individual cohort in accordance with the harmonization principles adopted by the PGC-PTSD Workgroup 4 . Briefly, diagnosis of current PTSD was based on the diagnostic criteria defined by each cohort's principal investigator (see cohort descriptions in Supplementary Methods for complete details). All control subjects were trauma-exposed, and, if assessed in the respective cohort, control subjects that had a prior history of PTSD were excluded. Covariates were age, sex, genetic ancestry, and smoking status. A total of 1,896 subjects (42% cases) with DNA methylation from whole blood measured using the Illumina HumanMethyla-tion450 BeadChip were selected for inclusion in this meta-analysis.
Quality control (QC) procedures. Data contributed to this study were part of 10 different studies, which were each designed to address a study-specific question. In developing the pipeline to analyze this existing data, the PGC-PTSD group issued recommendations for studies to balance potential confounders. All cohorts were instructed on importance of balancing plate layouts, and the potential for residual sources of technical variation was further assessed by examining the association of methylation principal components with chip, position, plate, and demographic characteristics within each study. Each study performed analyses at their site. To ensure consistent QC procedures across all participating cohorts, a set of common scripts were developed and implemented uniformly 49 . ß-values, representative of the proportion of methylation at each probe, were caculated for all CpG sites in each sample. Samples with probe detection call rates <90% and those with an average intensity value of either <50% of the experiment-wide sample mean or <2000 arbitrary units (AU) were excluded. Probes with detection p-values > 0.001 or those based on less than three beads were set to missing as were probes that cross-hybridized between autosomes and sex chromosomes 50 . CpG sites with missing data for >10% of samples within cohorts were excluded from analysis. Normalization of probe distribution was conducted using Beta Mixture Quantile Normalization (BMIQ) 50 after background correction. Density plots of Type I and Type II probes before and after normalization were examined to confirm probe distributions between types were similar. ComBat was used to account for sources of technical variation including batch and positional effects, while preserving variation attributable to study-specific outcomes and covariates that would be used in downstream analyses (e.g., case status or sex). Proportions of CD8, CD4, NK, B cells, monocytes, and granulocytes were estimated using each individual's DNA methylation data based on the approach described by Houseman and colleagues and publicly available reference data (GSE36069) 51,52 . For cohorts without GWAS data (DNHS, PRISMO, WTC), the method described by Barfield and colleagues was used to generate ancestry PCs directly from the methylation data using CpG sites from the HumanMethyla-tion450 array that are located within 0-100 bp of 1000 Genomes Project variants with minor allele frequency > 0.01, and PCs 2-4, which were the components that correlate most with ancestry, were used to adjust for ancestry as recommended by the authors 53 . Using the empirical Bayes method in the R package limma 54 , moderated t-statistics were calculated for each CpG site, converted first into onesided p-values then converted into z-scores to account for direction of effects. Meta-analysis across cohorts was performed by weighting each cohorts' z-scores by its sample size relative to the total meta-analysis sample, summing weighted zscores across cohorts from which two-sided p-values were calculated. We used the genome-wide significance threshold suggested for HumanMethylation450 Bead-Chip to adjust for multiple testing (p < 3.6E-8) 55 . The iMETHYL database (http:// imethyl.iwate-megabank.org/) was used to evaluate whether methylation of AHRR CpGs associated with AHRR expression in PBMCs, the available tissue that most closely resembles whole blood 56 .
Metabolite analysis in the Marine Resiliency Study (MRS). Targeted, broadspectrum metabolomics was performed as previously described 57,58 with minor modifications. Lithium heparin plasma samples were collected and stored at −80°C until used for analyses in the MRS metabolomics study. Samples were analyzed from 116 participants exposed to military combat; 53 diagnosed with PTSD, and 63 controls without. Ninety microliters of plasma and 10 µl of added stable isotope internal standards were combined, extracted, and analyzed by hydrophilic interaction chromatography, electrospray ionization, and tandem mass spectrometry on a SCIEX 5500 QTRAP HILIC-ESI-MS/MS platform. From the 486 metabolites measured, we extracted and analyzed data for kynurenine and kynurenic acid to test our hypothesis that AHRR methylation associates with trypophan breakdown products. Cotinine, a metabolite of nicotine with a plasma half-life of about 18 h, was measured to provide independent information about tobacco exposure. A cotinine arbitrary unit (AU) of 1 × 10 6 was used as the upper NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19615-x ARTICLE NATURE COMMUNICATIONS | (2020) 11:5965 | https://doi.org/10.1038/s41467-020-19615-x | www.nature.com/naturecommunications limit for a non-smoker, light smoker, or second-hand tobacco exposure, with ≥1 × 10 6 indicative of regular or heavy smoking. Analyte AU data was log-transformed, and analyte z-scores were used for further statistical analysis. To determine if trytophan metabolite levels differed between PTSD cases and controls, we performed linear regressions of metabolites on PTSD status, including GWAS PCs 1-3 and cell type proportions as covariates. To determine if tryptophan metabolite levels were associated with AHRR methylation, we performed linear regressions of metabolites on each AHRR CpG site, including PCs and cell type proportions as covariates. Correlation between cotinine and kynurenine was measured using Pearson's correlation coefficient.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The main summary statistics data that support the findings of this study are available within Supplementary Data 2. Individual-level data from the cohorts or cohort-level summary statistics will be made available to researchers following an approved analysis proposal through the PGC Post-traumatic Stress Disorder group with agreement of the cohort PIs. The raw data for the GTP cohort is available in the Gene Expression Omnibus database with the accession code GSE72680. Owing to military cohort data sharing restrictions, data from the VA, VA MIRECC, MRS, Army STARRS, and PRISMO cannot be publicly posted. However, such data can be provided in de-identified from a data repository through a data use agreement following applicable guidelines on data sharing and privacy protection. For additional information on access to these data, including PI contact information for the contributing cohorts, please contact the corresponding author.