Single-cell transcriptomic and chromatin dynamics of the human brain in PTSD

Abstract

Post-traumatic stress disorder (PTSD) is a polygenic disorder occurring after extreme trauma exposure. Recent studies have begun to detail the molecular biology of PTSD. However, given the array of PTSD-perturbed molecular pathways identified so far¹, it is implausible that a single cell type is responsible. Here we profile the molecular responses in over two million nuclei from the dorsolateral prefrontal cortex of 111 human brains, collected post-mortem from individuals with and without PTSD and major depressive disorder. We identify neuronal and non-neuronal cell-type clusters, gene expression changes and transcriptional regulators, and map the epigenomic regulome of PTSD in a cell-type-specific manner. Our analysis revealed PTSD-associated gene alterations in inhibitory neurons, endothelial cells and microglia and uncovered genes and pathways associated with glucocorticoid signalling, GABAergic transmission and neuroinflammation. We further validated these findings using cell-type-specific spatial transcriptomics, confirming disruption of key genes such as SST and FKBP5. By integrating genetic, transcriptomic and epigenetic data, we uncovered the regulatory mechanisms of credible variants that disrupt PTSD genes, including ELFN1, MAD1L1 and KCNIP4, in a cell-type-specific context. Together, these findings provide a comprehensive characterization of the cell-specific molecular regulatory mechanisms that underlie the persisting effects of traumatic stress response on the human prefrontal cortex.

Main

Post-traumatic stress disorder (PTSD) is a common mental health disorder that occurs in the aftermath of serious trauma. PTSD occurs in the general population at a rate of approximately 6–8%^2,3, is moderately heritable (ranging from 24% to 40%)⁴ and is highly polygenic^5,6,7,8. Recent studies have begun to elucidate the molecular biology of the post-mortem brain in individuals who had PTSD^9,10,11,12. These studies have predominantly focused on the prefrontal cortex (PFC) and the amygdala and have identified molecular changes in several gene pathways, including GABAergic signalling, immunity and neuroinflammation, and glucocorticoid signalling. Given the highly diverse nature of the affected biological processes, it is unlikely that one particular cell type is responsible for PTSD pathophysiology.

Gene expression in the human brain is orchestrated by a combination of cis-regulatory elements (CREs), chromatin modifications and transcription factor binding to enhancers and promoters of genes. The role of CREs in specific cell types is critically important in understanding how disease risk variants regulate gene expression. Previous work has identified significant changes in DNA methylation in PTSD^12,13; however, to our knowledge, there have been no studies comprehensively examining the single-cell chromatin landscape of the brain in PTSD. Recent advances in genomic technologies now enable the interrogation of chromatin assemblies in individual cells, and these assays can be coupled with gene expression analysis to provide the resolution needed to identify how risk variants regulate transcription within individual cells^14,15,16.

Here we present a multi-omic analysis of more than two million individual nuclei from the dorsolateral PFC (DLPFC) of three diagnostic groups: individuals with major depressive disorder (MDD), individuals with PTSD and control individuals with neither condition (CON). We performed paired single-nucleus RNA sequencing (snRNA-seq) with single-nucleus assay for transposase-accessible chromatin sequencing (snATAC-seq). We annotated and censused all major cell types, including excitatory and inhibitory neurons and non-neuronal cell types, by integrating transcriptomic and epigenomic profiles. We identified both cell-type-specific differentially expressed genes (DEGs) in PTSD and converging and diverging expression changes between PTSD and MDD. We further performed high-resolution single-cell spatial transcriptomics on tissue from 18 donors, confirming DEGs and disruption of cell-to-cell communication (CCC) patterns in somatostatin (SST) interneurons and microglia in PTSD. Next, we constructed the gene expression regulatory landscape of PTSD by integrating RNA and ATAC modalities to define CREs and link them to specific genes and to fine-map eight credible risk loci within these putative enhancers. We discovered selective changes in the glucocorticoid system that were, unexpectedly, most pronounced in endothelial cells. In addition, we identified vulnerability of SST interneurons in PTSD and global shifts in the transcriptome, reflecting decreases in SST signalling and neurotransmission. Overall, this work enabled characterization of gene pathways and their dynamics in diverse cortical cell types and prediction of cis-regulatory logic and associated factors underpinning the molecular aetiology of PTSD.

Single-cell multi-omic analysis of the PFC in PTSD

To dissect the single-cell-type gene expression and regulatory changes in PTSD (Fig. 1a), we performed snRNA-seq and snATAC-seq on the DLPFC of 111 donors. We profiled 36 donors with PTSD, 36 with MDD (psychiatric control) and 39 with neither (CON) (Extended Data Fig. 1a). We report results from a discovery cohort of 935,371 nuclei for snRNA-seq (Fig. 1b and Supplementary Fig. 1), 473,033 nuclei for snATAC-seq (Fig. 1c and Supplementary Fig. 2) and 119,431 nuclei for single-nucleus multiomics (snMultiome) (Fig. 1d and Supplementary Fig. 3). The inclusion of a small snMultiome dataset was used to match barcoded mRNAs directly with chromatin from the same nuclei. Four snMultiome samples overlapped with the discovery dataset (Extended Data Fig. 1b and Supplementary Table 1).

Fig. 1: Multimodal genomic taxonomy of cell types in the human PFC and cell-type-specific gene expression changes.

a, Schematic of the analyses. TF, transcription factor. b, Uniform manifold approximation and projection (UMAP) of snRNA-seq (n = 935,371 nuclei) across 14 subtypes (top) and cell proportion of subtypes across diagnostic conditions (CON, MDD and PTSD; bottom). *FDR < 0.05. c, UMAP of snATAC-seq (n = 473,033 nuclei) across seven cell types (top) and proportion of cell types across conditions (bottom). d, UMAP of snMultiome (n = 119,431 nuclei) across 14 subtypes (top) and proportions across conditions (bottom). e, Significant DEG counts in both directions for each major cell type (left), EXN subtypes (right top) and IN subtypes (right bottom) coloured by the number of DEGs. f, Binary plot indicating occurrence of n = 1,184 PTSD snDEGs across cell types. g, log₂FC values of the top DEGs for each cell type from MAST (top) and bulk RNA-seq (bottom). h, Overlap between n = 1,184 PTSD snDEGs and n = 1,918 MDD snDEGs. i, Top biological process (BP) and molecular function (MF) Enrichr Gene Ontology terms for the n = 502 PTSD-specific DEGs from panel h. j, Log-normalized mean expression of significantly discordant genes (CTNNA3 (top) and HSPA1A (bottom)) in each subtype for PTSD (blue) and MDD (orange). The asterisk indicates significance from rank–rank hypergeometric overlap (RRHO).

We identified excitatory neurons (EXNs; four subtypes: CUX2, RORB, FEZF2 and OPRK1), inhibitory neurons (INs; five subtypes: LAMP5, KCNG1, VIP, SST and PVALB), oligodendrocytes (OLGs), oligodendrocyte progenitor cells (OPCs), endothelial cells (ENDs), astrocytes (ASTs) and microglia (MG; Fig. 1b and Extended Data Fig. 2a) with expression profiles of cell-type-specific marker genes (Extended Data Fig. 2d). The number of nuclei per cell type and dataset is reported in Supplementary Table 2. Differential abundance of cell types was calculated using scCODA¹⁷ based on the number of nuclei isolated across modalities. We found decreases (false discovery rate (FDR) < 0.05) in the number of OLGs in PTSD versus CON and MDD (Fig. 1b). We identified all seven canonical brain cell types using snATAC-seq (Fig. 1c and Extended Data Fig. 2b) with chromatin accessibility profiles of each marker (Extended Data Fig. 2e) and snMultiome (Fig. 1d and Extended Data Fig. 2c). To identify more specific transcriptomic cell subtypes (Extended Data Fig. 2f), we used previous annotations¹⁸ to identify 61 subtypes: 18 EXNs and 30 INs and 13 non-neuronal subtypes, including two OLGs, three ASTs, macrophages, T cells (T SKAP1 CD247), smooth muscle cells and pericytes (Extended Data Fig. 2f,g).

PTSD-specific changes in gene expression

We performed differential gene expression analysis across all 14 cell-type clusters using two common cell-based DEG analyses: MAST¹⁹ (Supplementary Table 3) and Wilcoxon (Supplementary Table 4). In addition, we utilized a mixed linear model of gene expression²⁰, which allows correction of sample-based bias (Supplementary Table 5). Finally, we performed a sample-based pseudobulk analysis using DESeq2 (ref. ²¹) (Supplementary Table 6). Here we define snDEGs as overlapping significant DEGs between MAST and Wilcoxon (Methods). For our DEG analysis, we analysed the canonical cell types and subtypes (Fig. 1e). We found 322 DEGs in EXNs and 66 in INs, with the most (557) occurring in CUX2 EXNs. Across the seven cell types, we identified 3,431 DEGs from MAST and 2,989 from Wilcoxon. Of those, 1,376 DEGs overlapped in the same direction and included 1,184 unique genes (PTSD snDEGs; Fig. 1f and Supplementary Table 7). We found that 9.92% of DEGs were significantly altered in both EXNs and INs (Fig. 1f), and no DEG was regulated across all cell types. Cosine similarity revealed the highest concordance of gene expression within neuronal subtypes (Extended Data Fig. 3a).

Gene Ontology analysis of PTSD snDEGs (Extended Data Fig. 3b,c and Supplementary Table 8) and PTSD cell-type-specific DEGs (Extended Data Fig. 3d and Supplementary Table 9) showed enrichments in ubiquitin binding, cell stress response and cadherin binding. Nebula revealed similar enrichments, including cadherin binding, phosphodiesterase signalling and axonogenesis (Supplementary Table 10) across 437 DEGs. Pseudobulk analysis revealed enrichment for transcription factor binding and oestrogen receptor signalling from additional significant DEGs for EXNs (73) and INs (19), but few in non-neuronal cell types (for example, MGs (8) and OPC (1); Supplementary Table 11). In addition, we observed sex-specific responses to PTSD and MDD, consistent with previous observations⁹ (Supplementary Tables 12–16).

Analysis of high-quality bulk tissue RNA-seq from a partially overlapping cohort (n of approximately 150 samples)⁹ showed consistent direction and magnitude for the top PTSD snDEGs (Fig. 1g), including previously identified genes FKBP5, YBX3 (ref. ⁹) and HDAC9 (ref. ¹¹). Globally, we found high correlation between both datasets (r = 0.69; Extended Data Fig. 3e). We were able to identify many novel subtype-specific genes including GRM5 (refs. ^22,23) in LAMP5 INs and CLU²⁴ in CUX2 and OPRK1 EXNs (Supplementary Table 7), which have been previously implicated in suicide and Alzheimer’s disease. We found increases in FKBP5 expression in ENDs (log₂ fold change (log₂FC) = 1.42, FDR = 3.42 × 10^−128) and in OLGs (log₂FC = 1.02, FDR = 7.63 × 10^−282; Fig. 1g). FKBP5 transcripts have been reported to be upregulated in the post-mortem frontal cortex in PTSD⁹ and inversely associated with mushroom spine density²⁵. We compared our DEG results to a partially overlapping smaller dataset (n = 10–11 per diagnostic group)²⁶ and found moderate overlap in the number and direction of nominally significant DEGs with both our PTSD (61.4%) and MDD (72.7%) cohorts across cell types (Extended Data Fig. 3f,g).

