Genetic variant in CACNA1C is associated with PTSD in traumatized police officers.

Posttraumatic stress disorder (PTSD) is a debilitating psychiatric disorder that may develop after a traumatic event. Here we aimed to identify epigenetic and genetic loci associated with PTSD. We included 73 traumatized police officers with extreme phenotypes regarding symptom severity despite similar trauma history: n = 34 had PTSD and n = 39 had minimal PTSD symptoms. Epigenetic and genetic profiles were based on the Illumina HumanMethylation450 BeadChip. We searched for differentially methylated probes (DMPs) and differentially methylated regions (DMRs). For genetic associations we analyzed the CpG-SNPs present on the array. We detected no genome-wide significant DMPs and we did not replicate previously reported DMPs associated with PTSD. However, GSE analysis of the top 100 DMPs showed enrichment of three genes involved in the dopaminergic neurogenesis pathway. Furthermore, we observed a suggestive association of one relatively large DMR between patients and controls, which was located at the PAX8 gene and previously associated with other psychiatric disorders. Finally, we validated five PTSD-associated CpG-SNPs identified with the array using sanger sequencing. We subsequently replicated the association of one common SNP (rs1990322) in the CACNA1C locus with PTSD in an independent cohort of traumatized children. The CACNA1C locus was previously associated with other psychiatric disorders, but not yet with PTSD. Thus, despite the small sample size, inclusion of extreme symptom severity phenotypes in a highly homogenous traumatized cohort enabled detection of epigenetic and genetic loci associated with PTSD. Moreover, here we showed that genetically confounded 450K probes are informative for genetic association analysis.

Introduction

Posttraumatic stress disorder (PTSD) is a debilitating psychiatric disorder that may develop after exposure to traumatic events. The disorder is characterized by recurrent, involuntary re-experiencing of the traumatic event, accompanied by avoidance of trauma-related stimuli; negative alterations in cognition and mood and in reactivity [1]. Its lifetime prevalence in Western countries is estimated at 5.7–7.4% [2, 3]. In all, the societal impact of PTSD is high, not only because of its negative impact on quality of life of affected individuals and their families, but also due to high economic costs resulting from health care utilization and productivity loss and the high percentage of patients not adequately responding to treatment [4, 5]. Additionally, PTSD is associated with increased risk for a wide range of severe chronic medical disorders and increased all-cause early mortality [6-8]. Longitudinal studies consistently show that only a subset of individuals will develop PTSD in the aftermath of apparently similar traumatic events, indicating differential susceptibility [9]. To date, the biological mechanisms and etiology underlying differential susceptibility for PTSD remain poorly understood, but accumulating evidence indicates that genetic and epigenetic factors are involved [10, 11].

Genetic factors contribute substantially to PTSD development: its heritability is estimated between 0.28 and 0.46 [12-14]. PTSD has a highly heterogeneous phenotype and can be characterized as complex trait, indicating that common genetic variation within hundreds to thousands of genes might contribute to its expression. For most of these genetic contributors, small effect sizes can be expected [15]. Moreover, PTSD development cannot be explained by genetic factors alone [12, 16, 17]. PTSD involves by definition exposure to a precipitating traumatic event suggesting an additional role for epigenetic processes that mediate gene regulation under the influence of the environment [18]. Recent insights show that epigenetic associations with psychiatric disorders can not only be detected in the most obvious target tissue, i.e. the brain, but also in peripheral blood [19]. With regard to PTSD, several candidate gene studies have detected genetic and epigenetic associations. However, the number of hypothesis-free, i.e. genome wide genetic and epigenetic studies on PTSD, remains limited [20, 21]. In this study, we aimed to identify epigenetic and genetic loci associated with PTSD in a Dutch cohort of highly trauma exposed police officers with and without PTSD. We selected our cohort to be homogenous, [also] with respect to trauma history, but showing extreme phenotypes regarding PTSD symptom severity (i.e. highly affected vs. resilient with minimal symptoms), maximizing likelihood to detect PTSD-related associations in a relatively small sample [22]. We used the Infinium HumanMethylation 450 BeadChip of Illumina, a commonly used platform for epigenome-wide association studies. Next to ~430K informative epigenetic markers, this array contains a substantial number of probes biased by common genetic variation. These probes involve common single nucleotide polymorphisms at the CpG site of interest (CpG-SNPs) directly affecting the readout of the methylation array at that position. Usually these probes are discarded from epigenetic surveys, despite the fact that some of the ~14K probes show well-defined genotypic clusters that may be used for genetic association analysis. Detailed information on promiscuous 450K probes, i.e. genetic confounded and cross hybridizing probes, and the design of the array is described by Chen et al. [23]. In this study, we used these 14K CpG-SNP probes for genetic association analyses and subsequently validated identified candidate PTSD-associated SNPs with additional sanger-sequencing. Consistent genetic findings between array and sequencing methods were replicated in a Dutch cohort of trauma exposed children with and without PTSD.

