Cellular Diversity in Human Subgenual Anterior Cingulate and Dorsolateral Prefrontal Cortex by Single-Nucleus RNA-Sequencing

Contact for reagent and resource sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by S.M. (marencos{at}mail.nih.gov).

Experimental design and statistical analysis

snRNA-seq data were generated from the DLPFC and sgACC of 8 male subjects (age = 44.8 ± 11.2 SD; sex = male; race = 5 White and 3 Black; postmortem interval = 24.3 ± 7 h; pH = 6.5 ± 0.17), each without a history of psychiatric (including substance use disorders) or neurologic disorders during their lifetime (for detailed demographics, see Extended Data Fig. 1-1). Technical replicates for all 16 unique specimens at the droplet-based nuclei capture step resulted in 32 total samples. Quality control steps (detailed in snRNA-seq data preprocessing) eliminated seven samples with low total cDNA yield resulting in 25 samples with 105,047 nuclei and 28,195 genes. Pearson correlations were performed to compare technical replicates.

Unsupervised clustering algorithms in the Seurat version 4.1.0 (Hao et al., 2021) pipeline were used to assign nuclei into general cell classes and to identify and eliminate doublets, clusters with high mitochondrial contamination, and artifact nuclei (for details, see Unsupervised clustering, quality control, and clustering) yielding a final dataset of 63,901 nuclei. Additional subclustering was performed on Ex and In cells using Seurat version 4.1.0 tools. The default nonparametric Wilcoxon rank sum test with Bonferroni correction was applied to identify cluster markers. Pearson correlation was used to compare our cell classification to published snRNA-seq datasets (for details, see Comparison with published human snRNA-seq datasets).

To define consensus cell class identities, our data were integrated with the Allen Cell Types Database Human Multiple Cortical Areas SMART-seq dataset via reciprocal PCA in Seurat version 4.1.0, and a nearest centroid classifier was applied to map each cell onto the cluster labels in the Allen dataset. Regional subpopulation proportional representation was assessed using the GLIMMIX (generalized linear mixed model) procedure with logit function and the option of β distribution with default covariance structure in SAS version 9.4 to test for effects of subpopulations, brain region, and their interaction (see Regional differences in subclusters). Differential expression between regions and subclusters was tested using a pseudo-bulk approach using the edgeR likelihood ratio test (for details, see Regional differences in subclusters).

Linkage disequilibrium score regression (LDSC) was implemented to examine cell type enrichment for GWAS signals from studies of ASD (Grove et al., 2019), attention-deficit/hyperactivity disorder (ADHD) (Demontis et al., 2019), BD (Stahl et al., 2019), MDD (Wray et al., 2018), and SCZ (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014) (for details, see Cell type-specific genes and psychiatric disorder enrichment).

Postmortem human specimen collection and dissection

Postmortem brains were collected at the Human Brain Collection Core (HBCC), National Institute of Mental Health Intramural Research Program, with permission from the legal next-of-kin according to the National Institutes of Health Institutional Review Board and ethical guidelines under protocol 17 M-N073. Cases were obtained from the Offices of the Chief Medical Examiner of the District of Columbia, Northern Virginia, and Central Virginia. Clinical, neuropathological, and toxicological characterization was performed as previously described (Deep-Soboslay et al., 2005; Lipska et al., 2006; Kunii et al., 2015).

Samples were dissected from frozen coronal slabs using a dental drill. For the DLPFC (BA9/46) dissections, gray matter tissue from the midpoint of the middle frontal gyrus was targeted. Subcortical white matter was carefully excluded. For the sgACC (BA24/25), the portion of the cingulum immediately ventral to the genu of the corpus callosum was dissected.

snRNA-seq nuclei extraction and isolation

The nuclei isolation protocol was adapted from Krishnaswami et al. (2016) with additional washes to reduce mitochondrial contamination. To maximize our ability to capture the entire cellular diversity in an anatomically defined brain region during dissection, we obtained two separate 50-80 mg aliquots from each frozen dissected block and combined them after nuclei isolation.