Gene expression comparisons between PTSD and MDD

We included a psychiatric control group with a condition (MDD) that is highly comorbid in patients diagnosed with PTSD (more than 50%)^27,28,29. We systematically compared the upregulated and downregulated genes in PTSD versus CON and MDD versus CON to identify common and divergent molecular mechanisms. We analysed the MDD cohort in the same manner as the PTSD cohort (Methods) and identified 1,918 MDD snDEGs (Extended Data Fig. 4a and Supplementary Table 17). Overall, we found high overlap (57.6%) in snDEGs between MDD and PTSD (Fig. 1h), with 502 PTSD-specific and 1,236 MDD-specific snDEGs (Supplementary Tables 18 and 19). Analysis of the PTSD-specific snDEGs revealed enrichment for pathways related to calmodulin signalling, cadherins and ubiquitin binding (Fig. 1i and Supplementary Table 20). We found high correlation in global transcript levels between snRNA-seq and bulk MDD datasets (r = 0.72; Extended Data Fig. 4b) and moderate correlation between the top MDD DEGs (Extended Data Fig. 4c). We observed the highest degree of overlap in DEGs between MDD and PTSD in non-neuronal cell types (Extended Data Fig. 4d).

To characterize the differences between PTSD and MDD, we examined DEGs that were divergent in expression (that is, up in PTSD and down in MDD, and vice versa) through threshold-free rank–rank hypergeometric overlapping³⁰. We found significantly high enrichment in concordant genes (87.8%) and low enrichment in discordant genes by P value (Extended Data Fig. 4e) and odds ratio (Extended Data Fig. 4f). Twelve genes were significantly discordant (down in PTSD and up in MDD) across cell types. Of note, CTNNA3 (encoding catenin α3) was downregulated in seven out of nine neuronal subtypes in PTSD and upregulated in MDD, and HSPA1A (encoding heat shock protein A1A) was downregulated in PTSD but upregulated in MDD endothelial cells (Fig. 1j and Supplementary Table 7).

Our large sample size allowed us to examine a subset of donors who had PTSD but not depression. We calculated DEGs between PTSD with no depression (n = 8; Supplementary Table 21) and PTSD comorbid with depression (n = 27; Supplementary Table 22) and each versus CON and found the highest proportion of overlapping DEGs between PTSD with depression and PTSD snDEGs (Extended Data Fig. 4g) for most cell types. Although our PTSD with no depression cohort is admittedly underpowered compared with the rest of our study, several findings emerged from this analysis. We found significant changes in BDNF in EXNs (log₂FC = −0.45, FDR = 3.86 × 10^−40), a gene long implicated in the effects of classical³¹ and rapid antidepressant³² actions (Supplementary Table 21). Although the strength of these findings is difficult to conclude with certainty, this preliminary analysis points to an unmet need to collect and analyse tissue from donors with PTSD and without comorbid mood disorders.

Spatial validation of PTSD gene expression changes

We identified many transcript expression changes across all cell types of the DLPFC. To confirm these changes, we used the 10X Xenium platform to generate spatial gene expression data at single-cell resolution for n = 18 samples (4 PTSD, 4 MDD and 10 CON, Supplementary Fig. 4) for a targeted panel of genes (366) selected on the basis of cell-type markers and expression changes in our snRNA-seq analysis (Fig. 2a, Extended Data Fig. 1c and Supplementary Table 23). We analysed our spatial data parsed by both soma (Supplementary Table 24) and nucleus (Supplementary Table 25), and subsequently clustered transcripts using Squidpy³³ and canonical cell-type annotations (Fig. 2a) from our snRNA-seq results (Fig. 2b). We observed 550,520 total cells (scXenium) with a median of 347 total transcripts per cell and 100 unique transcripts (Extended Data Fig. 5a, top). Similarly, we observed 523,314 total nuclei (snXenium) with a median of 87 total transcripts per nucleus and 42 unique transcripts (Extended Data Fig. 5a, bottom). We found significant enrichment for END markers (putative blood vessels) near the pial surface and OLG markers localized to the white matter (Extended Data Fig. 5b). Visualization of cell-type-annotated nuclei and individual transcripts revealed the detailed structure of the DLPFC with high resolution (Extended Data Fig. 5c, top). We did not find major anatomical differences among samples after post-Xenium haematoxylin and eosin staining (Extended Data Fig. 5c, bottom).

Fig. 2: Spatial transcriptomic analysis of the PTSD PFC.

a, Sample description and UMAP of snXenium (n = 523,314 nuclei) across seven cell types (top) and proportion of cell types across conditions (bottom). b, Expression of canonical markers for cell-type annotation (diagonal). The circle radius represents the fraction of cells, and the shade represents mean log-normalized expression. c, PTSD snXenium DEG counts in both directions using MAST. d, Overlap between PTSD snRNA-seq and snXenium MAST DEGs (lighter shade denotes snRNA-seq, and darker shade indicates snXenium). e, Volcano plot showing PTSD versus CON MAST DEGs in END from snXenium. Genes labelled and outlined in black overlap with snRNA-seq END MAST results. Horizontal dashed line indicates FDR < 0.01 and vertical dashed line indicates absolute fold change >1.1. f, Mean FKBP5 log-normalized counts comparing n = 10 CON versus n = 4 PTSD samples by cell type. g, CON slide with FKBP5 transcripts indicated in dark purple and each nucleus in its corresponding cell-type colour (END in light purple). The zoomed-in image in the inset shows a local blood vessel with high FKBP5 expression. One experiment was conducted. h, PTSD slide with FKBP5 transcripts indicated in dark purple and each nucleus in its corresponding cell-type colour (END in light purple). One experiment was conducted.

To validate our snRNA-seq dataset, we explored the expression changes across cases (n = 4 MDD and n = 4 PTSD) versus controls (n = 10 CON) for all cell types of the DLPFC using MAST (Fig. 2c). We observed most gene expression changes in neurons. We compared these results with our snRNA-seq MAST results and found the highest DEG overlap by fold-change direction (53%), with neurons (73% EXNs and 42% INs), OPCs (64%) and ENDs (62%; Fig. 2d). Globally, we found high correlation of both nuclear transcripts of snRNA-seq and snXenium (r = 0.62; Extended Data Fig. 5d) and scXenium (r = 0.61; Extended Data Fig. 5e). The high correlation between snRNA-seq and scXenium suggests that soma transcript levels are accurately reflected in the nucleus.

We observed a significant change for 88 DEGs in ENDs (Fig. 2e), with 31 overlapping with our snRNA-seq DEGs. For example, we found that FKBP5 was highly expressed in PTSD non-neuronal cells, and showed a significant increase (log₂FC = 0.14, FDR < 0.01) in expression in ENDs (Fig. 2f). We observed high FKBP5 expression across the entire DLPFC (Fig. 2g,h) and conspicuously higher levels near blood vessels (Fig. 2g,h). In addition, we highlight changes in PTSD INs, where we observed 102 DEGs with a moderate overlap of 42% with our snRNA DEGs (Extended Data Fig. 5f). One of the most upregulated transcripts that we found was CXCL14 (log₂FC = 1.4, FDR < 0.01; Extended Data Fig. 5g–i), a chemokine expressed in interneurons throughout the brain that downregulates GABAergic transmission between neurons^34,35.

PTSD alters interneuron and microglia communication

To understand the function of expression changes in the PTSD transcriptome, we leveraged empirically derived ligand–receptor pairs to perform a CCC analysis by creating a ligand–receptor communication network based on transcript expression levels (Methods). We found decreased microglial sending patterns in PTSD that were increased in MDD (Extended Data Fig. 6a) and driven by upregulation of the SPP1 transcript in microglia and by neuronal expression of integrin receptors in the SPP1 (osteopontin) pathway (Extended Data Fig. 6b). The osteopontin receptor subunit integrin α4 was differentially regulated in several neuron subtypes in PTSD (up) compared with MDD (down; Supplementary Table 26).

Our findings also suggest differences in neuronal communication in PTSD, which may be caused by alterations in synaptic transmission. Because neurotransmitter levels cannot be measured by gene expression, we used NeuronChat³⁶ to approximate neurotransmitter levels by expression levels of synthesis and transport transcripts. We calculated the difference in neurotransmitter sending cell signalling expression between PTSD and CON samples. We found that, in PTSD, SST INs have significant decreases in sender communication compared with other neuronal cell types (Fig. 3a), which confirms a previous report⁹. We also found that the differential strength of communication from SST INs to every other neuronal cell type was downregulated in PTSD, with additional decreases with ASTs, ENDs and OPCs (Fig. 3b and Supplementary Table 27). In addition, we observed significant downregulation of the SST transcript (Fig. 3c) in INs across the DLPFC (Fig. 3d,e). Noticeably, the most downregulated communication occurs from SST INs to KCNG1 INs. Compared with other inhibitory cell types, SST INs underutilize GABA, with decreases in the expression of the GABA transport transcript SLC32A1 (CON third quartile expression = 0.196 TPM (transcripts per million); PTSD third quartile expression = 0 TPM) in particular, resulting in reduced GABA output to GABRA5, GABBR1, GABRB1 and GABRG1 receptors (Fig. 3f).

Fig. 3: CCC alterations in PTSD.

a, Barplot of the log differential output of all sender cell types. The SST IN was downregulated in its output signalling (log differential output = −0.518). b, Circos plot showing differential strength of all cell-to-cell interactions from SST INs between PTSD and CON. The shading indicates their quantitative values, with higher numbers indicating higher differential strengths. c, Mean SST log-normalized counts comparing 10 CON versus 4 PTSD samples by cell types. d, CON slide with SST transcripts indicated in dark green and each nucleus in its corresponding cell-type colour in a lighter shade (IN in light green). The CON inset zooms into layer 5, showing higher expression of SST in IN nuclei. One experiment was conducted. e, PTSD slide with SST transcripts indicated in dark green and each nucleus in its corresponding cell-type colour in a lighter shade (IN in light green). The PTSD inset zooms into layer 5, where SST expression is decreased. One experiment was conducted. f, Circos plots of GABA–GABRA5, GABA–GABBR1, GABA–GABRB1 and GABA–GABRG1 showing a decrease in differential output strength of all cell-to-cell interactions from SST INs in individuals with PTSD. Edges are coloured by interaction strength in PTSD cells (red denotes stronger, and blue indicates weaker). g, Schematic of the slice electrophysiology experimental setup. ChR2 was selectively expressed in SST IN in the mouse mPFC. IPSCs were evoked with 470-nm light pulses and recorded in mPFC layer 5 pyramidal neurons (Pyr). h, Representative eIPSCs (baseline in black, MRX 016 in red, and after washout in grey) recorded in control versus SPS mice i, Quantification of the percentage of IPSC amplitude change in the GABARα5 antagonist MRX 016 compared with baseline. Two-sided Student’s t-test was used, P = 0.01. n = 10 control samples and n = 18 SPS samples. Lines indicate mean, and error bars denote s.e.m. j, Schematic of the experimental setup. k, Representative traces from control (top) and SPS (bottom) held at 0 mV with Baclofen (a GABA_B receptor agonist) and washing in CGP 55845 (a GABA_B receptor antagonist). CGP 55845 induced a decrease in GABA tonic current in a pyramidal neuron in control but not in SPS. Red dotted lines indicate average baseline holding current before and after CGP 55845 bath application. l, Quantification of the size of GABA tonic current in control versus SPS. Two-sided Welch’s t-test was used, P = 0.02. n = 28 control samples and n = 22 SPS samples. Lines indicate mean, and error bars denote s.e.m. *P < 0.05.