Subjects and methods

Participants

This study is part of a larger neuroimaging study on the effects of a single oxytocin administration on neural activity, in which trauma-exposed police officers with and without PTSD between 18 and 65 years old were included [24]. For the current study, blood samples were collected from 73 participants: 34 PTSD patients and 39 controls. PTSD patients fulfilled DSM-IV criteria for current PTSD, with a Clinician Administered PTSD Scale (CAPS) score ≥ 45 [25]. Exclusion criteria for PTSD patients were current severe MDD (with psychotic features or suicidal intent), bipolar disorder, suicidal ideation, alcohol/substance abuse, and psychotic disorders, as measured with the Mini International Neuropsychiatric Interview or Structured Clinical Interview for DSM–IV [26, 27]. Controls reported at least one traumatic event (DSM-IV A1 criterion), with a CAPS score <15. Controls were matched to PTSD patients based on sex, age, education, and years of service. Exclusion criteria for controls were lifetime MDD or PTSD or any current DSM-IV Axis-I psychiatric disorder. Exclusion criteria for all participants consisted of daily use of psychoactive medication or systemic glucocorticoids, contraindications for MRI or oxytocin administration, serious medical conditions or history of neurological disorders. All participants provided written informed consent prior to study initiation. The study was approved by the Institutional Review Board of the Academic Medical Center, Amsterdam, the Netherlands. Patients were recruited through a diagnostic outpatient center for police personnel with psychological trauma (PDC, Diemen, the Netherlands). Additionally, PTSD patients and controls were recruited through advertisements. Group differences in demographic and clinical characteristics were tested using independent sample t-tests (log transformed in case of non-normally distributed data) for continuous variables and Chi-square tests or Fisher’s exact tests for categorical variables, performed in SPSS (version 20.0.0, IBM® SPSS® Statistics).

Procedures

Sample preparation and quality control of the 450K DNA-methylation array are described in appendix 1. The 450K DNA-methylation array comprises a substantial number of SNPs located in or near CpG sites at single base extensions (SBE) sites of probes. As was previously shown by others [23], a substantial number of these CpG-SNPs show 3 or 2 well divided clusters of beta-values, reflecting genotypes. Therefore, we divided the probes into two data sets, one for epigenetic and one for genetic association analysis (see Figure S1 for schematic overview). For both data sets promiscuous 450K probes, i.e. cross hybridizing probes, were excluded from analysis [23]. In addition, probes located on both the X and Y chromosomes were not included in further analysis. The epigenetic probe set included 432.222 probes, excluding probes with a common SNP with a CEU minor allele frequency (maf) > 0.01 located at the SBE site. This probe set was used to identify differentially methylated probes (DMPs) and differentially methylated regions (DMRs). The CpG-SNP genetic probe set included the 13.344 probes with a common SNP located at SBE site with a CEU maf > 0.01. Unprocessed genome wide epigenetic and meta data are freely available NCBI Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) Accession ID:

Epigenetic association analyses

A complete description of DMP and DMR analyses is provided in appendix 1. For both DMP and DMR analyses, we included the following covariates in the model: age, leukocyte cell distribution (CD8+ and CD4+ T cells, Natural killer cells, B cells, Monocytes, and Granulocytes), and sex. For DMP association analyses lmfit (“R”package limma) function was used. We assumed a q-value < 0.05 (Benjamini–Hochberg adjusted alpha) as genome-wide significant. In order to replicate previously reported epigenetic associations (DMPs) with PTSD we additionally applied a hypothesis-driven approach; we performed a NCBI PubMed search using keywords “450K”, “Illumina” and “PTSD”. Heterogeneity between studies complicates interpretation and power of the replication analyses. We therefore applied strict inclusion criteria on PTSD characterization, technical platform, abuse history, and mean age of the sample and the type of tissue of which DNA was obtained. From the identified studies we selected the candidate DMPs that were previously reported to be genome-wide significantly associated with PTSD, and for these DMPs data was extracted from our association data. Subsequently, we performed gene-set enrichment (GSE) analyses on our top 100 DMP hits, excluding DMPs located at intergenic regions (IGRs, >1 kb distance from the gene, 26 DMPs removed, 74 DMPs remained). Enrichment analyses were done via the publicly available consensus pathway database web tool, developed by the Max Planck institute, version 31, 2015 [28].

DMR association analyses were performed using the minfi function bumphunter. DMRs were defined as including two or more probes (L > 1) in the bump. Genome wide significant DMRs (adjusted P-value < 0.05 according the familywise error-rate over the area of the bump, FWER) or DMRs showing consistent differential methylation over more than 75% of the probes within the cluster were considered promising. Technical validation was performed targeting one promising DMR (PAX8) and one randomly chosen region located at the SOSTDC1 gene. The latter region was included to achieve a broader range of β-values that favors proper correlation analysis between the two methods. SOSTDC1 expressed average β-values < 50% within that region. Technical validation was performed through targeted amplicon sequence analysis using the Illumina MiSeq platform (details see appendix 1).

Genetic association analyses

In order to detect genetic associations between PTSD and the 450K CpG-SNPs, we performed association analyses that conform the above described statistical methodology (epigenetic association analyses). We included only CpG-SNPs located at the SBE site or subsequent guanine base with a maf <0.01. Secondly, the dataset was filtered on CpG-SNPs showing a delta mean beta difference exceeding 15% and a nominal P-value < 0.05, which yielded 15 CpG-SNPs showing well divided clusters of beta-values (representing the genotypes). Since our study is limited in statistical power, only five CpG-SNPs, showing the lowest nominal P-value and most well divided clusters, were followed up by a technical validation using sanger sequencing. Validation (and replication, described below) P-values < 0.01 were assumed significant according a Bonferroni correction, based on five independent tests. Since the amount of collected DNA was limited, we first applied a whole genome amplification for all available samples using the Qiagen REPLI-g Mini Kit™ conform the manufacturers protocol, yielding >4 μg of DNA for each participant. Primer sets for candidate SNPs were designed using Primer3 [29], all primers were M13 extended for standardized Sanger sequencing (Table S1). For validation association analyses genotypes were recoded into ordinal variables followed by Chi-square tests assessing genotype differences between PTSD patients and controls, performed in SPSS (version 20.0.0, IBM® SPSS® Statistics). P-values at which significance was assumed were Bonferroni corrected for the number of tests performed (0.05/number of tests).

Replication analyses of validated CpG-SNPs were performed in a cohort of trauma-exposed children with (N = 46) and without PTSD (N = 50), see appendix 1 for details on participant characteristics and genotyping procedure. Epigenetic profiles may differ in a greater or lesser extent between tissues. Since DNA of replication cohort was extracted from saliva and DNA of the discovery cohort was obtained from whole blood we considered only replication of genetic associations in the discovery cohort appropriate. Statistical analyses were performed as described above for validation analyses. P-values at which significance was assumed were Bonferonni corrected for the number of tests performed (0.05/number of tests). Genotypes and metadata of the five reported SNPs for all individuals of the replication cohort are available in Table S3—tab Genotypes_Replicationcohort

