Discovery of 95 PTSD loci provides insight into genetic architecture and neurobiology of trauma and stress-related disorders

Data collection and GWAS

The PGC-PTSD²¹ Freeze 3 data collection includes 1,307,247 individuals from 88 studies (Supplementary Table 1). Data in this freeze were assembled from three primary sources (Fig. 1A): PTSD studies based on clinician administered or self-reported instruments (Freeze 2.5^11,12 plus subsequently collected studies), MVP release 3 GWASs utilizing the Posttraumatic Stress Disorder Checklist (PCL for DSM-IV),¹³ and 10 biobank studies with electronic health record (EHR)-derived PTSD status. We included 95 GWASs, including EA (N=1,222,882; effective sample size (N_eff)=641,533), AA (N=51,034; N_eff=42,804) and LAT (N=7,017; N_eff=6,530) participants (Supplementary Table 2).

Comparisons of PTSD data subsets

Population, screening, and case ascertainment differences between datasets led to the assumption that there would be substantial cross-dataset variation in PTSD genetic signal. We investigated this possibility using the software MiXeR, which incorporates an LD score regression (LDSC) genetic correlation model plus an extension to estimate sets of causal variants and model their shared and unshared polygenic overlap.^22,23 In univariate analyses, we identified comparable genetic architecture between PGC2.5 and EHR datasets in regards to h²_SNP, polygenicity, and discoverability. However, the uniformly assessed MVP had a higher h² SNP, was more discoverable and less polygenic than the other datasets (p < 0.05; Supplementary Table 3). Bivariate analyses identified high genetic correlations between the three subsets (r_g range = 0.79 - 0.87) (Extended Data Fig. 1, Supplementary Table 4). The MiXeR model did not provide a substantially better fit than the r_g model (all AIC < 0), indicating that there is no evidence for subset-specific genetic causal variation.

European ancestry PTSD GWAS

Given the similarities of the PTSD subsets, we performed a sample-size weighted fixed-effects meta-analysis of GWAS. For the EA meta-analysis (137,136 cases and 1,085,746 controls), the GC lambda was 1.55, the LDSC²⁴ intercept was 1.0524 (SE = 0.0097) (Supplementary Table 5), and the attenuation ratio was 0.0729 (SE=0.0134), indicating that 92.7% of the observed inflation in test-statistics was due to polygenic signal; thus artifacts produced only minimal inflation.

The EA meta-analysis identified 81 independent genome-wide significant (GWS) loci, including 5 GWS loci on the X chromosome (Extended data Fig. 2, Supplementary Figs. 1 and 2, Supplementary Table 6, regional association plots in Supplementary Data 1, forest plots in Supplementary Data 2, Supplementary Text). Relative to recent prior PTSD GWAS, 67 loci are novel^11-13(Supplementary Table 7). No region exhibited significant effect size heterogeneity (Supplementary Fig. 3).

We next sought to gain insights into whether loci harbor multiple independent variants. While FUMA²⁵ annotations reported independent lead SNPs within risk loci based on pair-wise LD (Supplementary Table 8), COJO²⁶ analysis of each locus conditional on the leading variants suggested that only one locus carried a conditionally independent GWS SNP (rs3132388 on chromosome 6, p=2.86 ×10^-9). This locus however, is in the MHC region, whose complicated linkage disequilibrium (LD) structure²⁷ may not be accurately captured by reference panels.

African and Native American ancestry PTSD GWAS meta-analyses

The AA meta-analysis included 51,034 predominantly admixed subjects (N=11,560 cases and 39,474 controls). There was minimal inflation of test statistics, with GC lambda = 1.031. No GWS loci were identified (Supplementary Fig. 4). The LAT meta-analysis was performed in 7,017 subjects (N=2,064 cases and 4,953 controls). There was minimal inflation of test statistics (GC Lambda=0.993) and no GWS loci were identified (Supplementary Fig. 5).

Multi-ancestry GWAS meta-analysis

A multi-ancestry fixed-effects meta-analysis of EA, AA, and LAT GWAS (150,793 cases, 1,130,197 controls) identified 85 GWS loci. Compared to the EA meta-analysis, 10 loci lost GWS, while 14 previously suggestive loci (p < 5 x 10^-7) became GWS (Fig. 2). In total, the present study identified 95 unique GWS PTSD loci between the EA and multi-ancestry meta-analyses (Table 1). Due to the complex local ancestry structure in individuals with African and Native American ancestry, which complicates LD modeling, we focused subsequent fine-mapping analyses (Fig. 1B) on data from the EA GWAS.

• Download figure
• Open in new tab

Figure 2: GWAS meta-analyses in European and multi-ancestry individuals identify a total of 95 PTSD risk loci.

Overlaid Manhattan plots of European ancestry (EA; 137,136 cases and 1,085,746 controls) and multi-ancestry meta-analyses (150,760 cases and 1,130,173 controls), showing 81 genome-wide significant (GWS) loci for the EA (full circles) and 85 GWS loci for the multi-ancestry (hollow circles) analyses. Circle colors alternate between chromosomes, with even chromosomes colored blue and odd chromosomes colored black. The y axis refers to −log₁₀ p-values. The horizontal red bar indicates the threshold for GWS associations (p < 5×10-8).

Gene-mapping