Source data

To validate the finding of reduced communication from SST INs to EXNs through GABRA5, we used a rodent model of traumatic stress, single prolonged stress (SPS)³⁷ (Extended Data Fig. 6c,d), and assessed synaptic function in the mouse PFC (mPFC) using patch-clamp slice electrophysiology. We expressed channelrhodopsin-2 (ChR2) selectively in SST INs and recorded the optically evoked SST IN-mediated inhibitory postsynaptic currents (eIPSCs) in EXNs from control and SPS mice (Fig. 3g). After bath applying the GABA receptor α5 subunit (GABRAα5)-specific antagonist MRK 016, we observed a significant reduction in eIPSC amplitude (Fig. 3h,i), confirming that GABRAα5-containing GABA receptors are major contributors to SST IN-to-EXN eIPSCs. By contrast, the average amplitude decreased slightly after MRK 016 in SPS mice, with the change significantly smaller than control, indicating a significant reduction in GABRA5-mediated SST IN-to-EXN synaptic transmission (Fig. 3i). In addition, we assessed the level of tonic GABA current in control versus SPS (Fig. 3j). Tonic inhibition is mediated by GABA_B receptors (containing GABA_BR1 subunit; encoded by GABABR1) directly and by extrasynaptic GABA_A receptors regulated by GABA_B receptors³⁸. SST INs have been shown to silence excitatory synaptic transmission through tonic activation of GABA_B receptors³⁹. We found that blocking GABA_B receptors with CGP 55845 caused a significantly larger reduction in the level of tonic current in EXNs in control than in SPS mice (Fig. 3k,l). Together, these results support the conclusion that SST IN synaptic transmission to EXNs is reduced in SPS mice.

We computed neurotransmitter–receptor interactions for the top 15 most differential interactions between PTSD and CON (Extended Data Fig. 6e), and found significant, but not complete, reduction of information flow involving glutamate (Glu)–GRM2 and CORT–SSTR2 pathways in PTSD. We also observed significantly reduced signalling in ten additional pairings across both glutamatergic and GABAergic neurons (Extended Data Fig. 6e). We observed increased PTSD communication in CRH (encoding corticotropin-releasing hormone)–CRHR1 signalling (generally from INs to EXNs; Extended Data Fig. 6f), consistent with the role of CRH in the stress response. We highlight the cell-type-specific ligand–receptor activities for several pairs: CRH–CRHR1, Glu–GRM5, CORT–SSTR2 and Glu–GRM7 (P = 0.0039, 1.8 × 10^−12, 0.06 and 1.8 × 10^−12, respectively). We observed the most increased communication for CRH–CRHR1 from VIP INs to CUX2 EXNs, for CORT–SSTR2 from LAMP5 INs to OPRK1 EXNs for Glu–GRM5 from CUX2 EXNs to FEZF2 EXNs, and for Glu–GRM7 from CUX2 EXNs to OPCs (Extended Data Fig. 6f). These findings suggest cell-type-specific gene expression changes altering Glu and GABAergic transmission in the PTSD cortex and provide a possible mechanistic function for PTSD pathophysiology via neuronal communication disruption.

PTSD alters cis-regulation of gene expression

We measured chromatin accessibility using snATAC-seq to identify cell-type-specific cis-gene regulatory mechanisms in PTSD. We found 913,717 chromatin peaks across all cell types, mostly within non-coding regions of the genome, such as introns (51.1%) and intergenic regions (26.2%) (Fig. 4a). Nearly half of all peaks are cell-type specific (Fig. 4b and Extended Data Fig. 7a) with the most overlap between EXNs and INs (57.1%) and between OLGs and OPCs (57.2%; Fig. 4b and Extended Data Fig. 7b). Comparisons with bulk-tissue brain open chromatin regions⁴⁰ (bCREs) from the DLPFC confirmed our findings, with the majority (88.6%) of peaks overlapping (fraction > 0.4) averaged across all seven cell types (Extended Data Fig. 7c). Integrating snATAC-seq with a subset of six high-quality snRNA-seq samples, we found 1,433,145 peak-to-gene links (correlation > 0.45, FDR < 1 × 10^−4) using the union peaks (Fig. 4c). Previous studies have shown that not all peaks control gene expression⁴¹, therefore we define a CRE as a peak with peak-to-gene linkage (correlation > 0.4, FDR < 0.05; Supplementary Table 28), and found that 89.18% of peaks are CREs across cell types (Fig. 4d and Supplementary Table 29).

Fig. 4: PTSD alters cis-regulation of gene expression across cell types.

a, Number of chromatin peaks separated by genomic category (intronic (intron), distal (dist), promoter (prom) and exonic (exon)) for the major cell types in snATAC-seq. b, Upset plot showing the number of peaks (first bar by cell type) and the overlap of peaks across different cell types (subsequent bars). c, Side-by-side heatmaps of linked ATAC (left) and gene (right) regions. Peak-to-gene links (n = 1,433,145) are shown using thresholds: FDR < 1 × 10^−4 and correlation > 0.45. d, Donut plot showing the number of CREs (n = 395,932) with significantly linked genes (correlation > 0.4 and FDR < 0.05) separated by genomic category. e, Venn diagrams showing the intersection of CLGs and PTSD DEGs for each cell type. f, UMAP plot of FKBP5 highlighted by snRNA-seq log-normalized gene expression (top) and snATAC-seq gene score (bottom). Non-neuronal cell types with high expression of FKBP5 are outlined by the dashed line. g, Chromatin accessibility signal tracks highlighting FKBP5 peak-to-gene links (correlation > 0.75) across cell types in region chromosome 6 (Chr. 6): 35200000–35900000. The signals for non-neuronal cell types are indicated by the dashed boxes.

To identify disease gene-regulatory mechanisms, we intersected our CRE-linked genes (CLGs) with DEGs (CLG–DEGs) for each cell type. We found the most overlap in CLGs in ENDs (544) and EXNs (313; Fig. 4e). We also found that in EXNs and ASTs, most CLG–DEGs corresponded with upregulated genes (84% EXNs and 66% ASTs) and that CLG-target DEGs overlapping between cell types were mostly linked to different CREs (Extended Data Fig. 7d). For example, YBX3 is upregulated in two different cell types (ENDs and ASTs; Fig. 1g) each with unique, cell-type-specific CREs (Supplementary Table 29). Similarly, FKBP5 is highly expressed in several cell types (Fig. 4f) and significantly upregulated in non-neuronal cells (Fig. 1g). We found that the peak-to-gene links that correlated highly (correlation > 0.75) with the transcription start site (TSS) of FKBP5 were all ENDs and MG CREs, probably influencing gene expression in these cell types (Fig. 4g, dotted boxes, and Supplementary Table 29).

To parse out the different epigenomic effects of PTSD and MDD on transcription, we called peaks for all cell types by condition (Supplementary Fig. 5a). After removing the overlapping peaks between conditions, and using the same definition of CREs, we found a total of 141,939 condition-specific CREs across all seven cell types and defined them by genic feature (Supplementary Table 30 and Supplementary Fig. 5b). With the PTSD-specific and MDD-specific CREs, we identified PTSD CLG–DEGs and MDD CLG–DEGs across all seven cell types (Supplementary Fig. 5c,d). We then overlapped the PTSD CLG–DEGs and MDD CLG–DEGs to investigate the cell-type-specific epigenomic differences between the two disorders. Globally, we found 383 PTSD-specific CLG–DEGs and 1,063 MDD-specific CLG–DEGs. We found high overlap between CLG–DEGs and snDEGs for both PTSD-specific (86%) and MDD-specific (72%) groups (Supplementary Fig. 5e), suggesting disease-specific chromatin regulation.

We confirmed our findings using snMultiome from a subset of samples. Similar to previous analyses, we found the peak-to-gene links from 14 subtypes for snMultiome (Extended Data Fig. 8a). In addition, we found 774,147 peaks for the seven canonical cell types (Extended Data Fig. 8b) and 1,614,118 peaks for all 14 subtypes (Extended Data Fig. 8c). Consistent with snATAC-seq, we found approximately 90% overlap among snMultiome and bulk ATAC-seq peaks for both cell types and subtypes at a minimum overlap threshold of 0.4 (Extended Data Fig. 8d,e). In addition, Jaccard matrices show high concordance in cell-type-specific peak sets among neuronal cell types, especially the EXNs (Extended Data Fig. 8f,g). We also found high concordance between snATAC-seq peaks and snMultiome peaks across the seven cell types (Extended Data Fig. 8h).

Transcription factor-binding is altered by PTSD

We inferred transcription factor-binding motifs⁴² within our CREs linked to DEGs to understand the upstream regulatory mechanisms governing gene expression. First, we calculated transcription factor motif variabilities in individual cell clusters (Extended Data Fig. 9a). Transcription factor motif enrichment patterns for all cell types revealed many shared patterns among cell types, consistent with previous findings⁴³. We hypothesize that chromatin accessibility differences harbour transcription factor-binding sites upstream of PTSD snDEGs. To test this, we probed for enrichment of the top 20 transcription factors with their binding sites in the cell-type-specific CREs (TF–CRE–gene) and identified all DEGs regulated by them for each cell type (Supplementary Tables 31–37). We found that EGR1 regulated the most DEGs among the top transcription factors enriched in EXNs (Supplementary Table 38) and a high degree of overlap among all neuronal transcription factor-regulated gene sets (Extended Data Fig. 9e). We highlight two transcription factors with high enrichment in EXNs (Extended Data Fig. 9b) and INs (Extended Data Fig. 9f) that are upstream of DEGs and PTSD genetic risk genes.