Results

Discovery cohort participant characteristics

PTSD patients and trauma-exposed controls did not differ on demographic or trauma history characteristics. As expected, PTSD patients had significantly higher PTSD symptom severity and more co-morbid major depressive disorder (MDD) than controls (Table 1).

Table 1

Demographic, trauma history and clinical characteristics of discovery cohort

	PTSD patients (n = 34)	Traumatized controls (n = 39)	P-value
Age	41.45 (8.96)	40.15 (10.12)	0.569a
Gender (females)	16 (47.1%)	19 (48.7%)	0.887b
Years of service	17.74 (13.43)	18.60 (9.93)	0.377c
Educational level			0.745d
Low	0 (0.0%)	0 (0.0%)
Middle	29 (87.9%)	33 (84.6%)
High	4 (12.1%)	6 (15.4%)
Current smoker	7 (21.9%)	8 (20.5%)	0.889b
PTSD symptom severity (CAPS total score)	70.21 (13.43)	4.59 (4.73)	<0.001a
Current co-morbid major depressive disorder	10 (29.4%)	0 (0.0%)	<0.001d
Work-related traumatic events (PLES total score)	19.13 (8.65)	19.74 (6.78)	0.740a
Childhood traumatic events (ETI total score)	5.64 (469)	3.95 (3.80)	0.116c

Continuous variables presented as mean (standard deviation); categorical variables presented as frequency

(percentage)

PTSD post-traumatic stress disorder, MDD major depressive disorder, CAPS clinician administered PTSD scale, PLES Police LifeEvent Scale, ETI Early Trauma Inventory

aIndependent t-test

bFisher’s exact test

cIndependent t-test on log-transformed data

dChi-square test

Sample preparation and quality control 450K array

Quantile normalization yielded proper data normalization (Figure S2). Table S2 lists the estimated relative blood cell count for which we did not observe significant group differences. Principal components (PC) one and two explained most of the variance in the 450K data set, 10.1 and 4.1%, respectively (Figure S3A). Pearson’s correlation analyses were performed between the first eight PCs with potential covariates/confounders. As several blood cell count estimates showed small to moderate correlations (r’s 0.2–0.5) with PCs 1 and 2, all blood cell estimates were included as covariates in further analyses. Additionally age and sex were included as covariates in further epigenetic association analyses, although these showed low correlations with PCs 1 and 2 (Figure S3B/C).

We observed no significant global DNA-methylation group difference. Furthermore, we found no genome-wide significant PTSD-associated DMPs (Table S3-tab DMPs, and Figure S4). Subsequently, we addressed a hypothesis-driven approach, evaluating previously reported PTSD associated genome-wide DMPs in our cohort. Our literature search yielded two genome-wide epigenetic surveys for PTSD, both using Illumina 27K arrays. Uddin et al. reported only on PTSD-associated regions or nearby genes, while Smith et al. reported on individual PTSD-associated DMPs [20, 21]. Therefore, only the latter was considered appropriate for DMP replication. None of the 21 significant loci reported by Smith et al. reached significance in our study (P-values < 2.4E−3 were assumed significant, Bonferroni-corrected). Nevertheless, for these 21 DMPs we did observe an overlap of 94.7% with regard to the direction of effect, i.e. hypo or hyper methylation in PTSD cases. Finally, we performed GSE analyses on the top 100 DMPs. DMPs located within the vicinity of 1 kb of a gene, annotated according to the Illumina manifest, were included. We observed 12 enriched pathway-based sets (49 of 74 genes entered GSE) of which six reached significance according to an adjusted P-value (q) < 0.05. The lowest P-value we detected involved the dopaminergic neurogenesis pathway (q = 0.012—wikipathways database), including three DMPs located at PITX3, NEUROD1, and the RET genes (Table S3—tabs top100 DMPs/GSE analyses and Figure S6). A secondary GSE analysis was based on the fifth category of the Gene Ontology (GO) database, which yielded only two significant enriched GO terms, both related to molecular functions involved in DNA binding (Table S3—tab top100 DMPs/GSE analyses).