To link GWS SNPs to relevant protein coding genes, we applied three gene mapping approaches implemented in FUMA: positional mapping, expression quantitative trait loci (eQTL), and chromatin interaction mapping (Supplementary Table 9). GWS SNPs within the 81 EA loci mapped to 415 protein coding genes under at least one mapping strategy. A total of 230 genes (55%) were mapped by two or more strategies, and 85 (20%) genes were mapped by all three strategies (Supplementary Fig. 6). Notably, some genes were implicated across independent risk loci by chromatin interactions/eQTL mapping, including EFNA5, GRIA1, FOXP2, MDFIC, WSB2, VSIG10, PEBP1, and C17orf58. Chromatin interaction plots are shown in Supplementary Data 3.

Functional annotation and fine-mapping of risk loci

Functional annotations were used to gain insights into the functional role of SNPs within the 81 risk loci (Supplementary Table 10): 72 loci contained at least one SNP with Combined Annotation Dependent Depletion (CADD)²⁸ scores suggestive of deleteriousness to gene function (≥12.37), 43 loci contained GWS SNPs with Regulome DB²⁹ scores likely to affect binding, and 23 loci contained at least one SNP in the exon region of a gene.

To narrow the credible window of risk loci and identify potentially causal SNPs, we fine-mapped loci using Polyfun+SUSIE³⁰, which identified a credible set for 67 loci. Credible set window lengths were on average 62% of the original set lengths (Supplementary Table 11) and contained a median of 23 credible SNPs (range 1-252). Only one contained a SNP with posterior inclusion probability > 0.95, a missense SNP in the exon of ANAPC4 (rs34811474, R[CGA]>Q[CAA]; Supplementary Table 12).

Gene-based, gene-set, and gene-tissue analyses

As an alternative approach to SNP-based association analysis, we tested the joint association of markers within genes using a gene-based association analysis in MAGMA,³¹ which is a 2-stage method that first maps SNPs to genes and then tests whether a gene is significantly associated with PTSD. The gene-based analysis identified 175 GWS genes (Supplementary Table 13, Supplementary Fig. 7). Of these, 52 were distinct from the genes implicated by the gene-mapping of individual SNPs within GWS loci. These notably include DRD2, which has been thoroughly investigated in the context of psychiatric disorders and is a significant GWAS locus for multiple psychiatric disorders including schizophrenia.³² Refer to the Supplementary Text and Supplementary Table 14 for further investigation of conditionally independent SNPs within these 52 genes.

MAGMA gene-set analysis of 15,483 pathways and gene ontology (GO) terms from MSigDB³³ identified 12 significant GO terms. Significant terms were related to the development and differentiation of neurons (e.g. go_central_nervous_system_development, p=2.0×10^-7), the synaptic membrane (e.g. go_postsynaptic_membrane, p=6.9×10^-7), regulation (go_positive_regulation_of_gene_expression 1.0×10^-6), and nucleic acid binding (p=1.52×10^-6) (Extended Data Fig. 3, Supplementary Table 15).

MAGMA gene-tissue analysis of 54 tissue types showed PTSD gene enrichment in the brain (most notably in cerebellum, but also cortex, hypothalamus, hippocampus and amygdala) and in the pituitary, with enrichment found across all 13 examined brain regions (Extended Data Fig 4). Cell type analysis conducted in midbrain tissue data³⁴ identified GABAergic neurons, GABA neuroblasts, and mediolateral neuroblast type 5 cell types as having enriched associations above other brain cell types tested (p<0.05/268) (Extended Data Fig 5). GABAergic neurons remained significant (p=4.4×10^-5) after stepwise conditional analysis of other significant cell types.

Multi-omic investigation of PTSD

To gain insights into which particular genes in enriched brain tissues were contributing to PTSD, we conducted a combination of a transcriptome-wide association study (TWAS)³⁵ and summary based mendelian randomization (SMR) analyses³⁶ using GTEx brain tissue data based on the EA GWAS summary data. TWAS identified 25 genes within 9 loci with Bonferroni-significantly different genetically regulated expression levels between PTSD cases and controls (p<0.05/14,935 unique genes tested) (Fig. 3A, Supplementary Fig. 8, Supplementary Table 16). SMR identified 26 genes within 4 loci whose expression levels were putatively causally associated with PTSD (p<0.05/9,003 unique genes tested) (Fig. 3B, Supplementary Table 17). Many of these genes have been previously implicated in PTSD37 and other psychiatric disorders (e.g., CACNA1E, CRHR1, FOXP2, MAPT, WNT3). Notably, the 3p21.31 (incl., RBM6, RNF123, MST1R, GMPPB, INKA1), 6p22.1 (incl., ZCAN9 and HCG17) and 17q21.31 (incl., ARHGAP27, ARL17A, CRHR1, MAPT, FAM215B, LRRC37A2, PLEKHM1, and SPPL2C) regions contained >10 putative causal genes each.

• Download figure
• Open in new tab

Figure 3: Manhattan plots of PTSD associations in multi-omic analyses.

Gene expression data from 13 brain tissue types and the pituitary were used to conduct a, Transcriptome-wide association study (TWAS) identifying 9 loci with differential expression between PTSD cases and controls and b, expression quantitative trait locus summary based mendelian randomization (eQTL SMR) identifying 4 loci where gene expression has putative causal effects on PTSD. c, Blood protein quantitative trait locus (pQTL) SMR identify 16 blood proteins whose abundance has a putative causal effect on PTSD. The y axis refers to -log₁₀ p-values of each respective analysis. The horizontal red bars indicate gene-wide significance (p < 0.05/14,935 for TWAS, p < 0.05/9,903 for eQTL SMR, and p < 0.05/1,209 for pQTL SMR). Significant findings are labeled.