To further investigate binding dynamics of these cell types, we performed transcription factor footprinting analysis for two transcription factors. EGR2 revealed stronger enrichment in PTSD EXNs than in CON EXNs (Extended Data Fig. 9c), consistent with its upregulation in CUX2 EXNs (log₂FC = 0.28, FDR = 1.51 × 10^−65). In EXNs, EGR2 was linked to 154 DEGs and several PTSD risk genes including CAMKV, EGR3 (upregulated), KCNB2 (downregulated) and OPCML (Extended Data Fig. 9d and Supplementary Table 31). Similar to the other top transcription factors in EXNs (Extended Data Fig. 9e), EGR2 is an immediate early gene. We also found significant enrichment of motifs for the neuron-specific chromatin modifier SMARCC1, which was upstream of 168 PTSD snDEGs and 26 PTSD risk variants in EXNs. To confirm downstream changes of our PTSD transcription factor networks, human induced pluripotent stem cell-derived excitatory neurons were infected with lentiviral RNA interference constructs of either EGR2 or SMARCC1 to assess transcript changes (Supplementary Table 39). We confirmed expression changes of 27 DEGs including CLU and UST (Extended Data Fig. 9d, pink circles). TFAP4 was significantly enriched in INs and OPCs (Extended Data Fig. 9a,f). In INs, TFAP4 was linked to several PTSD genome-wide association study (GWAS) genes including MAD1L1 and 36 DEGs including ANKRD55, SPRED2 and PDE10 (all downregulated), consistent with its role as a transcription repressor⁴⁴ (Extended Data Fig. 9g,h and Supplementary Table 32).

Cell-type-specific cis-regulation at PTSD risk loci

To assess the role that genetic risk variants have in specific cells of the brain, we performed linkage disequilibrium score regression (LDSC)⁴⁵ for PTSD, its clinical symptoms (hyperarousal, avoidance and re-experiencing), MDD, additional psychiatric disorders, and other clinical and non-clinical traits GWAS (Fig. 5a and Supplementary Table 40). We found significant enrichment (FDR < 0.05) for PTSD GWAS signals in ATAC peaks of both EXNs and INs for PTSD (total PCL) and in the PTSD case–control groups. These findings were mirrored in the symptom GWAS as well. We also found significant enrichment (FDR < 5 × 10^−4) for MDD risk in these same cell types. We compared these findings to the bCREs (grey dots) and found higher LDSC enrichment across all traits (Extended Data Fig. 10a). To confirm, we performed the same analysis using snMultiome to identify PTSD risk gene enrichments within neuronal subtypes. Comparing snMultiome with the same bCREs (in grey), we found higher LDSC enrichment (Fig. 5b). We identified an EXN subtype (OPRK1; −log₁₀P = 6.72) for the total PCL GWAS and an IN subtype (KCNG1; −log₁₀P = 6.15) for the PTSD case–control, suggesting possible differences in the two populations.

Fig. 5: Cell-type-specific cis-regulation at PTSD disease genetic risk loci.

a, LDSC enrichment of various GWAS traits including PTSD, other psychiatric and non-psychiatric disorders in the snATAC-seq cell types (*FDR < 0.05, **FDR < 0.005 and ***FDR < 0.0005). IBD, irritable bowel disease. b, Lollipop plot showing LDSC enrichment of GWAS traits comparing bCREs (grey) versus snMultiome peaks (coloured). The colour of the dot represents the subtype in snMultiome with the highest enrichment value for the trait. c–f, Cis-regulatory architecture at the following GWAS loci and cell types: ELFN1 in IN; and KCNIP4, EGR3 and LRFN5 in EXN. IN and EXN CREs were used to fine-map SNPs for total PCL⁶. For the fine-mapping track, we only plotted SNPs with PIP > 0.0001 for visualization purposes, and only highlight peak-to-gene loops for SNPs with PIP > 0.05, outlined in black if they are within a CRE (lead SNP in red, and MVP SNP in purple). The top two credible SNPs and MVP SNP are labelled. The loops are coloured by cell type (EXN in red and IN in green). Chromosome coordinates are: ELFN1 (chr. 7: 1450000–2350000; c); KCNIP4 (chr. 4: 20600000–22600000; d); EGR3 (chr. 8: 22193302–23193302; e); and LRFN5 (chr. 14: 41105000–4230000; f). NeuN⁺ TADs are shown above each ideogram, and the coloured triangles correspond to the TAD with significant inter-chromosomal looping in the area.

Next, we fine-mapped PTSD-associated polymorphisms within cell-type-specific CREs in EXNs and INs (Methods). We analysed all risk loci from two previous GWAS^5,7 and identified eight that fell within cell-type-specific CREs and focused on these for fine-mapping (Fig. 5c–f and Extended Data Fig. 10c–h). We evaluated the posterior inclusion probability (PIP) for the PTSD risk single-nucleotide polymorphisms (SNPs) and found high correlation of PIP values versus negative log-transformed GWAS P values (mean correlation = 0.769) of SNPs that lie within CREs for 7 out of 8 genes, suggesting that fine-mapping captured significant lead SNPs with strong GWAS associations (Extended Data Fig. 10b). For the gene MAD1L1, we found that the leading causal variant for total PCL was rs10235664 (PIP = 0.01), identified by Million Veteran Program (MVP), in INs (Fig. 5c). However, after incorporating IN CREs, we detected additional variants with higher PIP scores. Of note, the variant rs62444919 (PIP = 0.36) exhibited a high degree of linkage disequilibrium with the MVP significant variant rs10235664 (correlation = 0.84) and was also predicted to loop to the TSS of the gene ELFN1. The peak-to-gene loops have high-correlation values of 0.87 (FDR = 4.29 × 10^−29) and 0.67 (FDR = 2.29 × 10^−13), suggesting strong enhancer activity of ELFN1 transcription. We also found a topologically associated domain (TAD) in this region using a NeuN⁺ Hi-C dataset track, confirming chromosomal looping between ELFN1 and MAD1L1. In addition, our analysis was able to confirm a previously identified causal SNP for KCNIP4, rs4697248 (PIP = 0.93), in EXNs (Fig. 5d). We highlight fine-mapped risk SNPs for six PTSD genes (EGR3, LRFN5, OPCML, CAMKV, TCF4 and CRHR1) in EXNs and INs, with their corresponding CREs (Fig. 5e,f, Extended Data Fig. 10c–h and Supplementary Table 41). By incorporating cell-type-specific ATAC data into the PIP score and identifying TADs in these loci, we have identified new putative causal risk variants and highlight candidates for additional functional genomic validation.

Discussion

Here we report the findings from over 2 million nuclei isolated from the PFC of the human brain. To date, this is the largest single-cell genomics analysis of a human brain subregion to detail the biology of PTSD and MDD. This dataset enabled us to identify genes and pathways associated with PTSD pathology including glucocorticoid^26,46, immune and neuroinflammatory mechanisms¹⁰, and neurotransmitters, especially GABA⁹. By tracking differences between cases and controls, we were able to pinpoint selective vulnerability of SST neurons with marked depletion of outgoing synaptic signalling to other neuronal cells. We also identified distinct differences between PTSD and MDD focused on inflammatory pathways of microglia and cell adhesion in neurons (Extended Data Fig. 6). Through cell communication analyses, we found increased signalling from MGs in MDD specifically through the pro-inflammatory SPP1 pathway that were decreased in PTSD, suggesting suppression of neuroimmune processes and microglial activity in the PTSD brain that are not occurring in depression.

We identified more robust gene expression changes in endothelial cells than any other cell type and found significant upregulation of FKBP5, a possible risk gene for PTSD⁴⁷ (Fig. 1g). Several studies have emphasized the role of FKBP5 in the stress response within neuronal cell types⁴⁸, but our finding and other single-cell transcriptomic studies of the PFC⁴⁹ have consistently identified high levels of FKBP5 transcript in endothelial cells relative to neurons, suggesting a yet unknown neurovascular role. One possibility is that FKBP5 changes reflect changes in peripheral glucocorticoids movement from the periphery to the parenchyma and its upregulation is a compensatory mechanism for the hypocortisolic state commonly reported in PTSD^50,51. In this scenario, more glucocorticoids would need to cross the blood–brain barrier where they would be taken up by glial cells, necessitating an increase in FKBP5 transcript before reaching destination neurons (Fig. 1g).

We identified sets of genomic peaks that probably regulate gene expression changes in PTSD. Although CREs can be identified from epigenomic (ATAC) data alone, integration of single-nucleus transcriptomic data allows us to link gene expression to candidate regulatory regions with open chromatin. We found a high degree of overlap between CREs and corresponding gene expression changes across all cell types (Fig. 4e). We have highlighted cis-gene mechanisms and linked them to cell-type-specific gene expression changes. Our work indicates neurons (both excitatory and inhibitory) as the central cell types harbouring genetic risk for PTSD and MDD. By overlaying chromatin accessibility maps with GWAS statistics of SNP locations, we have unraveled the cis-regulatory relationships disrupted by credible disease variants that map to MAD1L1, KCNIP4, EGR3, CRHR1, TCF4 and LRFN5 (refs. ^5,7). In the case of MAD1L1, we found that causal variants within peaks were putative enhancers for the neighbouring gene ELFN1. ELFN1 was previously identified as an interneuron-specific synaptic protein whose transcript expression was altered in the PTSD frontal cortex⁹. CCC analysis of receptor–ligand expression patterns revealed significantly decreased output of SST INs where ELFN1 expression is highest. We have described decreases in SST IN output to most major cell types with the strongest decreases observed with other neuronal cell types. In particular, we found decreased GABAergic signalling to major GABA receptors including GABRA5, GABBR1, GABRB1 and GABRG1 (Fig. 3f) and verified these findings in vivo (Fig. 3h,k). We observed increased glutamatergic transmission through mGluR7, a presynaptic binding partner of ELFN1 that suppresses synaptic release⁵². Elevation of mGluR7 transmission is consistent with reduced excitatory drive to SST INs and subsequent reduction in their output. Given that ELFN1 is tightly coupled with mGluR7 at the EXN–IN synapse, we speculate that disruption in ELFN1 function may lead to dysregulated mGluR7 signalling and impaired ability to fine-tune excitatory–inhibitory balance in the circuit. Together, these findings suggest PTSD-specific effects at this gene locus that may impart SST interneuron vulnerability and reduced GABAergic communication.

This study provides novel insight into the biology of PTSD and MDD but is limited by several factors. We were able to assemble a demographically well-balanced cohort, matched for sex, age, post-mortem interval and psychiatric disease. We believe that the inclusion of MDD as a psychiatric control is the best possible, as more than 50% of patients with PTSD are diagnosed with some form of depression²⁹, and our study allowed us to identify convergent and divergent molecular effects of both diseases. However, we lack a well-powered PTSD group with no depression phenotype (n = 8), limiting our ability to identify PTSD-specific molecular changes. As is the case with most post-mortem brain studies, there is also a moderately high rate of illegal drug use and drug-related death in both patient cohorts, which could also obscure PTSD-specific or MDD-specific molecular effects. Future efforts will be directed at increasing the number and types of donors to account for these effects.

Methods

Human tissue donors