There were no DMRs between PTSD patients and controls that reached genome-wide significance, however for one DMR we did observe consistent hyper methylation over more than 75% of the probes within the cluster (Table 2), located at the PAX8 gene (Fig. 1a). Technical validation of 450K data, i.e. including PAX8 (Fig. 1a and SOSTDC1 Fig. 1b) using MiSEQ technology was successful (r2 = 0.85, Fig. 1c) and described in detail in appendix 1.

Table 2

DMR association analysis PTSD vs. traumatized controls

Gene	chr: start-end	Value	Area	FWER	L	Cluster (L)	% L
PAX8	chr2:113992762–113993142	0.07	0.42	0.82	6	8	0.75
LDHC	chr11:18433500–18433564	0.07	0.20	0.98	3	9	0.33
FAM71F1	chr7:128354862–128354967	−0.07	0.20	0.98	3	7	0.43

Detected DMRs (L > 2) using minfi’s “bumphunter” function in whole sample, female and male stratified analyses. chr: chromosome and position (hg19); Value: average difference in methylation in the bump with respect to controls, <0 represents hypo methylation in PTSD and >0 represents hypermethylation in PTSD; Area: area of the bump with respect to the 0 line; FWER: multiple test adjusted P-value with respect to the area, according the familywise error rate; L: number of probes in DMR; cluster(L): number of probes in cluster; %L: percentage of probes within the cluster contributing to the DMR—in bold ≥ 75%

Fig. 1

Technical validation of 450K data. a DMR PAX8 using the bsr amplicon MiSEQ method and b 450K data SOSTDC1 using the bsr amplicon MiSEQ (lower plot) method. For both a and b X-axis plots represents chromosomal position (hg19), Y-axis plots represents β-values of 450K data and MiSEQ data, respectively. Pink line represents PTSD cases and blue line represents traumatized controls. Shading represent 95% C.I. c Overall correlation plot of 450K and MiSEQ detected CpG sites