Among the GTEx tissues with the most TWAS and SMR signals was the dorsolateral prefrontal cortex (dlPFC). To gain insight into cell type resolution, we conducted MAGMA for cell-type-specific markers of dlPFC and cell-type-specific SMR. MAGMA showed a significant enrichment of dlPFC inhibitory and excitatory neurons, but also of oligodendrocytes and oligodendrocyte precursor cells (Supplementary Table 18), while the SMR analyses identified cell-type-specific SMR signals for 8 genes (KANSL1, ARL17B, LINC02210-CRHR1, LRRC37A2, ENSG00000262633, MAPT, ENSG00000273919, PLEKHM1) over 3 loci (6 out of 8 from 17q21.31) and all cell types (p<0.05/1,885 unique genes tested) whose expression levels were potentially causally associated with PTSD (Supplementary Table 19). The top-gene, KANSL1, was significant in all cell types.

Given previously reported associations between blood-based protein levels and PTSD,^38,39 we performed protein quantitative trait loci (pQTL) SMR³⁶ analysis for PTSD using data from the UK Biobank Pharma Proteomics Project⁴⁰ (N=54,306 samples and N=1,209 proteins). We identified 16 genes within 9 loci whose protein levels were significantly associated with PTSD (p<0.05/1,209 and p HEIDI > 0.05) (Fig. 3C, Supplementary Table 20), including members of the TNF superfamily (e.g., CD40, TNFRSF13C) implicating TNF-related immune activation in PTSD.

Gene prioritization

One research objective was to identify the genes with the greatest evidence of being responsible for the associations observed at each identified PTSD locus. Following recent research methods,41 we prioritized genes based on weighted sum of evidence scores taken across the functional annotation and post-GWAS analyses (Fig. 1B). Based on the absolute and relative scores of genes within risk loci, we ranked genes into Tier 1 (greater likelihood of being the causal risk gene) and Tier 2 (prioritized over other GWAS-implicated genes, but lower likelihood than Tier 1 of being the causal gene). 75% of loci contained prioritized genes (Tier 1 or Tier 2), the remaining loci did not contain any genes over the minimum threshold of evidence (score ≥ 4) to suggest prioritization. The prioritized genes for the top 20% of loci (ranked by locus p-value) are shown in Fig 4. A complete list of scores and rankings for all 415 protein coding genes mapped to risk loci is available in Supplementary Data 4.

• Download figure
• Open in new tab

Figure 4: Gene prioritization in PTSD loci.

Summary of evidence categories of prioritized genes (Tier 1 or 2) for the top 20% of PTSD loci (as ranked by leading SNP p-value). Locus number, prioritized genes within locus, gene locations (in terms of cytogenic band), and gene tier ranks (Tier 1, orange; Tier 2, blue) are indicated on the left. Categories of evidence are grouped and colored according to the domain they belong to. CADD scores, pLI scores and fine-mapping PIPs are written within their respective squares. The total weighted scores taken across all 9 evidence categories are shown on the rightmost squares. Abbreviations: eQTL, expression QTL; CI, chromatin interaction; CADD, combined annotation dependent depletion; RDB, regulome DB; pLI, predicted loss of impact; PIP, posterior importance probability; TWAS, transcriptome-wide association study; SMR, summary Mendelian randomization; pQTL, protein QTL;

We performed pathway enrichment analysis of the Tier 1 genes in SynGO. From Tier 1, 11 genes mapped to the set of SynGO annotated genes (CACNA1E, DCC, EFNA5, GRIA1, GRM8, LRFN5, MDGA2, NCAM1, OLFM1, PCLO, and SORCS3). Relative to other brain-expressed genes, Tier 1 genes were significantly overrepresented in the synapse (p=0.0009, qFDR=0.003), pre- and post-synapse (p=0.0086, qFDR=0.0086 and p=0.003, qFDR=0.004, respectively), and four subcategories (Extended Data Fig. 6). By contrast, there was no significant overrepresentation of genes when we applied this test to the entire set of 415 protein coding genes. Other notable Tier 1 genes included PDE4B related to synaptic function and TNF-related immune-regulatory genes, including TANK and TRAF3.

Genetic architecture of PTSD

SNP-based heritability (h²_SNP) estimated by LDSC was 0.053 (SE=0.002, p=6.8×10^-156). Whereas previous reports suggested sex-specific differences in PTSD,¹¹ no significant differences were found (p=0.13), and r_g between male and female subsets was high (r_g=0.98,SE=0.05, p=1.2×10^-98; Supplementary Table 5). MiXeR estimated 10,863 (SE=377) influential variants and a discoverability of 7.4×10^-6 (SE=2.2×10-7) (Supplementary Table 3), indicating a genetic architecture comparable to other psychiatric disorders.⁴²

Partitioned heritability across 28 functional categories identified enrichment in histone markers (H3K9ac peaks: 6.3 fold enrichment, SE = 1.12, p=3.11×10^-6; H3K4me1: 1.5 fold enrichment, SE=0.14, p=3.3×10^-4; Supplementary Table 21), and in evolutionary constrained regions across 29 Eutherians (18.37 fold enrichment, SE = 1.18, p=1.29×10-17). This is consistent with findings for multiple other psychiatric disorders, but has not been previously identified in PTSD.⁴²

Contextualization of PTSD among psychiatric disorders