Human post-mortem brain tissue samples were obtained from the National Center for PTSD Brain Bank (with consent of the next of kin) and the University of Pittsburgh Tissue Donation Program. Individuals were a mix of European and African-American descent. Male and female individuals were group matched for age and post-mortem interval (PMI). Sociodemographic and clinical details are listed in Extended Data Fig. 1a and include comorbidities such as the use of tobacco and substances of abuse. Inclusion criteria for PTSD, MDD and control cases were as follows: PMI < 48 h, age range of more than 18 to less than 85 years. A total of 111 individuals (36 PTSD: 19 female and 17 male individuals; 36 MDD: 18 female and 18 male individuals; and 39 healthy controls: 18 female and 21 male individuals) were used in this study. The brain tissue was fresh-frozen and samples from the DLPFC (Broadman areas 9/46)⁵³ were collected (approximately 25 mg) from each post-mortem sample.

Psychiatric history and demographic information were obtained by psychological autopsies performed post-mortem, as well as a review of medical records and toxicology reports. Trained clinicians conducted structured interviews with informants (usually the next of kin) with knowledge of the deceased individuals. To avoid systematic biases, PTSD, MDD and control cases from each source were characterized by the same psychological methods. Consensus DSM-IV diagnoses for each participant were made by trained clinicians who did not conduct the psychological autopsies. Of the 36 individuals in the PTSD cohort, 75% (27 cases) were also comorbid for MDD.

Isolation of nuclei from brain cells for RNA, ATAC and multiome assays