For the genetic survey we analyzed 13.344 CpG-SNP probes (maf > 0.01), which yielded 15 CpG-SNPs potentially associated with PTSD. Since the statistical power in this study is limited due to the small sample, only the top 5 CpG-SNPs were investigated further (Figure S5). These five candidate CpG-SNPs were located in or nearby the OR2V2, EIF3B, SDK1, SKA2, and CACNA1C genes (Table 3, column discovery P-value). Next we aimed to validate these five candidate CpG-SNPs using sanger sequencing (Table S4). We did not detect any genotype discrepancies between 450K-based and sanger-sequencing-based genotypes. Three of five candidates, located at the genes SDK1, SKA2, and CACNA1C, were significantly associated with PTSD in this validation (Table 3, column validation P-value). Next we performed replication analysis for the five candidate CpG-SNPs in an independent cohort of trauma-exposed children with and without PTSD (see appendix 1 and Table S5). We significantly replicated the association of one CpG-SNP (rs1990322) located at the CACNA1C locus with PTSD (Table 3, column replication P-value. In addition, for this replicated CpG-SNP at the CACNA1C locus we searched for evidence for genetic interaction with flanking 450K probes (cis meQTLs) within a window of 500 bp upstream and 500 bp downstream in the discovery dataset. We observed no indications for meQTL effects of rs1990322 and its neighboring CpG sites (Figure S7).

Table 3

Top 5 CpG-SNPs genetic association analysis PTSD vs. traumatized controls

rs ID	HGVSa (chr)	Probe ID	Gene	Discovery (P-valueb)	Validation (P-valuec)	Replication (P-valuec)
rs2546424	g. 180581761C>T (5)	cg22931151	OR2V2	0.0009	>0.01	>0.01
rs7808948	g. 2425741C>T (7)	cg03122926	EIF3B d	0.0040	>0.01	>0.01
rs13230345	g. 4118584G>A (7)	cg22535849	SDK1	0.0027	0.005	>0.01
rs7208505	g. 57187729G>A (17)	cg13989295	SKA2	0.0020	0.009	>0.01
rs1990322	g. 2760970G>A (12)	cg24393317	CACNA1C	0.0006	0.006	0.001

chr chromosome

aHuman Genome Variation Society nomenclature (hg19), the probe position and position of the SNP might show a difference of one bp since some of our reported variants are exactly located on the cytosine of the CpG while other variants are located on the guanine of the CpG which results in exactly the same methylation pattern

bDetection P-value-based linear model association

cValidation and replication P-value-based on Chi-square tests, replication P < 0.01 was assumed significant (Bonferroni, five tests)

dThe rs7808948 is located in the intergenic region between the EIF3B and CHST12 genes

Discussion

In this study we aimed to detect epigenetic and genetic loci associated with PTSD in a homogeneous cohort of traumatized police officers. Both a genome-wide and hypothesis-driven replication approach did not result in DMPs between PTSD patients and trauma-exposed controls. GSE analysis on the top 100 DMPs showed, however, a plausible association of the dopaminergic neurogenesis pathway with PTSD. Furthermore, we observed one DMR located at the PAX8 gene suggesting consistent hypermethylation in PTSD patients. Genetic analyses yielded three CpG-SNPs significantly associated with PTSD. Of these, one CpG-SNP, located at the CACNA1C locus, was also significantly associated with PTSD in an independent replication sample of trauma-exposed children. Notably, this result shows that the Illumina 450K array is not restricted to epigenetic surveys but can provide informative genetic data as well.

Although our sample was small, it was highly homogenous as all participants were former or current police officers, and cases and controls were matched for sex, age, education, and years of police service. All participants reported multiple prior traumatic events, without significant group differences in reported types of traumatic experiences. PTSD patients fulfilled current diagnostic criteria for PTSD, while our trauma-exposed controls had minimal PTSD symptoms and did not report lifetime PTSD or other trauma-related psychiatric disorders. Thus our controls were apparently resilient to adverse mental health outcome of trauma. This study design, including extreme phenotypes following similar trauma load, was considered to favor detection of PTSD-associated loci, as also suggested by others [22]. Nevertheless, our genome-wide survey clearly remains limited in statistical power. Post hoc evaluation of the statistical power of our genome wide epigenetic survey was based on the report of Tsai and Bell [30]. If we would assume that our effect sizes ranged between 5 and 10% and we would aim for a 80% power, we would need approximately a sample size of 100 per group. The absence of genome-wide significant DMPs is therefore not unexpected. With regard to replication of literature reports on genome wide significant DMPs associated with PTSD, our search resulted in one previous study [20] reporting genome-wide detected DMPs, albeit in a different ethnic group and with substantial lower probe density of the array (27K) [20]. Such methodological differences may explain why we were unable to replicate these DMPs although we did observe a 94.7% consistent direction of effect between our study and the study of Smith et al.

GSE analysis using the top 100 DMPs revealed three genes, PITX3, NEUROD1, and RET, annotated in the pathway of dopaminergic neurogenesis pathway (wikipathways), which could very plausibly be involved in PTSD. Considering its role in stress reactivity and anhedonia, the dopamine system has previously been implicated in PTSD, and dopamine dysfunctions have been observed in PTSD patients [31]. Interestingly, these genes are involved in neuron differentiation and/or maturation but not neuron regionalization and specification. Differentiation and maturation processes might be more susceptible to adaptation under the influence of the environment, which is clearly relevant for PTSD. Further translational research should be performed in order to unravel the role of this pathway in PTSD.

Compared to DMP analyses, DMR analyses suffer less from statistical power issues, since these follow a multi-probe association methodology implying a reduced multiple test penalty [32]. We detected no DMRs with genome-wide significance after family-wise error rate correction. Such corrections however may be considered too conservative in genome wide screens, given that no false positives are allowed [33]. We did observe one suggestive associated DMR in the PAX8 gene that we consider a promising candidate for further studies on PTSD, since it involved >75% concordant differential hypermethylated probes within the cluster window in the PTSD patients. The PAX8 gene encodes a transcription factor involved in thyroid and central nervous system (CNS) development. Interestingly, thyroid hormone alterations have been observed in PTSD patients [34]. Additionally, PAX8 has been repeatedly implicated in sleep disturbances, which constitute an important symptom of PTSD [35]. Furthermore, Siegmund et al. reported hypomethylation of PAX8 in schizophrenia patients, suggesting a broader role for this gene in psychiatric disorders [36].

To our best knowledge, we are the first to use the ~14K CpG-SNPs with maf > 0.01, present on the 450K DNA-methylation array of Illumina for genetic association analyses. Although this CpG-SNP set is far too small for imputation and for a substantial number of these CpG-SNP genotypes cannot be reliably scored due to interfering DNA methylation effects, we were still able to detect five CpG-SNPs potentially associated with PTSD, of which three (SDK1, SKA2, and CACNA1C) were significantly associated with PTSD in the validation analyses. A statistical power estimation (www.osse.bii.a-star.edu.sg) based on the sample size of our discovery cohort showed that in order to achieve a power ≥ 80% at an alpha of 5%, we should observe a maf difference between cases and controls of ~30%. Since we have observed maximum maf differences of 22% (for the CACNA1C and SKA2 annotated CpG-SNPs) between cases and controls, this indicated a limited statistical power and thus the importance of replication.

Importantly, one of these CpG-SNPs (rs1990322, cg24393317, CACNA1C) was also significantly associated with PTSD in an independent cohort of traumatized children with and without PTSD. This locus is located in an intronic region of the CACNA1C gene. This gene encodes for an alpha-1C subunit of the L-type voltage-dependent calcium channel, crucial for generating neuronal action potentials. Several genome-wide and candidate-gene studies have reported associations between this gene and susceptibility for psychiatric disorders, including MDD [37], and bipolar disorder and schizophrenia [38]. Additionally, CACNA1C SNPs repeatedly have been found to affect brain structure and function, including amygdala responsivity, commonly considered a biological hallmark of PTSD [39]. However, in a candidate-gene study, another common SNP at CACNA1C locus did not moderate between hurricane exposure and subsequent PTSD in Hispanic Blacks [40].

Although our PTSD-associated SNP (rs1990322) is located in the CACNA1C gene itself, we cannot rule out a long-range interaction of this SNP with the upstream neighboring gene FKB4. FKBP4 is a functional antagonist of FKBP5, modulating glucocorticoid receptor function. Although psychiatric studies of FKBP4 are limited, FKBP5 SNPs, methylation and gene expression have been implicated in several stress-related psychiatric disorders including PTSD [18], and pre-trauma FKBP5 gene expression was found to predict PTSD symptom development upon military deployment [41]. Further studies should elucidate the actual affected gene contributing to the expression of PTSD.

Although the association between PTSD and the SKA2 gene (cg13989295, rs7208505) in the discovery cohort was only technically validated but not replicated, this observation supports the accumulating evidence for the role of this gene in PTSD and PTSD-related clinical and biological phenotypes. SKA2 modulates glucocorticoid receptor function [42], which was previously found to predict PTSD symptom development after military deployment [41]. Within that same cohort it was recently observed that decreasing SKA2 methylation was associated with PTSD symptom development during deployment, depending on the exact genotype carried [43]. In addition, in another military cohort, genotype carrier status and genotype-adjusted methylation status of this locus was also associated with PTSD [44]. In highly trauma-exposed civilians, similar associations emerged for this locus after controlling for childhood trauma [43]. SKA2 genotype status, genotype-unadjusted and adjusted methylation, and gene and protein expression have also been associated with suicide risk, either directly or in interaction with trauma-exposure or anxiety level [45, 46]. As several studies reported differential SKA2-methylation between (prospective) psychiatric patients and controls when adjusting for genotype, we additionally investigated this in the discovery cohort conform methods described by Guintivano et al. [43, 46]. We found no evidence for epigenetic effects additional to the CpG-SNP association (data not shown). This may be explained by the fact that epigenetic effects at the SKA2 locus can only occur in C allele carriers and our sample includes a limited number of homozygous SKA2 C allele carriers, limiting statistical power to detect such an epigenetic effect.

We also observed a technical valid, but not replicated association between PTSD and CpG-SNP rs13230345 (cg22535849). This SNP is located in an intronic region (c.3203–511) of the SDK1 gene [47], encoding for the sidekick cell adhesion molecule 1. In silico evaluation of this variant showed that it might introduce a novel donor and acceptor splice site (Alamut Visual version 2.7, Interactive Biosoftware, Rouen, France), which in turn might affect transcription, leading to the absence or a dysfunctional protein. Interestingly, the SDK1 gene has been implicated in schizophrenia [48] and autism [49], and altered SDK1 methylation was observed in post-mortem brain tissue of MDD patients [47]. Furthermore, Bagot et al. investigated transcriptional regulation of stress-susceptibility in rats and found that SDK1 was associated with social interaction and anxiety behaviors [50].

It should be noted that despite our careful considerations with regard to the phenotypic character, i.e. PTSD cases and resilient controls, defining both our discovery and replication cohort, differences related to the mean age between these two cohorts might be involved in the absence of replication of the genetic variants in SKA2 and SDK1. While our controls in the discovery cohort typically reported multiple prior traumatic experiences (including various types of potentially traumatic experiences during childhood), the majority of controls in the replication sample reported a single traumatic experience. As PTSD risk increases dose-dependently with increasing number of traumatic experiences, it could be that the non-replicated SNPs are involved in PTSD susceptibility in the context of increasing trauma load. Such an interaction between SNPs and increasing trauma load was previously observed by Kolassa et al. [51] Such interaction would also explain why previously reported findings on SKA2 and SDK1 in relation to PTSD or related psychological traits in adults, is in line with our observations in our (adult) discovery cohort.

In summary, while our statistical power to detect genome-wide significant DMPs was low, we detected an enriched gene set in the dopaminergic neurogenesis pathway that plausibly plays a role in PTSD. In addition, we observed a suggestive association of one DMR located at the PAX8 gene that seems a good candidate for future studies in PTSD. Furthermore, genetic association analyses based on 450k methylation data, and subsequent validation and replication analyses revealed that the CpG-SNP located at the CACNA1C locus was consistently associated with PTSD both in our discovery cohort and replication cohort. Despite the fact that we were unable to replicate the CpG-SNPs located at SDK1 and SKA2 loci, previously reported findings of others and our data do suggest a role of these genes in PTSD.

Additional research has to be performed to elucidate whether the identified DMRs and SNPs are consistently and causally associated with development or expression of PTSD, and importantly with which specific symptoms and biological mechanisms they are associated. For this purpose, transdiagnostic approaches among different types of psychiatric disorders seem especially fruitful.

If replicated, these findings may eventually contribute to the development of PTSD biomarker sets that can potentially support or improve an early clinical diagnosis of PTSD or aid in identifying individuals at risk for PTSD.

Web Resources

Consensus pathway database web tool, developed by the Max Planck institute, www.consensuspathdb.org, version 31, 2015 [28].

Accession Numbers

Raw and processed data are freely available on reasonable request but under the strict condition that meta-data remains confidential. Summary statistics on DMPs are given in Table S3.

Supplementary Figures and tables

Supplementary file 1

EJHG video

blurb text EJHG tube video