Abstract
Despite experiencing a significant trauma, only a subset of World Trade Center (WTC) rescue and recovery workers developed posttraumatic stress disorder (PTSD). Identification of biomarkers is critical to the development of targeted interventions for treating disaster responders and potentially preventing the development of PTSD in this population. Analysis of gene expression from these individuals can help in identifying biomarkers of PTSD.
We established a well-phenotyped sample of 371 WTC responders, recruited from a longitudinal WTC responder cohort, by obtaining blood, self-reported and clinical interview data. Using bulk RNA-sequencing from whole blood, we examined the association between gene expression and WTC-related PTSD symptom severity on (i) highest lifetime Clinician-Administered PTSD Scale (CAPS) score, (ii) past-month CAPS score, and (iii) PTSD symptom dimensions using a 5-factor model of re-experiencing, avoidance, emotional numbing, dysphoric arousal and anxious arousal symptoms. We corrected for sex, age, genotype-derived principal components and surrogate variables. Finally, we performed a meta-analysis with existing PTSD studies (total N=1,016), using case/control status as the predictor and correcting for these variables.
We identified 66 genes significantly associated with highest lifetime CAPS score (FDR-corrected p<0.05), and 31 genes associated with past-month CAPS. Our more granular analyses of PTSD symptom dimensions identified additional genes that did not reach statistical significance in our overall analysis. In particular, we identified 82 genes significantly associated with lifetime anxious arousal symptoms. Several genes significantly associated with multiple PTSD symptom dimensions and lifetime CAPS score (SERPINA1, RPS6KA1, and STAT3) have been previously associated with PTSD. Geneset enrichment of these findings has identified pathways significant in metabolism, immune signaling, other psychiatric disorders, neurological signaling, and cellular structure. Our meta-analysis revealed 10 genes that reached genome-wide significance, all of which were down-regulated in cases compared to controls (CIRBP, TMSB10, FCGRT, CLIC1, RPS6KB2, HNRNPUL1, ALDOA, NACA, ZNF429 and COPE). Additionally, cellular deconvolution highlighted an enrichment in CD4 T cells and eosinophils in responders with PTSD compared to controls.
The distinction in significant genes between lifetime CAPS score and the anxious arousal symptom dimension of PTSD highlights a potential biological difference in the mechanism underlying the heterogeneity of the PTSD phenotype. Future studies should be clear about methods used to analyze PTSD status, as phenotypes based on PTSD symptom dimensions may yield different gene sets than combined CAPS score analysis. Potential biomarkers implicated from our meta-analysis may help improve therapeutic target development for PTSD.
Introduction
Posttraumatic stress disorder (PTSD) is a complex psychiatric disorder that can develop after experiencing a traumatic event. The attacks on the World Trade Center (WTC) on September 11, 2001 and their aftermath had a substantial impact on the physical and mental health of WTC rescue, recovery and clean-up workers, but only a subset developed PTSD. These differing clinical outcomes after experiencing trauma imply a role for biological and genetic influence in PTSD. Our cohort provides an unprecedented opportunity into PTSD insights, because they have been deeply phenotyped for a shared, specific trauma.
Understanding the biological mechanisms underlying PTSD will require careful dissection and analysis of many constituent symptoms and risk factors. PTSD is uniquely heterogeneous among psychiatric disorders, with complex and detailed diagnostic criteria that allow for 636,120 different combinations of symptoms1, and 79,794 different symptom combinations. Additional heterogeneity in PTSD stems from the type and extent of trauma. It has been well established that PTSD is a heterogenous disorder and that trauma type plays a role in differential outcomes. In the field of WTC exposures, some work has already been done to elucidate gene expression and clinical outcomes.2–10 Further, work by our group and others has demonstrated differential genetic heritability of PTSD according to trauma-type11, 12.
Identification of biomarkers will be critical to the development of targeted interventions for treating disaster responders and potentially preventing the development of PTSD in this population. Gene expression analysis from WTC responders is uniquely useful to deduce the biological heterogeneity in PTSD, given their exposure to a similar and well-documented trauma. Data on WTC-related traumatic exposures of responders analyzed here, in combination with their heterogenic clinical outcomes, makes this a critical study to understand PTSD development and chronicity after shared traumatic events. Although candidate gene expression and methylation studies have explored genes involved in canonical stress signaling pathways in PTSD, regulated by the hypothalamus-pituitary-adrenal (HPA) axis, and immune and sympathetic nervous systems, few have been able to control for length of time since exposure, nor so specifically delineate trauma type and secondary exposures such as dust cloud severity. While the WTC-related exposures experienced by rescue, recovery and clean-up workers in this cohort ranged in severity, the traumatic event –encompassing the 9/11 attacks and their aftermath– happened in a discrete, specific time window. Further, this sample is highly phenotyped with in-person clinical psychiatric evaluations, also including medical examination and laboratory testing.
The existence of this cohort and the generous participation of many responders to the WTC disaster enabled us to generate a large gene expression data set of 355 donors, to our knowledge the largest single traumatic event expression data set to date. We used the Clinician-Administered PTSD Scale (CAPS)13 score as a quantitative measure of PTSD symptom severity rather than a case/control definition, thus substantially increasing statistical power in this study.14–16 To our knowledge, ours is the first gene expression study to incorporate total CAPS scores and PTSD symptom dimensions as outcomes.
Methods
Participants
The WTC Health Program (WTC-HP) is a regional consortium of five clinical centers established in the greater New York City area by the Centers of Disease Control and Prevention in 2002, with the goal of providing health monitoring and treatment to WTC responders, comprising the WTC-HP General Responder Cohort17. We recruited participants from the WTC-HP Responder Cohort who had completed at least three periodic health monitoring visits at one of the four WTC-HP clinical centers participating in this study – Mount Sinai Medical Center, New York University, Northwell Health, and Rutgers/The State University of New Jersey – and who had provided signed consent to be contacted for research studies. Stratified random sampling was employed to ensure selection of WTC responders spanning the full range of WTC-related PTSD symptom severity, from no/minimal symptoms to severe/chronic PTSD symptom levels on the PTSD Checklist – Specific Version (PCL-S)18 completed during periodic health monitoring visits to the WTC-HP. Individuals with a lifetime history of chronic psychotic disorder or bipolar disorder type I, substance abuse/dependence or alcohol dependence over the prior three months, current pregnancy, acute medical illness or exacerbation of chronic medical illness, history of significant head injury or cerebrovascular accident, changes in medications or medication dosages over the prior month, or who were taking oral or regularly injected steroid medications were excluded from the study.
The study, conducted between April 2013 and September 2017, was approved by the Institutional Review Board of the Icahn School of Medicine at Mount Sinai, and all participants provided written informed consent. A total of 471 WTC responders completed in-person clinical assessments, yielding a final sample of 371 participants who met study eligibility criteria and completed study procedures, and 355 of those who had viable RNA-sequencing data (Figure 1). The mean age of participants was 54.1 (SD=8.3) years, 82% were male; ethnicity proportions are given in Table 1. The sample was composed of 40.8% police responders and 59.2% non-traditional responders (e.g., construction workers).
Assessments
Data on 10 WTC-related exposures 19(e.g., exposed to human remains, received treatment for an illness/injury during WTC recovery work) was obtained from interviews and self-report questionnaires completed by participants during their first health-monitoring visit to the WTC-HP, an average of 4.3 (SD=2.7) years following 9/11/2001. In-person clinical assessments were conducted an average of 13 (SD=2.3) years following 9/11/2011. Trained Masters- or PhD-level clinical interviewers administered the Structured Clinical Interview for DSM-IV (SCID)1 to assess current and lifetime Axis-I psychiatric diagnoses, and CAPS13, lifetime and past-month versions, to assess lifetime and past-month WTC-related PTSD diagnostic status and WTC-related PTSD symptom severity. Lifetime and past-month PTSD diagnosis was defined as meeting DSM-IV criteria for WTC-related PTSD and a total score ≥ 40 on the lifetime and past-month CAPS, respectively.
On the same day as the clinical assessment, participants also completed the Childhood Trauma Questionnaire (CTQ)20, assessing physical, sexual, and emotional abuse, and physical and emotional neglect experienced in childhood; the Traumatic Life Events Questionnaire (TLEQ)21, assessing lifetime exposure to a range of traumatic events (e.g, crime, natural disaster, assault); a checklist of 15 stressful life events they might have experienced since 9/11/2011 (e.g., “lost a job/laid off/lost income”, “divorced from spouse”, “had debt trouble”), modified from the Diagnostic Interview Schedule (DIS) Disaster Supplement22; and a health questionnaire asking which medical conditions they had ever been diagnosed with23, modified to add common WTC-related conditions (asthma or chronic respiratory condition, chronic rhinitis or sinusitis, sleep apnea, or acid reflux). Participants additionally completed a history and physical examination conducted by a licensed nurse practitioner, as well as routine laboratory testing, to rule out acute medical illness or exacerbation of chronic medical illness.
Among the 355 participants, 108 were determined to have met DSM-IV criteria for lifetime WTC-related PTSD, with 53 of them still meeting past-month PTSD criteria. The heterogeneity of PTSD confers some problems when attempting to analyze the disorder by case/control status alone. Case/control analyses do not fully capture the symptom complexity of the disorder, resulting in poor genomic modeling. Similarly, while overall PTSD symptom severity is a better quantitative measurement, it does not fully capture variability across PTSD symptomatology on a useful clinical level.24 To address this variability, we examined five symptom dimensions (re-experiencing, avoidance, emotional numbing, dysphoric arousal and anxious arousal symptoms), assessed with the CAPS to more accurately examine the heterogeneity of PTSD symptomatology (Figure 2).25
Blood sample collection and RNA extraction
Participants were instructed to fast after midnight and underwent collection of blood samples between 8:00 and 10:00 am. Total RNA was purified from whole blood collected in PAXgene blood RNA tubes (Qiagen, Germantown, MD, USA) using a PAXgene blood RNA kit IVD (Qiagen, Germantown, MD, USA). Total RNA concentration and quality were estimated using a NanoDrop 200c spectrophotometer according to the manufacturer instructions (Thermo Fisher Scientific, Waltham, MA, USA). Samples with an optical density ratio 260/280 superior or equal to 1.8 passed the quality control. Total RNA concentration and quality were also estimated using an Agilent RNA 6000 nano kit and an Agilent 2100 bioanalyzer according to the manufacturer instructions (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA Integrity Number superior or equal to 7 passed the quality control. RNA samples (derived from blood) were processed for RNA Seq with polyA selection and sequenced on Illumina HiSeq High Output mode with a sequencing configuration of 2×150 paired-end reads (GENEWIZ, South Plainfield, NJ). A total of 10M paired reads per sample was set as a threshold to account for high globin reads; 29 samples were re-sequenced to meet threshold.
Gene expression quality control analysis
We processed whole-blood gene expression data using the RAPiD.19 RNA-sequencing pipeline, and calculated normalized TPM counts from RSEM.26 We performed quality control analysis on the counts to verify sequencing and residual contributions to variance using VariancePartition27. We corrected each sequencing batch for sex, age, and genotype-derived principal components using Limma/voom weighted least-squares linear regression.28 We rank normalized and combined the residuals from the linear regression of each batch, and these values were used in all subsequent association tests for CAPS total score and five PTSD dimension scores (re-experiencing, avoidance, numbing, dysphoric arousal and anxious arousal).
Differential gene expression analysis
We used whole blood RNA-sequencing to test for associations between gene expression and WTC-related PTSD symptom severity on (i) highest lifetime CAPS score (CAPSL), (ii) past-month CAPS score (CAPSPM), and (iii) PTSD dimension scores including re-experiencing, avoidance, numbing, dysphoric arousal and anxious arousal, correcting for batch and surrogate variables (Eqn. 1). In addition, study participants had a wide range of psychiatric and somatic comorbidities, including some with substantial shared genetic aetiology and overlapping symptom profiles (e.g., major depressive disorder); comorbidities that may represent systemic manifestations of of PTSD (e.g., cardiovascular disease29); and exposure to the dust cloud during and following 9/11. We expect all of these factors to have substantial impacts on gene expression. Significant co-linearity between some of these measures and CAPS scores preclude including these variables as covariates within our analysis, and testing directly for their effect on gene expression (in particular, due to high correlations between length of time at the WTC site, CAPS score, and dust-cloud exposure; and between CAPS score and co-morbid medical disorders potentially constituting systemic manifestations of PTSD). Instead, in order to test whether these comorbidities and exposures might account for some of the CAPS-expression associations we observe, we also tested for gene expression associations with (i) an index of dust cloud exposure30; and (ii) number of medical comorbidities. Next, we tested for (i) interaction effects between each of these measures and CAPS score; (ii) enrichment of genome-wide significant associations between these measures and CAPS score; and (iii) genome-wide correlations in association statistics. For all gene expression analyses, we established significance using a Benjamini-Hochberg31 FDR correction < 5%.
Gene-set enrichment of PTSD-associated genes
We tested for gene set enrichment among our genes associated with CAPSL, CAPSPM, and PTSD dimension scores by (i) analyzing the significant genes from the association tests for pathway enrichment by gene permutation testing and (ii) analyzing all genes from the ranked association test gene lists to subject permutation using the R versions of GSEA32 and fgsea33, 34. For gene permutation testing, we included all nominally significant genes (p<0.05), and tested for association with 105 gene sets using Kyoto Encyclopedia of Genes and Genomes (KEGG) database35 for pathway enrichment.
We applied phenotype permutation testing rather than gene set permutation to keep the correlations between the genes in the dataset and the genes in the gene set pathways. For the subject permutation testing, each test was run with 10000 permutations, and pathways that passed Benjamini-Hochberg multiple testing correction were considered significantly enriched. To synthesize this large amount of gene set information, we generated comparative PTSD symptom plots using Clusterprofiler36 in R. Comparative gene set plots contained pathways which passed FDR<0.05 significance threshold.
Meta-analysis with existing gene expression analyses
To replicate our gene expression results, we meta-analyzed our data with association statistics from five other genome-wide gene expression analyses:
WTC responder study Stony Brook University (SB WTC) 3; N=282. Data are divided into discovery (WTC-d) and replication (WTC-r) cohorts.
Trauma Mega-Analysis study (TMA) 12. TMA combines 7 different PTSD studies and transformed the expression data into three categories: combat (N=169), male interpersonal (N=112), and female interpersonal trauma (N=259)12. We analyzed these data as separate trauma studies for the purposes of our meta-analysis (Table 2).
Since the majority of studies focus on PTSD case/control status, rather than associations with continuous CAPS scores (as here), we repeated our analysis to compare gene expression between PTSD cases (defined as meeting DSM-IV criteria for PTSD and a total CAPS ≥ 40) and controls (all others in our sample).
We meta-analyzed PTSD case-control association statistics using a sample-size based meta-analysis approach in METAL37. We included all genes from our analysis that reached p<0.05 (N=9,380 for past-month and N=1,016 for lifetime) and all available genes from other studies (N=27-806; Table 2).
Cellular deconvolution associated with CAPS scores
We applied CIBERSORT to our raw counts matrix to deconvolute immune cell types in our patients using the immune cell matrix reference panel LM2238. We tested for association between cell type proportions and CAPSPM, CAPSL.
Results
Gene expression is associated with PTSD symptom levels in WTC first responders
We tested for association between expression of 12,220 genes and total CAPS score in a sample of 355 WTC first responders. We identified 31 genes significantly associated with total past-month CAPS score (CAPSPM; Figure 3), and 66 genes associated with lifetime (highest) CAPS score (CAPSL, Figure 3). Of these, 42/66 genes are associated only with CAPSL, (and not CAPSPM), while 7/31 genes were associated only with CAPSPM, (and not CAPSL). Genome-wide associations with CAPSPM and CAPSL were significantly correlated (Beta ρ = 0.82, p<2.2×10-16; FDR-adjusted P-values ρ = 0.79, p<2.2×10-16) (Table 3).
We tested for enrichment of 59 well-studied PTSD candidate genes11 (Sup. Table 1) within our association statistics. We did not observe enrichment of these genes within our past-month or lifetime association analyses (p=0.174, 0.245). However, of these candidate genes, SERPINA1 was significantly associated with CAPSL score.
Environmental exposure to the WTC dust cloud
Next, we tested whether our genes associated with CAPS scores are specific to PTSD, or are driven by spurious associations with comorbid diagnoses or environmental exposure to the dust cloud at Ground Zero. In particular, a number of individuals within our study have comorbid conditions and complex medical histories (Table 4), including disorders with substantial shared genetic and environmental etiology with PTSD, and disorders and traits that may present as systemic manifestations of PTSD.
First, we tested whether these medical comorbidities alone may account for the associations we observe; we identified 175 gene associations with an aggregate summary score of medical comorbidities. Notably, these genes do not include any of our significant associations with CAPS scores. We identified only one gene, STX10, with significant interaction between CAPS scores and comorbid conditions; that is, gene expression was elevated specifically among individuals with both high CAPSL and a large number of comorbid diagnoses (Figure 4).
Next, we tested for association between gene expression and exposure to the dust cloud at Ground Zero30. We identified 561 genes significantly associated with this exposure. We tested for, but did not find, any significant interactions between CAPS scores and dust cloud exposure (p>0.05), and did not observe a large correlation of expression results in genome-wide significant genes between the two analyses for CAPSPM or CAPSL (R2=0.1835, R2=0.2466, p<0.05) (Figure 4). Together, these analyses imply that our gene-CAPS score associations are specific and relevant to PTSD, rather than due to confounding by comorbid diagnoses or dust cloud exposure.
Gene expression analysis reveals PTSD dimension-specific associations
Next, we tested for gene expression associations with past-month and lifetime PTSD symptom dimensions (re-experiencing, avoidance, dysphoric arousal, anxious arousal, numbing; Sup. Table 2).
Our analysis revealed overlapping and unique genes for each symptom dimension, and significant correlation of genome-wide association statistics between symptom dimensions. In particular, both our past-month and lifetime analyses identified a large number of genes (61 and 82, respectively; Sup. Table 2) associated with anxious arousal, including many genes not associated with any other phenotype in our analysis (16, 27, respectively). By contrast, although we identified a substantial number of genes significantly associated with dysphoric arousal (20, 16 for past-month and lifetime respectively), only two genes were uniquely associated with this phenotype (1 gene in past-month; 2 in lifetime). Only one gene, COPE, was significantly associated with every phenotype tested in our analysis. A second gene, EIF4A1, was also significantly associated with every lifetime (highest) phenotype (Figure 5).
PTSD pathway enrichment demonstrates immune, psychiatric, and metabolic relationships
Our genetic enrichment and pathway analyses identified well-established PTSD mechanisms and pathways, including pathways associated with inflammation, neurological signaling pathways, structural remodeling within and between cells, and HPA-axis and signaling39–44. Based on our KEGG pathway enrichment analysis, we revealed a set of genes significantly associated with psychiatric disorders pathways that were significantly enriched in our results: FURIN, PPP2R1A, GNAI2, PCLB2, and GNB2. The two symptoms that we discovered were associated with FURIN and their FDR-corrected p-values were avoidance (p=0.028, B=2.55) and numbing (p=0.034, B=-0.009). PPP2R1A was associated with anxious arousal (p=0.009, B=2.60). PCLB2 was only associated with numbing (p=0.046, B=-0.952). GNAI2 was associated with anxious arousal (p=0.042, B=-0.780) and numbing (p=0.047, B=-1.03), and GNB2 was only associated with avoidance (p=0.036, B=1.130) (Figure 6).
Immunological and metabolic gene enrichment was consistent across CAPSL score and numbing, but was less pronounced in anxious arousal. For past-month analyses, immune function was most associated with anxious arousal, whereas metabolic function was enriched in CAPSPM, and numbing, to a lesser extent. For both lifetime and past-month scores, neurological signaling pathways were most significantly pronounced in numbing, and were less prevalent in the overall total CAPS score analysis. For lifetime scores, structural pathway enrichment was significantly higher in total CAPS score, anxious arousal, and numbing, whereas for past-month scores, structural enrichment was most associated with anxious arousal (Figure 6).
Meta-analysis prioritizes 10 genes associated with PTSD
We sought to replicate our associations with previous PTSD studies. Since the majority of publicly available PTSD gene expression analyses follow a case-control, rather than quantitative measure (CAPS score) analysis, we converted our continuous CAPSPM and CAPSL values to PTSD case/control using (CAPS≥ 40) and DSM-IV PTSD-criteria, and repeated our analysis. Our case/control and CAPS score association statistics were significantly correlated (ρ=0.72, p<2.2×10-16, =0.79, p<2.2×10-16); however, we note a substantial decrease in the number of significantly associated genes when using a case-control design, compared to our initial quantitative analysis, as we would expect16. Twelve genes were significant for past-month PTSD (PTSDPM) case/control and 22 genes were significant for lifetime PTSD (PTSDL) case/control, versus 31 genes for CAPSPM and 66 genes for CAPSL.
We meta-analyzed our results with five publicly available cohorts (N=739 cases, 438 controls); two including WTC responders, and three including combat and interpersonal trauma (Table 2). For PTSDPM we identified 5 significant genes –COPE, CIRBP, FCGRT, NACA, and ZNF429 (p<5.33×10-6)–, and for PTSDL 8 significant genes –COPE, CIRBP, TMSB10, FCGRT, CLIC1, RPS6KB2, HNRNPUL1 and ALDOA (p<4.92×10-5)–, including genes associated with inflammation and immune response (Figure 7). Of these 10 genes, only one (COPE) had significant heterogeneity of effect size between cohorts: our study and TMA-combat. Three further genes were unique to our study; NACA p=3.34×10-6, CLIC1, p=1.9×10-5, and HNRNPUL1 p=4.08 x10-5 (Figure 7). The remaining six genes were significant across multiple studies in our meta-analysis, with highly consistent (negative) direction of effect (i.e., consistently decreased expression in cases compared to controls): ZNF439 (p=4.78×10-6), CIRBP (p=1.29 x10-6), TMSB10 (p=6.31 x10-6), FCGRT (p=1.12 x10-5), RPS6KB2 (p=3.47 x10-5), and ALDOA (p=4.66 x10-5) (Table 5).
Cellular deconvolution identifies differences in cell populations between responders with PTSD and controls
Since many of our PTSD-associated genes are related to immune function, we tested whether immune cell type proportions were correlated with CAPS scores in individuals in our sample. We performed cell-type deconvolution to identify cell-type proportions for 22 cell types across all 355 individuals in our sample. We found significant increase of CD4 naïve T cell (p<0.0049) proportions with CAPSPM, and significant increase of eosinophils (p< 0.042) and CD4 memory resting T cells (p<0.044) associated with CAPSL. In addition, we found significant decrease of activated natural killer cells (p<0.040) associated with CAPSL. (Figure 8).
Discussion
Although trauma is ubiquitous as a human experience, the types of traumatic experiences vary greatly among individuals. Our study sample and design present a unique opportunity to examine gene expression and PTSD symptoms related to a particular traumatic experience, the 9/11 WTC disaster and its aftermath. To our knowledge, ours is the largest single-traumatic event gene expression dataset to date. Using CAPS score as a quantitative measure of WTC-related PTSD symptom severity, we found 66 genes associated with lifetime (highest) CAPS score and 31 genes associated with past-month CAPS score. We additionally found interaction between CAPSL and one gene, STX10, in medical comorbidities. Pathway analysis links our associated genes to metabolic, immunological, structural, and neurological pathways. In addition, we found 10 genes associated with lifetime and past-month PTSD: COPE, CIRBP, and FCGRT shared between the two; NACA and ZNF429 specific to past-month PTSD; and TMSB10, CLIC1, RPS6KB2, HNRNPUL1 and ALDOA specific to lifetime PTSD. Furthermore, cellular analysis of our cohort demonstrated an enrichment in CD4 T cells and eosinophils in responders with PTSD. Additionally, natural killer cells were decreased in these patients. Our findings support previous studies which tie together the immune system and PTSD as a systemic disease. Ultimately, these results represent an additional wealth of knowledge to help understand the genetic expression and biological etiology of PTSD and uncover potential biomarkers for the disease.
Examining both lifetime and past-month CAPS scores allowed us to look at not only chronic effects of PTSD but longevity of the outcomes. While CAPSPM is considered the standard for case/control analysis, looking at both CAPSL and CAPSPM allows us to ask more specific questions about the role and relevance of elevated gene expression in PTSD. For example, genes associated with CAPSPM might represent current expression changes in PTSD, while those associated with CAPSL (lifetime CAPS score representing, for each responder, the highest WTC-related PTSD symptom levels ever reached since 9/11/2011) may represent long-lasting expression changes resulting from lifetime PTSD.
We identified genes associated with past-month and lifetime CAPS scores that have been associated with other psychiatric disorders, such as major depressive disorder, schizophrenia and autism: FURIN, PPP2R1A, PLCB2, GNAI2, and GNB2 (Figure 6)45–57. We additionally identified a group of genes as significantly associated with several lifetime and past-month PTSD symptom dimensions (SERPINA1, RPS6KA1, and STAT3) that have been previously linked to PTSD pathophysiology1, 27–31, and genes associated specifically with anxious arousal that have been previously associated with anxiety disorders: DCTN1, and FLI158, 59. One gene in our study, COPE, was associated with every PTSD phenotype in this analysis. Until now, COPE has not been well studied in the context of psychiatric disorders60, but it has been implicated in the context of Alzheimer’s disease61, 62. Canonically, the COPE protein is the epsilon subunit of the coatomer protein complex I (COPI) which regulates endocytosis from the plasma membrane and is involved in Golgi-to-lysosome transportation. Other subunits of COPI have been implicated in hereditary diseases that cause microcephaly and in autoimmune disorders63–65.
Importantly, our analytical design allowed us to test for potential confounding effects of comorbidities and environmental exposures66, 67. We expected many of these comorbidities and exposures to have a significant impact on gene expression, particularly as some comorbidities may be more recently occurring than, for example, the highest CAPSL measure. Therefore, we tested directly for genes associated with (i) comorbidities and (ii) dust cloud exposure. Although we identified a large number of genes associated with each phenotype, we note that our associations do not overlap with our top PTSD genes; we did not observe significant enrichment of shared associations, and only observed one gene with significant interaction effect between comorbidities and CAPSL. Therefore, we conclude that our results are not confounded by these exposures; our gene associations are specific to PTSD rather than broadly corresponding to exposure or general ill-health.
KEGG enrichment of our associated genes determined genetic changes in metabolic, immunological, structural, and neurological pathways associated with total CAPS, numbing, and anxious arousal phenotypes (Figure 7). Anxious arousal and numbing tend to have few significantly associated genes in common, while both past-month and lifetime CAPS scores display a few unique genes but many shared ones. GNAI2, while associated with many numbing-related pathways (long-term depression, oxytocin signaling, etc.), is also significantly associated with some peripheral anxious arousal-related and CAPSL–related pathways (leukocyte migration, and chemokine signaling, respectively). On the other hand, anxious arousal is uniquely associated with genes such as HGS, ARFGAP2 and RAB7A, which are linked to endocytosis and immunity. Similarly, our past-month and lifetime CAPS phenotypes are uniquely associated with glycolysis/gluconeogenesis for genes ALDOA, GAPDH, ENO1, TPI1, and PKM, which may play roles in HPA-axis or metabolic dysregulation (Table 4).
Our analyses identified 10 genes reaching genome-wide significance. Of these, 3/10 (NACA p=3.34×10-6, CLIC1, p=1.9×10-5, and HNRNPUL1 p=4.08 x10-5) are specific to our study; all three are downregulated in WTC responders with PTSD compared to controls. NACA has been previously studied as a potential biomarker for depression in mice under stress conditions68. The function of HNRNPUL1 is relatively unknown. Studies have suggested it might play a role in DNA damage repair and nucleocytoplasmic RNA transport. CLIC1 has been implicated in other psychiatric disorders, but its primary function is inflammasome regulation69–71. In addition, 6/10 genes had highly consistent directions of effect. These include FCGRT (p=1.12x-05), which has been characterized as an immune regulator of dendritic cell cross-presentation of IgG immune complexes, necessary to activate a cytotoxic T-cell response and clear antigens72. Our analysis demonstrates a down-regulation of FCGRT in WTC responders with PTSD, suggesting reduced IgG immune complex clearance in these patients (Figure 6).
There is evidence that macro-and micro-level physiological damage is a fundamental component of PTSD, as well as cytoskeletal restructuring for fear-based memory formation in the amygdala43, 44. In this study we observed decreased expression in CIRBP (p= 1.29x-06), a protein that traditionally regulates stress and apoptosis under conditions of extreme cold (Figure 6). Its role as a potential biomarker has been previously explored for different psychiatric disorders73, 74. The decrease observed in TMSB10 expression also contributes to the dysregulation of apoptosis. TMSB10 is a pro-apoptotic protein that has been previously associated with downregulation of gene expression after trauma exposure3, 75, 76. In our meta-analysis, the directionality of effect (p<2.2x-16) for each gene was decreased, consistent with PTSD pathophysiology77, 78.
Since many of our PTSD-associated genes are related to immune function, we tested for the enrichment of immune cell types in our study. We found an overall enrichment of CD4-positive T cells for both past-month and lifetime CAPS scores (Figure 7), consistent with previous studies79, 80, including among WTC responders8. In addition to CD4 T cell enrichment, our study also found enrichment in eosinophils and a decrease in natural killer cells for CAPSL (Figure 7). These cellular diversities may point to a higher inflammatory signature in PTSD, particularly in the case of CD4 T cell enrichment. It has been well demonstrated that dysregulation of CD4 T cells leads to autoimmune activation81, and in combination with an increase in eosinophils can lead to an inflammatory cascade in patients. There is strong evidence that PTSD is associated with a pro-inflammatory state, which our findings support82–86.
While our study provides an in-depth look at the genetic expression and outcomes related to a specific traumatic experience, we note some significant caveats. Our expression analysis was limited to blood, but should be expanded to other tissues in the future, such as the brain. Similarly, our analysis was restricted to whole blood, but a more in-depth single cell analysis will be critical to determine gene expression in individual cell types. In addition, we note that our cohort includes a significant proportion of individuals who have self-selected into high-risk professions. As such, we expect a higher lifetime exposure to stressful situations, including potentially many other life-threatening scenarios. It is likely that the PTSD symptoms observed here are at least partially accounted for by other traumas and stressors, even though upon CAPS administration, study clinicians specifically inquired about WTC-related PTSD symptoms. Conversely, the high-risk nature of these individuals’ occupations may also mean increased exposure to resilience training for a sizable subsample, and greater access to social support networks of peers with similar experiences, potentially providing protective mechanisms.
In conclusion, this study has identified a vast number of biomarkers that will be potentially useful tools after independent validation. In particular, ten of these genes stand out as reproducible across multiple studies and should be considered as high priority. In combination with pathway and cellular deconvolution results, these findings highlight a strong connection with immune dysregulation and other psychiatric illnesses. We believe that future studies should focus on validation of our PTSD-associated genes and also single-cell RNA-sequencing approaches to delineate the role of immune cell types in PTSD.
Data Availability
Gene expression summary statistics will be made available upon request.
Supplemental Table 2. (A) Correlations of beta values and FDR-adjusted p values for CAPS lifetime versus past-month and all associated symptom dimensions. (B) Correlations of symptom dimensions against all other symptom dimensions for past-month and (C) lifetime. All correlations were significant (p<2.2×10-16).
Footnotes
This updated version includes all tables within our main manuscript, instead of as supplemental items, to enable readers to more easily refer to all our data.
References
- 1.↵
- 2.↵
- 3.↵
- 4.
- 5.
- 6.
- 7.
- 8.↵
- 9.
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.
- 41.
- 42.
- 43.↵
- 44.↵
- 45.↵
- 46.
- 47.
- 48.
- 49.
- 50.
- 51.
- 52.
- 53.
- 54.
- 55.
- 56.
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.
- 84.
- 85.
- 86.↵