We measured the genetic overlap between PTSD and other psychiatric disorders using the most recent available datasets.^32,43-52 We observed moderate to high positive r_g between PTSD and other psychiatric disorders (Extended Data Fig. 7A). To gain further insights into this overlap, we used MiXeR to quantify the genetic overlap in causal variation between PTSD and bipolar disorder (BPD), MDD, and schizophrenia (SCZ) (Extended Data Fig. 7B). The strong majority (79-99%) of the variation influencing PTSD risk also influenced these disorders (Extended Data Fig. 7B, Supplementary Tables 22 and 23). Similar to r_g, PTSD had the highest fraction of concordant effect directions with MDD (among the shared variation) (87% concordant, SE=2%), significantly higher than the directional concordance with BPD (67%, SE=1%) and SCZ (65%, SE=0.5%).

While our results indicate an overall strong r_g between PTSD and MDD (r_g=0.85, SE = 0.008, p< 2×10^-16), the correlation between PTSD and MDD varied significantly across PTSD subsets, with the most homogeneously assessed subset, MVP, showing the lowest correlation, and the biobank subset being most strongly associated (Supplementary Table 24). Further, to evaluate if specific genetic regions differ substantially from genome-wide estimates we used LAVA⁵³ to estimate the local h² SNP and r _g of PTSD and MDD across the genome, as partitioned into 2,495 approximately independent regions (Supplementary Table 25). Local h² SNP was significant (P<0.05/2,495) for both PTSD and MDD in 141 regions. Of these, local r_g was significant (p < 0.05/141) in 40 regions, all in the positive effect direction, where the mean local rg² was 0.57 (SD=0.24). In addition, we assessed the local r_g between PTSD and MDD specifically for the 76 autosomal GWS EA loci (Supplementary Table 26). While LAVA identified 20 significantly correlated loci (r_g<6.58×10^-4), there was also evidence for PTSD loci lacking evidence for correlation with MDD (Supplementary Figures 9 and 10 showcase 6 selected loci with low and high r_g).

Contextualization of PTSD across other phenotype domains

Considering all 1,114 traits with SNP-based heritability z>6 available from the Pan-UKB⁵⁴ analysis, we observed Bonferroni-significant r_g of PTSD with 73% of them (Supplementary Table 27). Examining the extremes of estimates observed, the top positive r_g was with sertraline prescription (r_g=0.88, p=3.25 x 10^-20), a medication frequently prescribed for PTSD and other internalizing disorders⁵⁵. Other leading associations included medication poisonings (e.g. “Poisoning by psychotropic agents” r_g=0.88, p=3.92×10^-20), which could support a link with accidental poisonings or self-harm behaviors.^56,57 Converging with epidemiologic studies, there were correlations with gastrointestinal symptoms⁵⁸ (e.g., “Nausea and vomiting” r_g=0.80, p=2.39×10-16), mental health comorbidities⁵⁹ (e.g., Probable Recurrent major depression (severe)” r_g=0.87, p=1.18×10^-18 ; “Recent restlessness” r_g=0.86, p=4.21×10^-54;), chronic pain⁶⁰ (multi-site chronic pain r_g=0.63, p=7.5×10^-301) and reduced longevity^61-63 (“Mother’s age at death” (r_g=-0.51, p=7.6×10^-27).

Drug target and class analysis

We extended MAGMA gene-set analysis to investigate 1530 gene sets comprising known drug targets (Supplementary Table 28). We identified one drug (stanozolol, an anabolic steroid) significantly enriched for targets associated with PTSD (p=1.62×10⁻⁵). However, stanozolol has only two target genes in our analyses (ESR1, JUN), and likely reflects the strong association of ESR1 with PTSD in gene-level analyses (p=8.94×10^-12).

We further examined whether high-ranking drug targets were enriched for 159 drug classes defined by Anatomical Therapeutic Chemical (ATC) codes. We identified two broad classes where drugs were significantly enriched for association in drug target analyses (Supplementary Table 29). These were opioid drugs (ATC code N02A, p=2.75×10^-4), and psycholeptics (ATC code N05, p=3.62×10^-5), particularly antipsychotics (ATC code N05A, p=3.55×10^-7). However, sensitivity analyses limited to drugs with 10 or more targets identified no significant drug target sets nor drug classes.

Polygenic predictive scoring

We evaluated the predictive accuracy of PRS based on PTSD Freeze 3 in a set of MVP holdout samples (Fig. 5). In EA holdouts, risk was significantly different across the range of PTSD PRS: For example, individuals in the highest quintile of PTSD PRS had 2.4 times the relative risk of PTSD (log relative risk SE=0.032; 95%CI = [2.25, 2.56]; p=1.16×10^-167) than individuals in the lowest quintile. PRS explained 6.6% of the phenotypic variation in PTSD (Nagelkerke’s R² transformed to the liability scale at 15% population and sample prevalence), representing a major improvement over PRS based on Freeze 2. In contrast, among AA holdout samples, PRS explained only 0.9% (liability scale) of the variation in PTSD, consistent with previous work suggesting that AA PRS based on EA data lag behind in prediction.⁶⁵

• Download figure
• Open in new tab

Figure 5: Polygenic risk score analysis for PTSD across different data sets and ancestries.

PGC-PTSD Freeze 2 and Freeze 3 European ancestry (EA) based genetic risk score (PRS) predictions into independent samples of different ancestries. The y axis represents PTSD odds ratios relative to the lowest quintile of PRS. For EA, predictions based on Freeze 3 training data (10,334 cases and 55,504 controls; blue circles) demonstrate a significant performance increase compared to predictions based on the previous Freeze 2 training GWAS (Nievergelt et al. 2019; yellow circles). Based on Freeze 3 EA training data, EA individuals in the highest quintile of PRS have 2.8 fold the odds of PTSD relative to individuals in the lowest quintile PRS (blue circles). Lower prediction accuracies are found for individuals of African (AA; 10,151 cases and 22,420 controls; red circles)) and Native American (Latinx; LAT; 5,346 cases and 10,821 controls; purple circles) ancestries, indicating poor PRS transferability across ancestries.