Regions of interest were dissected on cryotome (leaflets of approximately 100–300 µm) from frozen DLPFC and stored at −80 °C. Cell nuclei isolation from brain sections were treated similarly to already established protocols^18,54,55 with some modifications needed to have nuclei suitable for two separate assays: gene expression and ATAC. To avoid experimental bias, nucleus isolation was done by the same person blinded to the metadata for the particular sample. All reagents were molecular biology grade and sourced from Sigma unless stated otherwise. Small amounts of tissue (10–20 mg) were added into 1 ml of ice-cold lysis buffer (‘buffer A’ is 250 mM sucrose, 25 mM KCl, 5 mM MgCl₂, 10 mM Tricine-KOH (pH 7.8), protease inhibitors without EDTA (Roche), RNAse inhibitor (80 U ml^−1; Roche), 1 mM dithiothreitol, 1% BSA (m/v; Gemini Bio-Products), 0.3% NP-40 (v/v), 0.15 mM spermine, 0.5 mM spermidine and water; dithiothreitol, RNAse Protector, protease inhibitors, spermine, spermidine and NP-40 were added immediately before use. The suspension was transferred to a 2-ml Dounce tissue homogenizer (Kimble) and lysed with constant pressure and without introduction of air with pestle A (2 × 5) and pestle B (4 × 5) in cycles of 5. The homogenate was strained through a pre-wetted 40-μm tube top cell strainer (Corning). All subsequent centrifugation was performed in a refrigerated, bench-top centrifuge with swing-out rotor (Eppendorf). Lysates were centrifuged at 1,000g for 10 min at 4 °C, pellets were saved and resuspended in 0.4 ml lysis buffer buffer A. The final 0.4 ml of solution was mixed with 0.4 ml (1:1) of Optiprep solution (buffer B was iodixanol 50% (v/v), 25 mM KCl, 5 mM MgCl₂, 20 mM Tricine-KOH (pH 7.8), protease inhibitors without EDTA, RNAse inhibitor (80 U ml^−1), 1 mM dithiothreitol and water. The suspension (25% iodixanol final) was mixed 10× head over tail. The 2-ml tube was filled with 0.6 ml of 40% iodixanol cushion (appropriate mix of buffer A and buffer B), then overlayed with 0.6 ml of 30% iodixanol (appropriate mix of buffer A and buffer B), and sample suspension was overlayed at the top. The tubes were then centrifuged at 3,000g for 30 min at 4 °C without brakes. After centrifugation ended, the interphase 30–40% ring was collected. Out of 300 µl collected, 100 µl was used for gene expression studies and 200 µl was used for ATAC studies. One thousand microlitre ‘RNA wash buffer’ was added to gene expression samples, and samples were mixed and left on ice for 10 min. RNA wash buffer was: RNase Protector (80 U ml^−1), 1 mM dithiothreitol, 25 mM KCl, 5 mM MgCl₂, 20 mM Tris-HCl (pH 7.5), 1% (m/v) BSA, 0.1% (v/v) Tween-20 and DPBS (14190, Gibco). Five hundred microlitres of ATAC lysis buffer was added to the ATAC sample, mixed well and left on ice for 5 min. The ATAC lysis buffer was: 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl₂, 0.01% (v/v) Tween-20, 0.01% (v/v) NP-40, 0.001% (v/v) digitonin, 1% (m/v) BSA and nuclease-free water. After 5 min, 500 µl of ATAC wash buffer was added to ATAC sample, mixed well and left on ice for 5 min. ATAC wash buffer was: 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl₂, 0.1% (v/v) Tween-20, 1% (m/v) BSA and nuclease-free water. Both samples were centrifuged at 1,000g for 10 min at 4 °C, supernatants were carefully and completely removed. Gene expression sample was counted on a haemocytometer, resuspended in ‘RNA run buffer’ (RNase Protector (80 U ml^−1), 1 mM dithiothreitol, 25 mM KCl, 5 mM MgCl₂, 20 mM Tris-HCl (pH 7.5), 250 mM sucrose and water) and concentration adjusted to 1 M per millilitre. The ATAC sample was counted on a haemocytometer and resuspended to 3.3 M per millilitre in ATAC run buffer (1× Nuclei Buffer (PN2000153, 10X Genomics) and water). The samples were brought to the Yale Center for Genome Analysis (YCGA) within 15 min and processed by following protocols for 10X Genomics Chromium Single Cell ATAC (CG000168 Rev B, 10X Genomics), and 10X Genomics Chromium Single Cell 3′ Reagent Kits v3 (CG000183 Rev C, 10X Genomics), respectively, for targeting 10,000 nuclei. For some assays that were done with 10X Genomics Chromium Next GEM Single Cell Multiome ATAC + Gene Expression (CG000338 Rev A, 10X Genomics) the 300 µl sample post-3,000g centrifugation was used in its entirety, the steps for ATAC were followed, and the buffers for ATAC had 80 U ml^−1 RNAse Protector. After the last centrifugation, the samples were resuspended in ‘diluted nuclei buffer’ (1× Nuclei Buffer (PN-2000207), 1 mM dithiothreitol, 1 U ml^−1 RNase Protector and water) counted on a haemocytometer at the concentration of 3.3 M ml^−1 and submitted to the YCGA where the protocol was followed for targeting 10,000 nuclei. Libraries were sequenced with paired-end 150-bp reads on an Illumina NovaSeq 6000 to a target depth of up to 500 million read pairs per sample for gene expression (RNA), and up to 500 million reads per sample for ATAC.

Xenium in situ transcriptomics

Gene panel design

The Xenium In Situ platform uses targeted panels to detect expression genes. Genes (n = 267) are included on the Xenium Human Brain Panel, and an additional 99 genes were selected and curated primarily based on our single-cell atlas of DEGs for the PTSD brain (366 genes total). The complete list of genes in our Xenium panel can be found in Supplementary Table 23.

Xenium sample preparation

Fresh-frozen tissue samples from control (n = 10), PTSD (n = 4) and MDD (n = 4) participants were obtained from the National PTSD Brain Bank. Frozen blocks were dissected from the DLPFC from the superior frontal gyrus approximately 1 cm anterior to the genu of the corpus callosum across all participants. Blocks were sectioned at 10 mm thickness and placed on a 10X Genomics Xenium slide as per protocol (10X Genomics protocol CG000580 Rev B). The Xenium slide with the thaw-mounted tissue was fixed in paraformaldehyde solution and then permeabilized in 70% methanol. Tissue slides were washed and then inserted into 10X Genomics Xenium slide cassettes and then processed through hybridization, ligation and amplification protocols (10X Genomics protocol CG000582 Rev C). In brief, the pre-designed Xenium Human Brain Gene Expression panel with 266 genes along with 99 custom probes were added to the tissue. Each circularizable DNA probe contained two regions that hybridized to target RNA and a third region that encoded a gene-specific barcode. Probes were loaded at 10 nM at 50 °C for overnight hybridization. For ligation, the two ends of the probes bound the target RNA and were ligated to generate a circular DNA probe at 37 °C for 2 h. Following ligation, the circularized probe was amplified for 2 h at 30 °C, producing multiple copies of the gene-specific barcode for each target.

Xenium imaging

Prepared tissue slides were then loaded for imaging on the Xenium Analyzer for in situ analysis. Fluorescently labelled oligos were bound to the amplified DNA probes. Cyclical rounds of fluorescent probe hybridization, imaging and removal generated optical patterns specific for each barcode, which were converted into a gene identity. Identified transcripts were then visualized using Xenium Explorer software before bioinformatics analysis.

Animals

All animal care and procedures were performed according to the ethical guidelines of the Institutional Animal Care and Use Committee at Yale University School of Medicine and Yale Animal Resources Center. Adult 2-month-old male and female C57BL/6 SST–Cre mice (n = 31 mice in total) were used for this study. Until stress, all mice were group housed and maintained in standard environmental conditions (23 °C; 12 h–12 h light–dark cycle) with ad libitum food and water.

Viral injections

Adult mice were anaesthetized using 1–2% isoflurane, placed in a stereotaxic apparatus, and maintained at 37 °C during surgery. Each mouse received a unilateral microinjection of AAV-Ef1a-DIOhChR2(123T/T159C)-eyfp (500 nl; Addgene) using a Nanoject III (Drummond) at 50 nl min^−1 with a 5–10-min diffusion period before removing the glass capillary. The coordinates used for mPFC were anteroposterior (AP) = +1.5 mm, mediolateral = +0.2 mm, and dorsoventral = +2.4 mm. Two weeks were allowed for recovery and viral expression after injections.

SPS

Mice for SPS were handled for 3 days before the stress protocol. On the day of SPS, mice (n = 16) underwent a 2-h restraint in a restrainer followed by forced swim for 20 min in 23–25 °C water. After 15 min of recuperation, mice were exposed to diethyl ether vapour until loss of consciousness. Mice were then single housed until euthanized for slice electrophysiology (7–14 days after SPS). Control mice (n = 15) were group housed throughout the experiment. All mice were weighed daily after SPS. Animals were randomly assigned to control or SPS groups. Analysis of behavioural data was performed blinded to condition. Behavioural data have been provided in Source Data File 1.

Slice electrophysiology

Brain slices containing the mPFC were prepared from control and SPS-treated male and female mice as previously described^56,57. After decapitation, brains were removed rapidly and placed in ice-cold (approximately 4 °C) artificial cerebrospinal fluid (ACSF) in which sucrose (252 mM) was substituted for NaCl (sucrose-ACSF) to prevent cell swelling. Coronal slices (300 μm) containing the mPFC were cut in sucrose-ACSF with an oscillating-blade vibratome (Leica VT1000S). Slices were then transferred to the fixed stage of an Olympus BX50WI microscope for whole-cell recordings. The chamber was continuously perfused with normal ACSF at a flow rate of 2–3 ml min^−1, with the temperature maintained at 33 ± 0.5 °C. The standard ACSF (approximately pH 7.35) was equilibrated with 95% O₂/5% CO₂ and contained: 128 mM NaCl, 3 mM KCl, 2 mM CaCl₂, 2 mM MgSO₄, 24 mM NaHCO₃, 1.25 mM NaH₂PO₄ and 10 mM d-glucose. Pyramidal neurons were visualized using a microscope (×40 IR lens) with infrared differential interference contrast (IR/DIC). Low-resistance patch pipettes (3–5 MΩ) were pulled from patch-clamp glass (Warner Instrument) using a horizontal micropipette puller (P-1000, Sutter Instrument). Patch pipettes (3–5 MΩ) were filled with a solution contained the following: 125 mM Cs gluconate, 10 mM HEPES, 5 mM BAPTACs, 10 mM Na2-phosphocreatine, 2.38 mM CaCl₂, 4 mM Mg-ATP and 0.3 mM Na2GTP pH 7.33. Whole-cell recordings were performed using an Axoclamp-2B amplifier (Axon Instruments). The output signal was low-pass-filtered at 3 kHz, amplified 100× and digitized at 15 kHZ and acquired by using Clampex 10.5/Digidata 1550A software.

Tonic GABA current was recorded in voltage-clamp configuration by holding the membrane potential at 0 mV. The difference between holding currents in the presence of the GABA_B receptor agonist baclofen (10 μΜ; Tocris Bioscience) and after the application of the GABA_B receptor antagonist CGP 55845 (20 µM; Tocris Bioscience) was measured. Cells were held for 8 min before recording baseline for 2 min, followed by 5–8 min of recording in each drug application. One minute of recorded traces at plateau of each condition was analysed. Holding current was measured by using a Gaussian fit of the histogram of the 1-min trace.

For experiments measuring eIPSCs by stimulating SST interneurons optogenetically, responses in pyramidal cells were recorded when blue light (470 nm) was applied using CoolLED pE-800 through the ×60 objective. The intensity was set to obtain optimal single-peak responses (approximately 300 pA at baseline, 0.1-ms stimulation). Ten traces were recorded and averaged for light responses. Average eIPSC peak amplitudes before and after bath applying MRK 016 (1 μΜ; MedChemExpress) were measured to calculate the percentage change in IPSC amplitude. Electrophysiology data have been provided in Source Data File 1.

iPS cell-derived glutamatergic neuron transfections

Control human induced pluripotent stem (iPS) cell-derived NPCs (neural progenitor cells) from two lines: NSB553-S1-1, (male, European ancestry) and NSB2607-1-4 (male, European ancestry, source: National Institute of Mental Health) were used. Both lines were mycoplasma negative at the time of experimentation. NPCs were cultured in hNPC medium (DMEM/F12 (10565, Life Technologies), 1× N2 (17502-048, Life Technologies), 1× B27-RA (12587-010, Life Technologies), 1× antibiotic–antimycotic and 20 ng ml^−1 FGF2 (Life Technologies)) on Geltrex (A1413301, Thermo Fisher). hiPS cell–NPCs at full confluence (1–1.5 × 10⁷ cells per well of a six-well plate) were dissociated with Accutase (Innovative Cell Technologies) for 5 min, spun down (5 min at 1,000g), resuspended and seeded onto Matrigel-coated plates at 3–5 × 10⁶ cells per well. Medium was replaced every 2–3 days for up to 7 days until the next split.

At day −1, confluent hiPS cell–NPCs were transduced with rtTA (Addgene 20342) and NGN2 (Addgene 99378) lentiviruses. At day 0, 1 µg ml^−1 dox (doxycycline) was added to induce NGN2 expression. At day 1, transduced hiPS cell–NPCs were treated with antibiotics to select for lentiviral integration (1 mg ml^−1 G-418). At day 3, NPC medium was switched to neuronal medium (Brainphys; 05790, Stemcell Technologies), 1× N2 (17502-048, Life Technologies), 1× B27-RA (12587-010, Life Technologies), 1 µg ml^−1 natural mouse laminin (Life Technologies), 20 ng ml^−1 BDNF (450-02, Peprotech), 20 ng ml^−1 GDNF (450-10, Peprotech), 500 µg ml^−1 dibutyryl cyclic-AMP (D0627, Sigma) and 200 nM l-ascorbic acid (A0278, Sigma) including 1 µg ml^−1 dox. Of the medium, 50% was replaced with fresh neuronal medium once every second day.

On day 5, young hiPS cell–NPC NGN2 neurons were replated onto Geltrex-coated plates. Cells were dissociated with Accutase (Innovative Cell Technologies) for 5–10 min, washed with DMEM, gently resuspended, counted and centrifuged at 1,000g for 5 min. The pellet was resuspended in neuron medium, and cells were seeded at a density of 1.2 × 10⁶ per well of a 12-well plate. At day 11, isogenic glutamatergic neurons were treated with 200 nM Ara-C (C6645, Sigma) to reduce the proliferation of non-neuronal cells in the culture, followed by half medium changes. At day 17, Ara-C (cytosine β-d-arabinofuranoside hydrochloride) was completely withdrawn by a full medium change while adding medium containing four short hairpin RNA (shRNA) vectors targeting a transcription factor (Sigma-Aldrich; targeting SMARCC1: TRCN0000015632, TRCN0000015631, TRCN0000015630 and TRCN0000015628; targeting EGR2: TRCN0000013841, TRCN0000013840, TRCN0000013839 and TRCN0000013838; total multiplicity of infection = 1.0) or a multiplicity of infection-matched scramble shRNA control vector SHC016 (Sigma-Aldrich). Medium was switched to non-viral medium 4 h post-infection. At day 19, transduced isogenic glutamatergic neurons were treated with corresponding antibiotics to the shRNA lentiviruses (1 μg ml^−1 puromycin) followed by half medium changes until neurons were harvested in TRIzol reagent (Invitrogen) at day 21 and bulk sequenced at the Yale Center for Genome Analysis for RNA integrity number and concentration, which were assessed using a Bioanalyzer (Agilent). Libraries were constructed using the SMARTer Stranded RNA-seq Kit (Takara Bio) preceded by rRNA depletion using 1 µg of total RNA. Samples were barcoded for multiplexing and sequenced at 75-bp paired-end on an Illumina HiSeq4000. Samples were pooled eight per lane and sequenced at a depth of 50 million reads.

snRNA-seq processing and analysis pipeline

The count matrix was generated by aligning reads to the hg38 genome using 10X Genomics CellRanger⁵⁸ (v6.1.1). For ambient RNA removal, to further eliminate technical artefacts and background noise in each sample, we used the remove-background command from CellBender (v0.2.2)⁵⁹, which outputs an HDF5 file in which the empty droplets were filtered out. For doublet detection, doublets were identified using a combination of two computational methods — Scrublet⁶⁰ (v0.2.3) and DoubletDetection⁶¹ (v4.2) — and removed from each sample. For sample aggregation, we clustered each sample in Pegasus⁶² (v1.5.0), a Python tool for analysing transcriptomes of single cells, then manually inspected all uniform manifold approximation and projection (UMAP) embeddings and removed six low-quality samples that displayed low heterogeneity of cell types, resulting in a total of 105 snRNA-seq samples. After aggregating these samples into a single data object, we filtered cells based on the following criteria: at most 10% mitochondrial genes, at least 200 genes and at least 500 unique molecular identifiers. Mitochondrial, sex and ribosomal genes were excluded, and only the robust genes were included in the final gene set. The final snRNA-seq data object had 935,371 cells and 27,982 genes. For clustering and annotation, after cell and gene filtration, the data were log-normalized. Next, the top 2,000 highly variable features were found, and principal component analysis was applied using the top 50 components. Batch correction was performed with Harmony. The top 50 components were then used to build a k-nearest-neighbours graph with k = 100 neighbours, and Leiden clustering was used to identify cell clusters. After the first round of clustering, marker genes were used to inspect and remove clusters that had mixed or unmatched marker gene expression, and another round of clustering was performed. Cell-type annotation was adapted using markers from Ma et al.¹⁸ with three levels of hierarchy: (1) 7 cell types, (2) 14 subtypes, and (3) 61 transcriptomic subtypes. The seven cell types include: EXN, IN, OLG, OPC, END, AST and MG. The 14 subtypes include: CUX2, RORB, FEZF2, OPRK1, LAMP5, KCNG1, VIP, SST, PVALB, OLG, OPC, END, AST and MG. Subclustering was performed on each of the 14 subtypes using markers from Ma et al.¹⁸ to identify the gene label for each subcluster, resulting in 61 distinct transcriptomic subtypes at the finest clustering level (Extended Data Fig. 2d,e).

snATAC-seq processing and analysis pipeline

The peak count matrix was generated by aligning reads to the hg38 genome using 10X Genomics CellRanger ATAC¹⁴ (v2.0.0). For quality control and filtering, using the outputs of CellRanger, we created Arrow files for each sample with cells filtered based on the quality control parameters filterTSS=4 and filterFrags=1000 in ArchR⁶³ (v1.0.2), an R package for analysing scATAC-seq data. ArchR considers three characteristics: (1) the number of unique nuclear fragments, (2) the signal-to-background ratio, and (3) the fragment size distribution. Owing to nucleosomal periodicity, we expected to see depletion of fragments that are the length of DNA wrapped around a nucleosome, around 147 bp. We used the addDoubletScores function to infer potential doublets that can confound downstream results. For sample aggregation, after manually inspecting the quality control plots, we removed several low-quality samples, the matching ATAC samples of the six low-quality RNA samples from above, and samples that did not have a matching RNA sample, which resulted in a total of 94 samples. We integrated these samples into an ArchR project, on which all downstream analyses for snATAC-seq were performed. The filterDoublets function was used to filter doublets for downstream analyses. For clustering and annotation, the addIterativeLSI function was used to implement iterative LSI dimensionality reduction. Then, addClusters, which uses Seurat’s graph clustering as the default clustering method, was used to call clusters in the reduced dimension subspace. We visualized the embedding using addUMAP and plotEmbedding. Then, we annotated the clusters using marker genes from Lake et al.⁶⁴. After inspecting the marker genes, we manually cleaned the UMAP and repeated this process twice to produce a final embedding of 473,321 cells annotated into seven cell types: EXN, IN, OLG, OPC, END, AST and MG.

snMultiome processing and analysis pipeline

For count matrix generation, we ran CellRanger ARC (v2.0.2) with the hg38 genome as the reference to process the initial reads from the raw FASTQ files. For quality control and filtering, we used Signac⁶⁵ (v1.11.0), a comprehensive R package for the analysis of single-cell chromatin data, to preprocess each snMultiome sample. Signac creates a Seurat object containing both RNA and ATAC assays. For quality control, we filtered cells using the following parameters: nCount_ATAC > 1,000, nCount_RNA > 200 and TSS.enrichment > 2. Then, we filtered doublets in the snRNA-seq data using DoubletDetection. For sample aggregation, we clustered each sample in Signac and aggregated a total of 25 snMultiome samples in the final data object after manually inspecting each UMAP embedding. For clustering and annotation, we normalized the gene expression data using SCTransform and reduced the dimensionality using principal component analysis. For chromatin accessibility data, we performed latent semantic indexing with functions FindTopFeatures, RunTFIDF and RunSVD to reduce sparsity, which is common in snATAC-seq data. We then computed a joint neighbour graph that represents both gene expression and chromatin accessibility measurements using the weighted nearest-neighbour methods in Seurat v4 using the function FindMultiModalNeighbors. Finally, we plotted the joint UMAP embedding with RunUMAP and annotated the 119,431 cells into 7 cell types and 14 subtypes (same as in snRNA-seq) using canonical marker genes. We converted the Seurat object into an ArchR project to call peaks on the 7 cell types and 14 subtypes.

Xenium processing and analysis pipeline

For segmentation, we analysed the Xenium data by both cell and nuclei. For cellular boundaries, we used the default settings from the output of Xenium Analyzer, which is based on a 15-μm outward expansion of DAPI-stained nuclei. For nuclear boundaries, we performed resegmentation using the resegment command from 10X Genomics Xenium Ranger (v2.0) with parameters --expansion-distance=0 and --resegment-nuclei=true. For quality control and filtering, we used Squidpy³³ (v1.6.1), a tool for analyzing spatial single-cell data, to create a UMAP embedding for each sample. For sample aggregation, after inspecting the individual sample UMAP embeddings and quality control metrics, such as the total and unique number of transcripts per cell or nucleus, we aggregated a total of 18 Xenium samples. For clustering and annotation, clustering was performed separately for scXenium and snXenium, resulting in 550,520 cells in the scXenium UMAP embedding and 523,314 nuclei in the snXenium UMAP embedding. Both had a fixed number of 366 genes. The clustered object was annotated into the seven canonical cell types, on which DEG analysis was performed. DEG analysis was performed using MAST for both scXenium and snXenium, and the results were compared with the snRNA-seq MAST results. scXenium and snXenium used thresholds |FC | > 1.1 and FDR < 0.01, and snRNA-seq used thresholds |FC | > 1.2 and FDR < 0.01. For visualization, SpatialData⁶⁶ (v0.1.2) was used to visualize the snXenium slides to highlight each cell-type-annotated nuclei and individual transcripts in their exact spatial coordinates. A single CON slide and a single PTSD slide were used to display differences in gene expression for a particular gene across the entire DLPFC, although the exact difference in gene expression was quantified using all CON samples versus all PTSD samples. This is shown as a barplot of log-normalized counts comparing the expression of a gene in CON versus PTSD across the seven cell types.

Differential abundance analysis

We calculated the differential abundance of cell types between conditions for each data modality using single-cell compositional data analysis (scCODA (v0.1.9))¹⁷. For RNA and Multiome, we ran the analysis across the 14 subtypes; for ATAC and Xenium, we ran the analysis across the 7 cell types. For each modality, we used the CompositionalAnalysis function with parameter reference_cell_type=‘automatic’, which automatically selects the cell type with the least amount of proportional change across conditions. We performed this analysis for PTSD versus CON, MDD versus CON, and PTSD versus MDD to compare across the three diagnostic cohorts. At FDR < 0.05, we found a significant decrease in OLGs for PTSD versus CON and PTSD versus MDD.

Differential gene expression analysis

snRNA-seq-based differential expression analysis was performed for each of the 14 subtypes using four tests: (1) MAST¹⁹ (v1.21.3) with covariate correction including covariates age death, sex, ancestry, PMI, RNA integrity number and library size; (2) Wilcoxon without covariate correction using the FindMarkers function from Seurat (v5.2.1)⁶⁷; (3) Nebula (v1.5.3)²⁰ for correcting sample-specific effects with covariate correction; and (4) DESeq2 (v1.46.0)²¹ with covariate correction. For MAST and Wilcoxon, we first log-transformed the CPM (counts per million)-normalized values of the raw count matrix. Only genes expressed in at least 5% of the cells in either condition, the robustly expressed genes, were included in the analysis. For Nebula, we only used genes expressed in at least 20% of the cells in either condition and removed samples that had less than 30 cells to reliably estimate the dispersion parameters. In the Nebula function, we used an offset term that combines n_counts and n_genes for each cell (0.8 × nCount + 0.2 × nFeature), the negative binomial gamma mixed model, kappa = 800 and cutoff_cell = 20. For DESeq2, we removed genes with either less than 20 reads per sample or a CPM value less than 0.5, and removed MDD samples with a very low number of read counts. We then utilized the default settings recommended for single-cell analysis such as test = ‘LRT’, minmu = 1 × 10^−6 and minReplicatesForReplace = Inf. We set the fold-change threshold for significant DEGs at |FC | > 1.2; we also constrained the significant DEGs to be FDR < 0.01 for MAST and Wilcoxon, FDR < 0.1 for Nebula, and FDR < 0.05 for DESeq2 pseudobulk.

For snRNA-seq, we first performed DEG analysis comparing PTSD versus CON and MDD versus CON using all samples from the diagnostic cohort. We identified PTSD snDEGs (n = 1,184) and MDD snDEGs (n = 1,918) by overlapping MAST and Wilcoxon DEGs that were concurrent in direction across the seven cell types. We then compared these snRNA-seq DEGs to the bulk RNA-seq DEGs from Girgenti et al.⁹ and Chatzinakos et al.²⁶ studies. These DEG sets were used in downstream analyses, including cross-disorder comparisons between PTSD and MDD, comparisons with CLGs in the ATAC analysis, and the identification of downstream targets of cell-type-specific transcription factors. We also investigated sex differences by first performing a baseline comparison of CON (male) versus CON (female), followed by PTSD (male) versus CON (male), PTSD (female) versus CON (female), MDD (male) versus CON (male), and MDD (female) versus CON (female). To examine PTSD samples comorbid with MDD, we performed PTSD (+MDD) versus CON and PTSD (−MDD) versus CON analyses and compared the results to PTSD versus CON in a three-way set relationship. For snXenium and scXenium, we performed DEG analysis using MAST and compared the results to snRNA-seq MAST results to validate the DEGs.

PTSD versus MDD cross-disorder comparison

To investigate the convergent and divergent properties between PTSD and MDD, we first overlapped the PTSD snDEGs and MDD snDEGs, which yielded 502 DEGs specific to PTSD and 1,236 DEGs specific to MDD. Next, we used the Rank–Rank Hypergeometric Overlap (v1.46.0)³⁰ method to identify concordant and discordant gene signatures between PTSD and MDD in a threshold-free manner. We identified 12 significantly discordant (down in PTSD and up in MDD) genes across cell types. We separated the PTSD cohort into those with co-existing MDD (n = 27 samples) and those without MDD (n = 8 samples) and performed DEG analysis against CON to identify several genes specific to PTSD. We also identified PTSD-specific and MDD-specific DEGs from the snXenium and scXenium MAST DEG analyses.

Gene Ontology enrichment analysis

We performed all Gene Ontology enrichment analyses using the Enrichr R package (v3.4)⁶⁸. We conducted a total of five Gene Ontology Enrichr analyses for: PTSD snDEGs, PTSD cell-type-specific DEGs, 502 PTSD-specific DEGs, PTSD Nebula DEGs and PTSD DESeq2 pseudobulk DEGs. We used the Enrichr function to query the Gene Ontology Molecular Function 2023 and Gene Ontology Biological Process 2023 databases. We highlighted the top FDR significant terms and removed duplicates for all Gene Ontology barplots.

CCC analysis

We applied the standard workflow of CellChat⁶⁹ (v1.5.0), utilizing single-cell gene expression of ligands and receptors to infer a CCC network. We used the normalized count matrix and cell-type annotations from the snRNA-seq dataset. Among all possible ‘sender’ (a cell expressing ligand) and ‘receiver’ (a cell expressing receptor) pairs, we constructed a network of communication strength based on the expression levels of ligand and receptor genes across cell types. High communication strength exists if the sender cell type has a high expression of ligand genes and the receiver cell type has a high expression of receptor genes, and vice versa. The number of inferred ligand–receptor pairs depends on the method for calculating the average gene expression per cell group. We utilized the robust method in CellChat for mean calculation called ‘truncatedMean’, with the cut-off ‘trim’ parameter set to the default value of 0.05. Two separate analyses were done for each diagnostic condition, namely, PTSD and CON. We then identified the differential CCCs across conditions. Most importantly, we normalized the 3D matrices (sender cell types × receiver cell types × ligand–receptor interactions) in each of the CellChat objects to have the same sum to account for batch effects of snRNA-seq. We then summed across the ligand–receptor interactions to reduce the 3D matrices into a 2D CCC network. A differential communication network was then constructed through the subtraction of these two networks between PTSD and CON groups. Nodes were coloured based on the aggregation of sending communication strength, or receiving communication volume, for output or input, respectively. We further highlighted specific cellular communication differences by using a subset of the 3D matrix (that is, MG cell type × receiver cell types × ligand–receptor interactions), effectively reducing it into a 2D object.

Neuronal-specific CCC analysis

We applied the standard workflow of NeuronChat (v1.1.0)³⁶, which utilizes single-cell gene expression of ligands and receptors to infer a neural-specific CCC network. We used the normalized count matrix and cell-type annotations from the snRNA-seq dataset. The number of inferred ligand–receptor pairs depends on the method for calculating the average gene expression per cell group. We utilized the default method in NeuronChat of mean calculation called ‘trimean’, with the cut-off ‘trim’ parameter set to 0.1. NeuronChat places greater emphasis on the second quartile, followed by the third and fourth quartiles, of gene expression. Two separate analyses were done for each diagnostic condition, namely, PTSD and CON, before moving on to comparative differential analyses. For the differential comparison analyses, we normalized the 3D matrices (sender cell types × receiver cell types × neuroligand–receptor interactions) in each of the CellChat objects to have the same sum to account for batch effects across snRNA-seq samples. We then summed across the neuroligand–receptor interactions to reduce the 3D matrices into 2D CCC networks. A differential communication network can then be constructed through the subtraction of these two networks between PTSD and CON groups. Sender (ligand) communication based on gene expression levels was aggregated to calculate the overall sending output for each cell type. The difference in sending cell signalling expression was calculated by subtracting between PTSD and CON groups. Specific pathways, or slices of the 3D matrices, were taken to explore which neurotransmitter–receptor interactions contributed to the difference between PTSD and CON.

Peak calling on snATAC-seq clusters

Pseudobulk replicates were created for each of the seven cell types with the addGroupCoverages function in ArchR to resolve inherent sparsity in the snATAC-seq data. MACS2 (v2.2.9.1)⁷⁰ was then used to call peaks on each pseudobulk replicate. ArchR uses an iterative overlap peak merging procedure with the addReproduciblePeakSet function to create a final peak set of fixed-width peaks using the following procedure: (1) peaks were first ranked by their significance, (2) the most significant peak was retained and any peak that directly overlapped with the most significant peak was removed from further analysis, and (3) of the remaining peaks, this process was repeated until no more peaks existed. ArchR analyses all the pseudobulk replicates from a single-cell type together, performing the first iteration of iterative overlap removal, then checks to see the reproducibility of each peak across pseudobulk replicates and only keeps peaks that pass a threshold indicated by the reproducibility parameter, resulting in a single merged peak set for each cell type. This procedure was repeated to merge the cell-type-specific peak sets by renormalizing the peak significance across the different cell types and performing the iterative overlap removal, resulting in a single merged peak set of fixed-width (501 bp) peaks, denoted the union peaks. We utilized both cell-type-specific peaks and the union peaks depending on the analysis that we performed. The addPeakMatrix function adds the PeakMatrix using the union peaks to the data object, which was used for downstream analyses such as adding peak-to-gene links and performing transcription factor motif analysis. To assess the quality of the peak calling procedure, we overlapped our cell-type-specific snATAC-seq peaks with the bulk adult bCREs, which were derived from adult brain tissues to study the regulatory landscape of the human brain⁷¹. From each of the 96 DNase-seq experiments, V4 cCREs (candidate cis-regulatory elements) with a Z-score > 1.64 were selected. V4 cCREs were obtained directly from ENCODE (https://screen.encodeproject.org)⁴⁰. Subsequently, cCREs with limited experimental support were removed and only those elements that were supported by at least five experiments were retained, which resulted in a collection of 253,638 adult bCREs. We used bedtools intersect (v2.30.0) with the minimum overlap fraction (-f), ranging from 0.2 to 1.0 to calculate the percentage of overlap.

snRNA-seq and snATAC-seq integration

We integrated our snATAC-seq with a subset of our snRNA-seq data (six high-quality samples that displayed high heterogeneity of cell types) using the addGeneIntegrationMatrix function. We used a constrained integration approach, which constrains the cells in the snRNA-seq data to the best-matched cells in the snATAC-seq data for the seven cell types. Once we had the GeneIntegrationMatrix, we identified peak-to-gene links using the addPeak2GeneLinks function, which leverages integrated snRNA-seq data to look for correlations between peak accessibility and gene expression. We set maxDist = 2 MB to consider long-range interactions. We used the plotPeak2GeneHeatmap function with k-means clusters = 20 to plot the side-by-side heatmaps of linked ATAC and gene regions.

Characterization of cell-type-specific CREs

Once we had the peak-to-gene links, we used thresholds correlation > 0.4 and FDR < 0.05 to identify peaks likely to regulate gene expression. The 395,932 cell-type-specific CREs, which are the unique peaks from this subset, are the backbone of all ATAC-related analyses. We identified CLGs and overlapped them with the cell-type-specific DEGs to identify CLG–DEGs for each cell type. These CREs were also used in the GWAS analyses to assess LDSC trait heritability and fine-map causal variants at PTSD risk loci in a cell-type-specific manner.

Characterization of condition-specific CREs

To investigate epigenomic differences between PTSD and MDD, we curated a set of condition-specific CREs in our ATAC data. First, we called peaks on each cell-type condition group (for example, EXN CON, EXN MDD and EXN PTSD) in the same manner as was done on the seven cell types. We identified unique peaks for each condition using bedtools intersect. After adding peak-to-gene links and identifying CREs that were likely to regulate gene expression, we identified a total of 141,939 condition-specific CREs across all cell-type condition groups. We then identified PTSD CLG–DEGs by overlapping the PTSD-specific CRE-linked genes and PTSD DEGs and MDD CLG–DEGs by overlapping the MDD-specific CRE-linked-genes and MDD DEGs for each cell type. Globally, we found 383 PTSD-specific CLG–DEGs and 1,063 MDD-specific CLG–DEGs that were specific to each condition.

snMultiome ATAC analysis

As we achieved subtype resolution in the neuronal clusters with the snMultiome data, we were able to call peaks on the seven cell types and 14 subtypes. We compared these peaks to the bCREs to assess peak calling quality and also compared them with the cell-type-specific ATAC peaks. The peak overlap ratio between ATAC and Multiome was calculated using bedtools intersect with a minimum overlap of one nucleotide. The ratio was then determined by dividing the number of intersecting peaks by the total number of ATAC peaks. With subtype resolution of neurons, we were able to identify the specific EXN or IN subtype that was most enriched in a particular trait using LDSC.

Transcription factor motif enrichment analysis and footprinting

Transcription factor motif enrichments were calculated to predict which regulatory factors were most active in a given cell type. The addMotifAnnotations function in ArchR indicates motif presence (P < 5 × 10^−5) in the union peaks using a binary matrix, and the peakAnnoEnrichment function tests either cell-type-specific marker peaks or differential peaks for enrichment of various motifs. Besides the previous transcription factor motif enrichments, the R package chromVAR predicts enrichment of transcription factor activity on a per-cell basis while controlling for known technical biases. The addDeviationsMatrix function, which adds the MotifMatrix to the data object, was used to compute deviation z-scores for each motif by comparing the number of fragments that map to peaks containing the motif to the expected number of fragments in a background peak set that accounts for confounding factors such as GC content bias, PCR amplification and Tn5 tagmentation. Then, we used the getPositions, addGroupCoverages and getFootprints functions to perform transcription factor footprinting analysis with pseudobulk aggregates of cells in the same cell type. In the plotFootprints function, ArchR performs normalization by subtracting the Tn5 bias from the footprinting signal. We plotted the footprinting results for the transcription factors EGR2 and TFAP4 comparing PTSD and CON in EXNs and INs, respectively.

Transcription factor gene-regulatory network construction

Transcription factor–CRE–gene links were created for each cell type with DEG (up or down) and PTSD GWAS gene (true or false) status indicated for each gene. We created GRNs (gene regulatory networks) using EXN transcription factor–CRE–gene links and IN transcription factor–CRE–gene links using three highly enriched transcription factors in each cell type. We used the R package igraph (v1.2.6) to plot the GRN, showing transcription factor-to-target gene links and overlaying additional information such as whether the target gene is a DEG or PTSD GWAS gene. For the GRNs, we used correlation cut-offs of 0.7 and 0.6 for EXNs and INs, respectively, for visualization purposes. For the EXN GRN, we indicated all genes that overlapped with the iPS cell DEGs for EGR2 and SMARCC1 in pink. Across all transcription factor–CRE–gene linkage analyses, we used the FDR value of the transcription factor–CRE linkage as a way to test multiple comparisons and indicate statistical significance.

Estimating GWAS enrichment using ATAC and Multiome peaks

To estimate heritability of various complex traits, we used LDSC (v1.0.1)⁷². GWAS summary statistics for PTSD-related traits (tota lPCL, case–control, re-experiencing, avoidance and hyperarousal), other psychiatric traits (schizophrenia, MVP depression, Alzheimer’s disease and alcoholism) and non-psychiatric traits (irritable bowel disease, diabetes, height and education) were correctly formatted with munge_sumstats.py and lifted to hg38 coordinates using UCSC liftover. Cell-type-specific peaks were formatted for LDSC using make_annotation.py, and linkage disequilibruim scores were computed for each set using ldsc.py. Benjamini–Hochberg multiple-testing correction was applied to the enrichment P values. We compared the results using single-cell peaks from snATAC-seq or snMultiome with the bulk bCREs, and saw enhanced enrichment using single-cell data.

GWAS fine-mapping

For PTSD GWAS loci screening, we first identified risk loci from two previous PTSD GWAS^5,7. We selected a neighbourhood of around 1–2 Mb that included the MVP-lead SNP. After fine-mapping these loci using the total PCL and re-experiencing summary statistics, we identified eight that fell within cell-type-specific EXN and IN CREs. For peak enrichment scores as prior weights, we performed genetic fine-mapping with the sum of single-effect (SuSiE, v0.12.35) regression model⁷³ to infer risk variants with epigenetic information from our snATAC peaks, such that the prior can improve the accuracy of detecting causal variants⁷⁴. For SNPs within ATAC peaks, we used the peak enrichment scores (−log₁₀P), which indicate the strength of each peak, as the prior weights for SuSiE to prioritize SNPs with strong ATAC signals. The enrichment scores in the union peakset ranged from 1 to 507 with a mean of 19.7 and standard deviation of 28.3. To reduce the effects of highly enriched peaks, we set the enrichment score values above the 95th percentile to the value of the 95th percentile (67.49) and rescaled the scores to 1–100. For SNPs outside ATAC peaks, we set the prior weights to 0.1. We used CREs (peak-to-gene correlation > 0.4 and FDR < 0.05) to perform fine-mapping in SuSiE. For fine-mapping with SuSie, we performed SuSiE regression with the susie_rss function using the PTSD GWAS summary statistics (total PCL and re-experiencing)⁵. SuSiE calculates the PIP, the probability for a given SNP being causally associated with the trait of interest, for each variant. The linkage disequilibrium information was derived from the 503 European samples of the 1000 Genomes Project (1KG). For input files, we used the 1KG Phase 3 dataset⁷⁴ as the reference panel for our analysis. This dataset can be accessed via the Resources section on the PLINK 2.0 website (Resources — PLINK 2.0 (cog-genomics.org)). We selected the European samples based on the super-population information provided by 1KG and excluded all duplicate and ambiguous SNPs. We applied quality control to the 1KG data with European ancestry using PLINK⁷⁵. To be more specific, we only included SNPs with a 99% genotyping rate and the minor allele frequency no less than 0.005. We also excluded samples with more than 5% missing genotypes and markers that failed to pass the Hardy–Weinberg test. The final reference panel included 503 non-overlapping European samples, genotyped at 8,190,311 SNPs. To obtain the coordinates for risk SNPs, we turned to the UCSC Genome Browser⁷⁶ and focused on the human genome version GRCh38/hg38. For fine-mapping figures, after running SuSie, we integrated all the findings including the ATAC signal track, the CRE–gene links that link to the TSS of the risk gene, the PIP output values from SuSie, correlation of each variant to the lead SNP, and the gene track into the final fine-mapping figure using the predetermined regions. We plotted all variants with a PIP > 0.0001, highlighted CRE–gene loops for variants with PIP > 0.05, and outlined each variant in black if it fell within a CRE. The top two credible SNPs and the MVP SNP were labelled. The lead SNP was coloured red, and the MVP SNP was coloured purple. The CRE–gene loops were colored by their CRE–gene correlation values and by cell type (IN in green, and EXN in red). For HiC TAD, for ELFN1, KCNIP4, EGR3 and LRFN5, contact domains were found from the NeuN⁺ Hi-C file⁷⁷ using the Juicer (v1.6)⁷⁸ arrowhead command with parameters r = 10,000 (the resolution of the sequencing depth of the Hi-C file) and m = 1,000 (the size of the sliding window along the diagonal in which contact domains will be found). The TADs highlighted the regions of chromosomal looping in the fine-mapping figures.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.