Frozen samples were physically dissociated in a homogenizing buffer using a Dounce homogenizer (DWK Life Sciences), five strokes with a narrow pestle followed by 15 strokes with a wide pestle. The resulting nuclear suspensions were filtered through a 35 μm mesh cell strainer. Nuclear suspensions were washed 3 times, pelleting the nuclei with centrifugation at 1200 × g for 9 min before resuspending the pellets in 0.5 ml of wash buffer. Nuclear suspensions from the same frozen block (described above) were combined at this step and thoroughly mixed.

For nuclei counting, a 10 μl aliquot of each nuclear suspension was stained 1:1 with 0.4% Trypan Blue. Nuclei were counted using a hemocytometer under a bright-field microscope. Concentrations ranged from 1 million to 1.5 million nuclei per milliliter. A final aliquot of 50,000 nuclei in 1 ml of storage buffer was taken for the 10× Genomics nuclei capture.

10× genomics nuclei capture and library construction

We followed the Chromium Single Cell 3′ Reagent Kits User Guide version 2 Chemistry (CG00052) and used Chromium Single Cell 3′ Library & Gel Bead Kit version 2 (PN-120237) to produce Single Cell 3′ libraries ready for Illumina sequencing (https://support.10xgenomics.com/). Droplet-based nuclei-capture and library preparation experiments were conducted using materials and support from the Single Cell Analysis Facility, National Cancer Institute.

First, we created technical duplicates from the nuclei suspensions and loaded them onto the Chromium Controller microfluidic device. We loaded volumes provided in the Cell Suspension Volume Calculator Table of the 10× user guide to target a recovery rate of 6000 nuclei per sample. The microfluidic device partitioned single nuclei into bar code-containing Gel Bead-In-Emulsions (GEMs). Samples were incubated, and poly-adenylated mRNAs were reverse-transcribed into cDNAs containing a unique molecular identifier (UMI) and a cell bar code unique to each GEM. Afterward, GEMs were broken, and the barcoded cDNAs were amplified by PCR. Next, a sample index, R2 (read 2 primer sequence), and standard P5 and P7 primers used in Illumina bridge amplification were added during library construction via End Repair, A-tailing, Adaptor Ligation, and PCR. The 16 bp 10× bar code and 10 bp UMI were encoded at the start of Read 1 while the 8 bp sample index information was incorporated into the i7 index read, and cDNA insert sequenced as Read 2.

Sequencing, bar code processing, and alignment

Sequencing was performed by the Center for Cancer Research Sequencing Facility at the Frederick National Laboratory for Cancer Research. Thirty-two 10× Genomics Chromium Single Cell 3′ version 2 paired-end libraries were multiplexed and sequenced on one run on the HiSeq 4000. Base-calling was performed with Illumina instrument runtime analysis software, RTA version 2.4.11. The output raw base call (BCL) files were demultiplexed using Bcl2fastq version 2.17, allowing one mismatch in barcodes, to generate FASTQ files. Sequence alignment and bar code processing were performed with 10× Genomics's processing software, Cell Ranger version 2.1.1 (Zheng et al., 2017). A customized human reference genome, GRCh38-1.2.0_premRNA (described below), was used for counting reads and the “expect-cells” parameter for estimating number of nuclei in each sample above a background expression level was set to 6000, consistent with the number of nuclei targeted during library preparation.

Because snRNA-seq captures precursor mRNA (pre-mRNA) transcripts that have not been alternatively spliced, these transcripts have retained intronic features. Recent RNA-seq and comparative studies of single nucleus and single-cell RNA-seq datasets have shown that including reads mapping to intronic regions as gene counts accurately estimates true gene expression (Gaidatzis et al., 2015; Lake et al., 2017; Bakken et al., 2018). Additionally, this method has an added benefit for downstream clustering analysis by increasing read depth for inherently low coverage snRNA-seq datasets. Therefore, a customized, pre-mRNA reference genome was created following recommendations from 10× Genomics. “Gene_biotype” tags in the original GTF file were modified such that Cell Ranger's count function included intronic features as gene counts.

snRNA-seq data preprocessing

Data processing, analysis, and handling were performed using the R package Seurat version 4.1.0 (Hao et al., 2021). Filtered gene count matrices were read into R using Seurat's Read10X function. After initial review of experimental and sequencing quality control (QC) characteristics, 7 of 32 total samples were identified as “dropout” samples and removed from the dataset. These samples displayed low total cDNA yield calculated from BioAnalyzer electropherograms. They contained high sequence duplication rate (>45%) calculated with the sequencing QC tool FastQC (Andrews, 2010) and high sequencing saturation rate (>85%) from Cellranger QC outputs. Together, these QC results suggested inadequate RNA capture and subsequent PCR overamplification of these “dropout” samples Twenty-five samples (14 DLPFC, 11 sgACC) remained for downstream analysis and were assessed for further quality control. Detailed 10× Genomics Library capture, sequencing, and alignment data for each sample of brain tissue is listed in Extended Data Fig. 1-2.

Filtering of noisy genes and outlier nuclei was performed; 5475 genes detected in fewer than three nuclei were filtered from the data. Nuclei with <200 genes do not contain enough transcriptional information for effective clustering and subsequent biological interpretation. This threshold is like those used in published snRNA-seq studies (Hu et al., 2017). No nuclei fell below this threshold in our data. Nuclei with large library sizes (>5000 genes detected and 15,000 UMIs) are likely to be “doublets” or “multiplets,” which are expected technical artifacts of the microfluidic droplet technology (Habib et al., 2017; Hu et al., 2017; Lake et al., 2018). We excluded 227 nuclei whose values exceeded these thresholds (0.2% of total nuclei). Next, the fraction of reads mapping to mitochondrially encoded genes was calculated for each barcoded nucleus. Mitochondrially encoded mRNA transcripts represent false transcriptional signals and thus should not be present in snRNA-seq data. Nuclei with >20% of reads mapping to mitochondrially encoded genes were filtered out (25,273 nuclei). Subsequently, mitochondrially encoded genes were removed from the count matrix. After these preprocessing QC steps, 105,047 nuclei and 28,195 genes were retained for analysis.

Correlation between technical replicates of samples

To assess similarity between technical replicates, we calculated the Pearson correlation coefficient between average gene expression of replicates using the cor function in R (Extended Data Fig. 1-3). Correlations of average expression across all genes and then across only highly variable genes were assessed. Correlation of gene expression (Pearson r > 0.998) and correlation of number of nuclei obtained (Pearson r = 0.992) across the technical replicates was high (Extended Data Fig. 1-3). Given the high correlation coefficient values between technical replicates and the added clustering power that is achieved with larger sample size (number of nuclei), technical replicates were combined for subsequent downstream analyses. Technical replicates could be treated as additional samples of the same subject because they were separate aliquots from nuclei suspensions, representing different nuclei from the same subject.

Unsupervised clustering, quality control, and clustering

Normalization of gene expression to account for library size and expression level was performed with the Seurat function NormalizeData using default settings. Next, highly variable genes (HVGs) in our dataset were identified using the Seurat function FindVariableGenes for unsupervised clustering. Briefly, genes were first binned into 20 bins based on their mean expressions. Z scores were calculated within each bin. This method of binning controlled for the correlation between gene expression and variability. After binning, genes with SD > 0.65 and mean log normalized expression between 0.015 and 4 were selected as HVGs. Next, the Seurat function ScaleData was used to compute z scores for each gene in the list of HVGs. We calculated z scores by fitting gene expression to a negative binomial model because it more accurately resembles gene expression distributions. In this regression, the total number of UMIs detected and fraction of mitochondrial genome encoded gene counts per nucleus were regressed out to further eliminate library size effects on the downstream clustering analysis. Next, for linear dimensional reduction, a principal component analysis (PCA) was performed on the z scores of the HVGs using the Seurat function RunPCA. The top 50 principal components (PCs) were selected after assessment of PCs with an elbow plot and were used as input to Seurat's unsupervised clustering algorithm function FindClusters. In this algorithm, Euclidean distances between nuclei were calculated in the PCA space (top 50 PCs) and a K-nearest neighbor (KNN) graph-based approach was used to group nuclei into clusters. The resolution parameter, which represents the granularity of the clustering analysis, was set at 3, and all other parameters were set at default settings. For two-dimensional t-distributed stochastic neighbor embedding (tSNE) visualization, the RunTSNE function from Seurat was performed at default settings on the same principal components that were used for clustering.

After clustering, a variety of methods were used to identify and disqualify clusters. The distribution of number of genes, number of UMIs, fraction of mitochondrial gene counts, and fraction of ribosomal protein gene counts were visually assessed for each unbiased cluster using violin plots. To assess similarity of these clusters to each other, hierarchical clustering analysis on the average expression for each cluster was performed using R function hclust, a pairwise correlation heatmap between all clusters was generated with R function cor to view similar mean expression profiles, and distances of tSNE coordinates were visually assessed. Doublet or mixed clusters, identified by high numbers of genes and UMIs detected in addition to the combined expression of multiple canonical cell class marker genes, and clusters with very low numbers of detected genes, low UMI counts, and lack of distinguishing expression profiles, were discarded. The remaining subset of nuclei were reclustered using new HVGs calculated from the subset. Cluster removal and reclustering were repeated until all clusters that formed could be identified by their specific expression of canonical cell-type marker genes from existing literature and could not be removed by quality control as described above. A total of 87,729 of 105,047 nuclei remained after cluster curation.

Finally, to remove artifact nuclei, which include collections of ambient RNA encapsulated during the droplet experiment or doublets/multiplets, we curated each cluster by calculating cell-type scores for each nucleus (barnyard curation). Cell-type scores were calculated as a nucleus's average expression of the top five specific genes for each broad cell class, listed in Extended Data Figure 1-4. Top five specific genes were determined as the five highest-ranked genes when significant differentially expressed genes (DEGs) for each broad cell cluster were ranked by specificity (percent of in-group nuclei expressing the gene/percent of out-group nuclei expressing the gene). An additional criterion for top five specific genes for broad cell classes was that the genes were required to overlap with marker genes from at least one other single-cell or nucleus RNA-seq dataset (Darmanis et al., 2015; Lake et al., 2016; Habib et al., 2017; Hu et al., 2017; Sathyamurthy et al., 2018). By assessing each nucleus's score for each of the broad cell cluster, nuclei with “promiscuous” cellular identity (i.e., having scores greater in cell classes other than the cell class to which it belongs) were filtered out from further analysis. A total of 63,901 of 87,279 nuclei remained after barnyard curation.

Following QC, 63,901 nuclear expression profiles remained from 25 samples, averaging 93 million reads per sample and 27,061 reads per nucleus, slightly above the range (19,000-26,000 reads per nucleus) typically required for sequencing saturation with 10× Genomics version 2 chemistry in human brain snRNA-seq (Habib et al., 2017). Using 10× Genomics' software Cell Ranger for read alignment and counting, 1641 median UMIs and 1100 median genes per nucleus were detected, which is consistent with other snRNA-seq studies in postmortem human brain samples (Habib et al., 2017; Lake et al., 2018). A modified human reference genome was used to incorporate intronic reads as gene counts (Gaidatzis et al., 2015; Lake et al., 2017; Bakken et al., 2018). Ninety percent of reads mapped confidently to the customized genome.

Determination of cluster marker genes

The FindMarkers and FindAllMarkers functions in Seurat were applied with default settings to identify marker genes for each cluster (i.e., genes enriched in a cluster compared with all other clusters). The default test used was a nonparametric Wilcoxon rank sum test. Only genes whose expression value were > 10% of the nuclei in the in-group, whose Bonferroni corrected p values < 0.05, and whose log fold-change > 0.25 were saved as marker genes. This approach was used to identify marker genes for broad cell clusters (Extended Data Fig. 1-5), and excitatory and inhibitory neuron subclusters (Extended Data Fig. 2-1).

Subclustering inhibitory and excitatory neurons

Transcriptional signals that differentiate neurons from glial cells may not be the same signals that differentiate subtypes of neurons, excitatory and inhibitory, from each other (Tasic et al., 2018). The signals that differentiate excitatory from inhibitory neurons, in turn, may also not be the same signals that differentiate canonical inhibitory neuron subtypes from each other. To effectively resolve the transcriptional diversity of inhibitory neuron subtypes without the influence of noninhibitory neurons, we performed subclustering on the restricted set of nuclei identified as inhibitory neurons. We implemented the same data-driven QC approaches described above for the broad cell class clustering analysis with an added step of merging nonunique clusters. From an initial 24 inhibitory neuron subclusters, we reassigned nuclei of one subcluster that did not have a uniquely identifiable expression signal (<10 unique DE genes) and one subcluster that mostly consisted of nuclei from one individual. This process is called “merging” and is performed by reassigning the identity of each nucleus individually to the subcluster with the highest average shared nearest neighbors. This process reassigned nuclei to subclusters that had very similar transcriptional profiles to its original cluster. Similarity of subclusters could be assessed by correlation heatmaps and relative position to one another in a hierarchical clustering tree. This resulted in 22 final inhibitory neuron subclusters. DE analysis was performed as described above for each subcluster against all other inhibitory neurons to identify potential marker genes. Subclustering was also performed for the excitatory neurons. From the initial 29 subclusters that formed, nine subclusters were identified as nonunique cell types and nuclei of these subclusters were reassigned to subclusters with highest average shared nearest neighbors. Merging these subclusters resulted in 20 final excitatory neuron subclusters represented by all 8 individuals. Statistics for the final list of subclusters are shown in Extended Data Figure 2-1.

Comparison with published human snRNA-Seq datasets

For comparison with published human brain snRNA-Seq datasets, we read in count matrices (Lake et al., 2016; Habib et al., 2017; Lake et al., 2017; Tasic et al., 2018; D. Wang et al., 2018) into R as Seurat objects. Average gene expression was computed for each cluster from each dataset. A Pearson correlation coefficient was calculated using R function, cor, on mean expression for all pairwise comparisons of clusters. We limited our genes for correlation comparisons to 3310 HVGs (disregarding 78 ribosomal protein genes) to reduce noise from gene expression signals that are less informative among brain cells and to avoid dataset-specific technical artifacts. All heatmaps were plotted with either the ggplot2 R package or the pheatmap R package with colors palettes generated from the RColorBrewer R package. Correlation heat maps are shown in Extended Data Figure 2-4.

Dataset integration in Seurat and nearest centroid classification

snRNA-seq count data from multiple cortical areas were downloaded from the Allen Cell Types Database (https://portal.brain-map.org/atlases-and-data/rnaseq/human-multiple-cortical-areas-smart-seq) and imported into Seurat (version 4) along with HBCC data. Nuclei with a missing cluster label were removed. Inhibitory neurons, excitatory neurons, and non-neuronal cells were analyzed separately. For each major subclass, cells were log-normalized and the top 2000 most variable genes across both datasets were selected using the SelectIntegrationFeatures function in Seurat. PCA was run separately on both HBCC and Allen Institute data for the selected genes. Integration anchors were obtained via the FindIntegrationAnchors function using reciprocal PCA on the top 30 principal components and a k.anchor default value of 5. Datasets were integrated via the IntegrateData function with a k.weight default value of 100.

For classification of HBCC data based on Allen Institute labels, a nearest centroid classifier was implemented via the lolR library and was trained on the Allen cluster labels with the integrated gene expression values from 2000 genes as described above. This classifier was then applied to predict the Allen cluster label of each cell in the HBCC data. A contingency figure of the original HBCC cluster ID and predicted Allen label was created, and the fraction of HBCC cells mapping to each Allen label was visualized in a heatmap (see Fig. 2E,F).

Regional differences in subclusters

To establish whether particular neuronal subpopulations were unequally represented in the DLPFC or the sgACC, we calculated the proportion of neuronal subpopulations using all neurons as the denominator, then used the GLIMMIX (generalized linear mixed model) procedure with logit function and the option of β distribution with default covariance structure in SAS version 9.4 to test for effects of subpopulations, brain region, and their interaction. Subject (id) was used as random effect, assuming subjects participating in the study are a random sample of all possible subjects. This procedure was repeated on absolute counts of each neuronal subpopulation (omitting the β distribution).

Differential expression analysis between subclusters and brain regions was conducted using a conservative pseudo-bulk approach via edgeR likelihood ratio test (Robinson et al., 2010), while adjusting for donor status using the Libra R package (Squair et al., 2021). DEGs were defined as those that passed false discovery rate (FDR)-corrected p < 0.05 and had an absolute logFC > 0.1. To explore the relevance of DEGs, these genes were used as input into gene enrichment analysis with the Gene Ontology (GO) (Gene Ontology Consortium, 2001), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), and DisGeNET (Piñero et al., 2015) as implemented in the ClusterProfiler (Yu et al., 2012) and disgenet2r R packages (Piñero et al., 2020).

Cell type-specific genes and psychiatric disorder enrichment

To test the association between cell type-specific genes and the heritability for psychiatric disorders, we used LDSC, as implemented previously (Finucane et al., 2015; Skene et al., 2018; Bryois et al., 2020). In brief, single-nucleus RNA expression data were scaled to a total of 1 million UMIs for each cell type, separately for each cortical region. A specificity metric was calculated as a ratio of the expression of each gene in a given cell type cluster over the total expression of that gene across all cell types (Skene et al., 2018; Bryois et al., 2020). The top 10% most specific genes for each cell type were used as input for LD score regression. Genome Wide Association Study (GWAS) data for the following psychiatric disorders were used: ASD (Grove et al., 2019), ADHD (Demontis et al., 2019), BD (Stahl et al., 2019), MDD (Wray et al., 2018), and SCZ (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014), accessed through the Psychiatric Genomics Consortium data sharing portal (https://www.med.unc.edu/pgc/). Gene coordinates were extended by 100 kb upstream and 100 kb downstream for each as done previously (Finucane et al., 2015; Bryois et al., 2020). p values were FDR-corrected for 7 broad cell types, and 47 subclusters (Benjamini and Hochberg, 1995).

The same top 10% of cell specific genes were used to test for enrichment with DEGs from a recently published large RNA seq study of ASD, SCZ, and BD samples (Gandal et al., 2018b). Hypergeometric tests were used to assess the enrichment or overlap between gene sets, and the p values were corrected as described above.

RNA ISH (RNAscope)

For RNAscope validation of cell types, 10-μm-thin cryosections were taken from the DLPFC and sgACC dissected from postmortem brain samples (Extended Data Fig. 2-7). The sections underwent post-fixation in 10% NBF at 4°C for 1 h. Tissue was dehydrated and stained following the Advanced Cell Diagnostics protocol for RNAscope 2.5 Duplex Chromogenic Assay (Anderson et al., 2016) with the following modifications: Protease IV treatment reaction time was reduced to 15 min; amplification reaction time was optimized to obtain the best resolution of puncta for the different target probes labeled with HRP-based green; counterstaining was done with 25% Hematoxylin 560 MX (Leica 042817) diluted in dH₂O, for 10 s; alternatively, Nissl staining was done with 1% thionine for 30 s, subsequent wash in dH₂O for 30 s twice, and differentiated in acidified ethanol for 5-15 s. Whole slide image acquisition was performed on a NanoZoomer S60 Digital Slide Scanner (Hammamatsu).

Code/software

Raw FASTQ files, final Seurat R Object containing nuclei expression data, and metadata are publicly available for download on the National Institute of Mental Health Data Archive (https://nda.nih.gov/experimentView.html?experimentId=1852&collectionId=3151) and Synapse (https://doi.org/10.7303/syn48374287.1). R code for this project is available for download on our GitHub repository: https://github.com/NIMH-HBCC. A shareable interactive user interface using iSEE (Rue-Albrecht et al., 2018) and Binder 2.0 (https://elifesciences.org/labs/8653a61d/introducing-binder-2-0-share-your-interactive-research-environment) is also accessible from the GitHub HBCC repository. Datasets and gene sets from referenced literature were obtained online from supplementary data.