Participants and studies

PTSD assessment and DNA collection for GWAS analysis were performed by each study following their protocols. A description of the studies included and the phenotypic and genotyping methods for each study Supplementary Methods and Supplementary Table 1. We complied with relevant ethical regulations for human research. All subjects provided written informed consent and studies were approved by the relevant institutional review boards and the UCSD IRB (protocol #16097×).

EHR Studies

A total of 10 EHR-based cohorts (not including the MVP, which also contributed data) provided GWAS summary statistics. These cohorts consisted of four US-based sites (Vanderbilt University Medical Center’s BioVu, the Mass General Brigham Biobank, Mount Sinai’s BioMe, and Mayo Clinic’s MayoGC) and six non-US sites (iPSYCH from Denmark, FinnGen, HUNT Study from Norway, STR-STAGE from Sweden, UK Biobank, and Estonia Biobank). More details on procedures at each site are provided in the Supplementary Text. At each site, a broad definition of PTSD cases was defined based on patients having at least 1 PTSD or other stress disorder code (see Supplementary Text for the list of corresponding ICD-9 and 10 codes). All other patients without such a code were defined as controls. From a total of 817,181 participants across all cohorts, this case definition resulted in 78,687 cases based on the broad definition (9.6%).

Data assimilation

Subjects were genotyped on Illumina (N=84 studies) or Affymetrix genotyping arrays (N=5 studies) (Supplementary Table 1). Studies which provided direct access to pre-quality control genotype data (N=64 studies) were deposited on the LISA server for central processing and analysis by the PGC-PTSD analyst. Studies with data sharing restrictions (N=24 studies) were processed and analyzed following their own site-specific protocols (Supplementary Table 28), and shared GWAS summary statistics for inclusion in meta-analysis.

Genotype quality control and imputation

Genotype data was processed separately by study. For genotype data processed by the PGC-PTSD analyst, quality control was performed using a uniform set of criteria, as implemented in the RICOPILI⁹³ pipeline version 2019_Oct_15.001. Modifications were made to the pipeline to allow for ancestrally diverse data and are noted where applicable. Quality control: using SNPs with call rates >95%, samples were excluded with call rates <98%, deviation from expected inbreeding coefficient (fhet < −0.2 or >0.2), or a sex discrepancy between reported and estimated sex based on inbreeding coefficients calculated from SNPs on X chromosomes. SNPs were excluded for call rates <98%, a > 2% difference in missing genotypes between cases and controls, or being monomorphic. Hardy-Weinberg equilibrium was calculated within only in the largest homogenous ancestry group found in the data. SNPs with a Hardy-Weinberg equilibrium P-value < 1 × 10−6 in controls were excluded.

After quality control, datasets were lifted over to the GRCh37/hg19 human genome reference build. SNP name inconsistencies were corrected, and genotypes were aligned to the strand of the imputation reference panel. Markers with non-matching allele codes or with excessive MAF difference (> 0.15) with the selected corresponding population in the reference data were removed. The pipeline was modified so that only the largest homogenous ancestry group in the data was used for the calculation of MAF. For ambiguous markers, strand was matched by comparing allele frequencies: if a strand flip resulted in a lower MAF difference between the study and the reference data, the strand was flipped. Ambiguous markers with high MAF (> 0.4) were removed. The genome was broken into 132 approximately equally sized chunks. For each chunk, genotypes were phased using Eagle v2.3.5 and phased genotypes were imputed into the Haplotype Reference Consortium panel⁹⁴ using minimac3. Imputed datasets were deposited with the PGC DAC and are available for approved requests.

Studies with data sharing restrictions followed similar criteria for quality control, as detailed in Supplementary Table 28 and in the references in the supplemental material. Studies were imputed to either the 1000G phase 3, HRC, SISu panel, or a composite panel. GWAS summary data were lifted to the GRCh37 reference build where required. As differences in the imputation panels and genome reference build can result in SNP-level discrepancies between datasets, each set of summary data was examined for correspondence to the centrally imputed data. Multi-allelic SNPs and SNPs with non-matching allele codes were excluded. Stand ambiguous SNPs with high MAF difference (>20%) from the average frequency calculated the PGC-PTSD data were flagged and examined for strand correspondence.

Ancestry determination

For studies where the PGC analyst had genotype data access, ancestry was determined using a global reference panel11 using SNPweights⁹⁵. The ancestry pipeline was shared with external sites to be utilized where possible. Subjects were placed into three large groupings: European and European Americans (EA; subjects with ≥90% European ancestry), African and African-Americans (AA; subjects with ≥5% African ancestry, <90% European ancestry, <5% East Asian, Native American, Oceanian, and Central-South Asian ancestry; and subjects with ≥50% African ancestry, <5% Native American, Oceanian, and <1% Asian ancestry), and Latinos (LAT; subjects with ≥5% Native American ancestry, <90% European, <5% African, East Asian, Oceanian, and Central-South Asian ancestry). Native Americans (subjects with ≥60% Native American ancestry, <20% East Asian, <15% Central-South Asian, and <5% African and Oceanian ancestry) were grouped together with LAT. All other subjects were excluded from the current analyses. For the MVP cohort, ancestry was determined using standard principal components analysis approach where MVP samples were projected onto a PC space made from 1000 Genomes Phase 3 (KGP3) samples with known population origins (EUR, AFR, EAS, SAS, and AMR populations). EHR cohorts followed their own site-specific ancestry classification protocols.

GWAS

GWAS was performed with stratification by ancestry group and study. Strata were only analyzed if they had a minimum of 50 cases and 50 controls, or alternatively 200 subjects total. Where noted (Supplementary Table 2), small studies of similar composition were jointly genotyped so that they could be analyzed together as a single unit. For GWAS, the association between each SNP and PTSD was tested under an additive genetic model, using a regression model appropriate to the data structure. The statistical model, covariates, and analysis software used to analyze each study is detailed in Supplementary Table 30. In brief, studies of unrelated subjects with continuous (case/control) measures of PTSD were analyzed using PLINK 1.9,⁹⁶ using a linear (logistic) regression model which included 5 PCs as covariates. For studies that retained related subjects, analyses were performed using methods that account for relatedness. QIMR was analyzed using GEMMA⁹⁷ v0.96, including the first five PCs as covariates. RCOG was analyzed using the generalized disequilibrium test.⁹⁸ UKBB was analyzed using Bolt-LMM⁹⁹ including 6 PCs, and batch and center indicator variables as covariates. VETS was analyzed using BOLT-LMM including 5 PCs as covariates. EHR based studies that included related subjects were analyzed using saddle point approximation methods to account for case/control imbalances. AGDS and QIM2 were analyzed using SAIGE¹⁰⁰ including 4 PCs and study specific covariates. BIOV was analyzed using SAIGE including 10 PCs and age of record. ESBB, FING, HUNT, and SWED were analyzed using SAIGE including 5 PCs. UKB2 was analyzed using REGENIE¹⁰¹ including 6 PCs, assessment center, and genotyping batch covariates. GWAS was additionally performed stratified by sex. For the X chromosome analysis, sex was added as a covariate.

Meta-analysis

Sample-size weighted fixed-effects meta-analysis was performed with METAL.¹⁰² Within each dataset and ancestry group, summary statistics were filtered to MAF ≥1% and imputation information score ≥0.6. Meta-analyses were performed within the EA, AA, and LAT ancestry groups. A multi-ancestry meta-analysis was performed as the meta-analysis of the three meta-analyses. Genome-wide significance was declared at P < 5 × 10⁻⁸. Heterogeneity between datasets was tested with the Cochran test. Markers with summary statistics in less than 80% of the total effective sample size were removed from meta-analyses. LDSC²⁴ intercept was used to estimate inflation of test statistics related to artifacts rather than genetic signal. The proportion of inflation of test statistics due to the actual polygenic signal (rather than other causes such as population stratification) was estimated as 1-(LDSC intercept-1)/(mean observed Chi-square-1).

Regional Association Plots

Regional association plots were generated using LocusZoom¹⁰³ with 1.5MB windows around the index variant (unless the locus region was wider than 1.5MB, in which case it was the locus region plotted plus an additional buffer to include data up to the recombination region). The LD patterns plotted were based on the 1000 Genomes Phase 3 reference data,¹⁰⁴ where a sample ancestry appropriate subpopulation (EUR, AFR, or AMR) was used.

Conditional analysis of significant loci

To determine if there were independent significant SNPs within risk loci, GCTA Conditional and Joint Analysis²⁶ was performed. Stepwise selection was performed using the --cojo-slct option and default parameters, where UKBB European genotype data was used to model LD structure.

SNP heritability

h² SNPof PTSD was estimated using LDSC. LD scores calculated within KGP3 European populations (https://data.broadinstitute.org/alkesgroup/LDSCORE/) were used for the input. Analyses were limited to HapMap 3 SNPs, with the MHC region excluded (chr6: 26–34 million base pairs). SNP-based heritability was also calculated as partitioned across 28 functional annotation categories (https://data.broadinstitute.org/alkesgroup/LDSCORE/) using stratified LDSC.¹⁰⁵

Comparisons of Genetic Architecture

We used univariate MiXeR (version 1.3)^22,23 to contrast the genetic architecture of phenotypes. MiXeR estimates SNP-based heritability and two components that are proportional to heritability: the proportion of non-null SNPs (polygenicity), and the variance of effect sizes of non-null SNPs (discoverability). MiXeR was applied to GWAS summary statistics under the default settings with the supplied European ancestry LD reference panel. The results reported for the number of influential variants reflects the number of SNPs necessary to explain 90% of SNP-based heritability. Bivariate MiXeR was used to estimate phenotype-specific polygenicity and the shared polygenicity between phenotypes. Goodness of fit of the MiXeR model relative to simpler models of polygenic overlap was assessed using AIC values. Heritability, polygenicity and discoverability estimates were contrasted between datasets using the z-test.

Local genetic correlation analyses

Local h² SNP and r _g between PTSD and MDD⁵⁰ were estimated using LAVA.⁵³ KGP3 European data was used as the LD reference. Local h²_SNP and r_g were evaluated across the genome, as partitioned into 2,495 approximately equally sized LD blocks. Local rg was only evaluated for loci where local heritability was significant (P < 0.05/2,495) in both phenotypes. Significance of local rg was based on Bonferroni adjustment for the number of r_g evaluated.

Polygenic risk scores (PRS)

PRS were calculated in ancestry-stratified MVP holdout samples, based on the EA Freeze 3 PTSD GWAS. GWAS summary statistics were filtered to common (MAF >1%), well-imputed variants (INFO > 0.8). Indels and ambiguous SNPs were removed. PRS-CS¹⁰⁶ was used to infer posterior effect sizes of SNPs, using the KGP3 EUR based LD reference panel supplied with the program, with the global shrinkage parameter set to 0.01, 1,000 MCMC iterations with 500 burn-in iterations, and the Markov chain thinning factor set to 5. PRS were calculated using the --score option in PLINK 1.9, using the best-guess genotype data of target samples, where for each SNP the risk score was estimated as the posterior effect size multiplied by the number of copies of the risk allele. PRS was estimated as the sum of risk scores over all SNPs. PRS were used to predict PTSD status under logistic regression, adjusting for 5 PCs. The proportion of variance explained by PRS for each study was estimated as the difference in Nagelkerke’s R2 between a model containing PRS plus covariates and a model with only covariates.

Functional Mapping and Annotation

We used the SNP2GENE module in FUMA²⁵ v1.4.1 (https://fuma.ctglab.nl) to annotate and visualize GWAS results. The complete set of parameters used for FUMA analysis are shown in the Supplementary Text. Independent genomic risk loci were identified (r² < 0.6, calculated using ancestry-appropriate KGP3 reference genotypes). SNPs within risk loci were mapped to protein coding genes using positional mapping (10KB window), eQTL mapping (GTEx v8 brain tissue,¹⁰⁷ BRAINEAC,¹⁰⁸ and CommonMind¹⁰⁹ data sources), and chromatin interaction mapping (PsychENCODE¹¹⁰ and HiC^111,112 of brain tissue types) methods. Chromatin interactions and eQTLs were plotted in circos plots. SNPs were annotated to functional annotation databases including ANNOVAR,¹¹³ CADD,28 and RegulomeDB.²⁹

Novelty of risk loci

The start and stop positions of independent risk loci were assessed for positional overlap with existing PTSD loci^11-13. Loci were declared novel if their boundaries did not overlap with a variant reported significant in prior GWAS.

MAGMA gene-based and gene-set analyses

Gene-based association analyses were conducted using MAGMA³¹ v1.08. SNPs were positionally mapped (0KB window) to 19,106 protein-coding genes. The SNP-wide mean model was used to derive gene-level p-values, with an ancestry appropriate KGP3 reference panel was used to model LD. Significance was declared based on Bonferroni adjustment for the number of genes tested. Gene-based association statistics were used in MAGMA for gene-set and gene-property analyses. Gene-set analysis used the MsigDB³³ version 7.0 including 15,483 curated gene-sets and gene-ontology (GO) terms. Gene-property analysis of tissues and tissue subtypes was performed using GTEx v8 expression data, with adjustment for the average expression of all tissues in the dataset. To evaluate cell type specific enrichment, the FUMA cell type module was used, selecting 12 datasets related to the brain (full list in Supplementary Text). Finally, MAGMA was used to estimate the enrichment of dlPFC cell types in PTSD risk based on the DER21 marker gene list from PsychEncode Consortium Phase 1 resource release.¹¹⁰

GWAS Fine-mapping

Polygenic functionally informed fine-mapping (Polyfun)³⁰ software was used to annotate our results data with per-SNP heritabilities, as derived from a meta-analysis of 15 UK Biobank traits. PTSD risk loci were fine-mapped using SUSIE,¹¹⁴ with these per SNP heritabilities used as priors, pre-computed UKB based summary LD information used as the LD reference, and locus start and end positions as determined by FUMA. The SUSIE model assumed a maximum of two causal variants.

Expression quantitative trait loci (eQTL) and blood protein quantitative trait loci (pQTL) analyses

To test for a joint association between GWAS summary statistics SNPs and eQTL, the SMR method,³⁶ a Mendelian randomization approach, was used. SMR software (version 1.03) was run using the default settings. The European samples of the 1000G were used as a reference panel. Bonferroni multiple-testing correction was applied on SMR P-value (P_SMR). Moreover, a post-filtering step was applied by conducting heterogeneity in dependent instruments (HEIDI) test. The HEIDI test distinguishes the causality and pleiotropy models from the linkage model by considering the pattern of associations using all SNPs significantly associated with gene expression in the cis-eQTL region. The null hypothesis is that a single variant is associated with both trait and gene expression, while the alternative hypothesis is that trait and gene expression are associated with two distinct variants. Finally, gene-trait associations based on SMR-HEIDI were defined as the ones for which P_SMR met the Bonferroni significance threshold and had P_HEIDI>0.05. We conducted a combination of SMR and HEIDI based on GTEx project latest (version 8) multi-tissue cis-eQTL databases¹⁰⁷ from 13 brain regions and pituitary tissue that showed significant enrichment in MAGMA/FUMA analyses (see above). We also used cell-type-specific eQTLs in dlPFC for SMR analyses.¹¹⁵ Finally, we used a blood UK Biobank pQTLs database of 1,463 plasma proteins⁴⁰ relying on a very large population (54,306) for SMR/HEIDI analysis to evaluate biomarker potential.

Brain focused TWAS

JEPEGMIX2-P¹¹⁶ software with default settings was used to conduct TWAS on 13 brain regions and pituitary tissue that showed significant enrichment in MAGMA/FUMA analyses using our PEC-DLPFC GReX model. JEPEGMIX2-P was applied on GWAS summary statistics to estimate gene-trait associations. This method was preferable since it relied on a covariance matrix based on 33K samples compared to other TWAS methods which use less than 3k samples.¹¹⁷ To determine significance, a Bonferroni correction threshold for the unique number of genes tested was applied) P < 0.05/14,935). As a less conservative approach, we also applied FDR at a q value threshold of 0.05.

Gene prioritization

Genes within risk loci were prioritized following the general approach previously described.⁴¹ Genes were given prioritization scores based on the weighted sum of evidence across all evidence categories: FUMA positional, eQTL, and CI mapping, variant and gene annotation scores (CADD, predicted loss of impact [pLI], and RDB scores), positional overlap in fine-mapping, significance in gene-based analyses, brain tissue TWAS, eQTL SMR, and pQTL SMR. Weights for each evidence category are provided in Supplementary Table 31. Within a given locus, the evidence scores were compared across genes to identify the most likely causal gene. Genes with scores >=4 were ranked as either Tier 1 (greater likelihood of being the causal risk gene) or Tier 2 (lower likelihood of being the causal risk gene) and genes with scores < 4 were left unranked. The ranking algorithm is as follows: For a given locus, if there was a gene whose evidence score >= 4 and this gene’s score was > 20% higher than all other genes in the locus, it was ranked as a Tier 1 gene (greater likelihood of being the causal risk gene). Within a locus with a Tier 1 gene, other genes with scores between 20% and 50% lower than the Tier 1 gene were labeled as Tier 2. For loci without a Tier 1 gene, all genes with scores >= 4 that were within 50% of the leading gene were ranked as Tier 2.

SynGO

PTSD related genes were tested for overrepresentation among genes related to synaptic terms in the SynGO¹¹⁸ web interface (https://www.syngoportal.org/). Brain expressed genes were selected as the background list for the overrepresentation tests. SynGO terms with FDR q < 0.05 were considered as being overrepresented.

Drug Targeting Analyses

Following a previously described approach,¹¹⁹ we analyzed the enrichment of gene-level associations with PTSD in genes targeted by individual drugs. We then examined the enrichment of specific drug classes among these drug target associations. We obtained gene-level associations using MAGMA³¹ v1.08. Variant-level associations were converted to gene-level associations using the “multi=snp-wise” model, which aggregates Z scores derived from the lowest and the mean variant-level P value within the gene boundary. We set gene boundaries 35 kilobases upstream and 10 kilobases downstream of the transcribed regions from build 37 reference data (National Center for Biotechnology Information, available at https://ctg.cncr.nl/software/magma).

We performed drug target analysis using competitive gene-set tests implemented in MAGMA. Drug target sets were defined as the targets of each drug from: the Drug–Gene Interaction database DGIdb v.4.2.0,¹²⁰ the Psychoactive Drug Screening Database Ki DB,¹²¹ ChEMBL v27,¹²² the Target Central Resource Database v6.7.0,¹²³ and DSigDB v1.0,¹²⁴ all downloaded in October 2020. We additionally used the drug target sets to identify targets of drugs of interest from gene-based analyses.

We grouped drugs according to the Anatomical Therapeutic Chemical class of the drug.¹¹⁹ Results from the drug target analysis were ranked, and the enrichment of each class in the drug target analysis was assessed with enrichment curves. We calculated the area under the enrichment curve and compared the ranks of drugs within the class to those outside the class using the Wilcoxon Mann-Whitney test. Multiple testing was controlled using a Bonferroni-corrected significance threshold of P < 3.27 × 10⁻⁵ for drug target analysis and P < 4.42 × 10⁻⁴ for drug class analysis, accounting for 1530 drug sets and 113 drug classes tested.

We initially limited drug target analyses to drugs with two or more targets. However, results suggested this low limit may lead to false positive findings. As a sensitivity analysis, we further limited these analyses to drugs with 10 or more targets. Multiple testing was controlled using a Bonferroni-corrected significance threshold of P < 5.42 × 10⁻⁵ for drug target analysis and P < 7.94 × 10⁻⁴ for drug class analysis, accounting for 923 drug sets and 63 drug classes tested.

Genetic correlations and causal associations with other phenotypes

Using LDSC, we assessed the r_g of PTSD derived from the PGC meta-analysis conducted in EUR cohorts with traits available from the Pan-UKB analysis conducted in EUR samples. Details regarding the Pan-UKB analysis are available at https://pan.ukbb.broadinstitute.org/. Briefly, Pan-UKB genome-wide association statistics were generated using the SAIGE and including a kinship matrix as a random effect and covariates as fixed effects. The covariates included age, sex, age x sex, age², age² x sex, and the top-10 within-ancestry principal components. We limited our analysis to data derived from UKB participants of European descent (N=420,531) because of the limited sample size available in the other ancestry groups. Initially, we calculated SNP-based heritability of phenotypes available from Pan-UKB, retaining only those with SNP-based heritability z>6 (Supplemental Table 25) as recommended by the developers of LDSC.¹²⁵ To define traits genetically correlated with PTSD, we applied a Bonferroni correction accounting for the number of tests performed.