Bulk and Single-Nucleus Transcriptomics Highlight Intra-Telencephalic and Somatostatin Neurons in Alzheimer's Disease.

Cortical neuron loss is a pathological hallmark of late-onset Alzheimer's disease (AD). However, it remains unclear which neuronal subtypes beyond broad excitatory and inhibitory classes are most vulnerable. Here, we analyzed cell subtype proportion differences in AD compared to non-AD controls using 1037 post-mortem brain samples from six neocortical regions. We identified the strongest associations of AD with fewer somatostatin (SST) inhibitory neurons (β = -0.48, p bonf = 8.98 × 10-9) and intra-telencephalic (IT) excitatory neurons (β = -0.45, p bonf = 4.32 × 10-7). Replication in three AD case-control single-nucleus RNAseq datasets most strongly supported the bulk tissue association of fewer SST neurons in AD. In depth analyses of cell type proportions with specific AD-related neuropathological and cognitive phenotypes revealed fewer SST neurons with greater brain-wide post-mortem tau and beta amyloid, as well as a faster rate of antemortem cognitive decline. In contrast, greater IT neuron proportions were associated with a slower rate of cognitive decline as well as greater residual cognition-a measure of cognitive resilience-but not canonical AD neuropathology. Our findings implicate somatostatin inhibitory and intra-telencephalic excitatory neuron subclasses in the pathogenesis of AD and in cognitive resilience to AD pathology, respectively.

Introduction

Late-onset Alzheimer’s disease (AD) is a neurodegenerative disease characterized by the gradual accumulation of specific neuropathologies, including beta amyloid and hyperphosphorylated tau proteins, followed by widespread brain atrophy and cognitive decline (Rombouts et al., 2000; Lerch et al., 2005, 2008). While these pathological hallmarks of AD are well established, a lack of clarity over which specific brain cell types are lost over the course of neurodegeneration and cognitive decline remains.

Recent advances in single-cell and cell type-specific gene expression profiling has revolutionized our knowledge of cell-type specific changes in neuropsychiatric disease (Mancarci et al., 2017; Gandal et al., 2018; Toker et al., 2018; Wang D. et al., 2018; Jew et al., 2020). By combining these datasets with bulk tissue gene expression data sampled from large numbers of well-characterized individuals, cellular deconvolution analyses have revealed important differences in AD, including fewer neurons and more astrocytes (Li et al., 2018; Patrick et al., 2020; Park et al., 2021). However, the majority of these analyses have focused on cellular differences at the level of broad cell types and comparatively less focus has been placed on resolving cellular differences among finer-resolution cell types such as subtypes of neurons (Cain et al., 2020). While emerging AD case/control single-nucleus RNA sequencing (snRNAseq) datasets offer an exciting opportunity to better resolve such cellular differences (Mathys et al., 2019; Cain et al., 2020; Zhou et al., 2020; Leng et al., 2021), technical constraints have limited the size of such datasets in terms of total numbers of cells and individuals sampled (Park et al., 2021), making it difficult to determine robust cellular differences in a disorder as heterogeneous as AD.

Here we performed a comprehensive analysis of brain bulk- and single-nucleus RNAseq datasets to quantify changes in cell-type proportions in AD. We quantified excitatory and inhibitory neuronal subpopulations and non-neuronal cell types by estimating relative cell-type proportions across three studies and six neocortical brain regions. We corroborated our bulk tissue-based findings by directly estimating cell-type proportions in three snRNAseq datasets collected from AD cases and controls. Finally, we explored how cell-type proportion differences relate to specific age-related neuropathologies, longitudinal cognitive decline, and an established measure of cognitive resilience. Together, our findings suggest a robust and specific loss of excitatory intra-telencephalic neurons and inhibitory somatostatin-expressing interneurons in AD.

Materials and Methods

Studies Used for Bulk Tissue RNA Sequencing Analyses

Post-mortem bulk-brain RNAseq data were processed from 1373 different individuals across three independent studies from the Accelerating Medicines Partnership–Alzheimer’s Disease (AMP-AD) consortium (summarized in Table 1), encompassing six brain regions:

TABLE 1

Summary of RNAseq datasets used in this study.

Dataset name	Study	Data type	Brain region	Total individuals	Controls	AD cases	Other
DLPFC	ROS/MAP	Bulk RNAseq	Dorsolateral Prefrontal Cortex	1092	138	285	669
FP	MSBB	Bulk RNAseq	Frontal Pole	134	44	90	0
IFG	MSBB	Bulk RNAseq	Inferior Frontal Gyrus	112	34	78	0
PHG	MSBB	Bulk RNAseq	Parahippocampal Gyrus	104	38	66	0
STG	MSBB	Bulk RNAseq	Superior Temporal Gyrus	117	35	82	0
TCX	Mayo Clinic	Bulk RNAseq	Temporal Cortex	147	67	80	0
Mathys	ROS/MAP	snRNAseq	Dorsolateral Prefrontal Cortex	47	9	16	22
Zhou	ROS/MAP	snRNAseq	Dorsolateral Prefrontal Cortex	32	8	8	16
Cain	ROS/MAP	snRNAseq	Dorsolateral Prefrontal Cortex	24	5	5	14
AIBS snRNAseq	Allen Institute	snRNAseq	Cingulate Gyrus and Medial Temporal Gyrus	3	3	0	0

Counts of individuals include only samples passing Quality Control.Column of individuals denoted as “Other” indicate individuals from the ROS/MAP cohort not meeting either of the defined control or AD case criteria used in mega-analyses of all cohorts.Note that some individuals from the Mount Sinai Brain Bank, MSBB, were sampled multiple times from different brain regions; therefore, not all samples are independent.

The Religious Orders Study and Rush Memory and Aging Project (herein ROS/MAP) cohort provided bulk RNAseq data for dorsolateral prefrontal cortex (DLPFC) from 1092 individuals. The mean age at death was 89.6 (standard deviation, SD = 6.6).

The Mayo Clinic study (herein Mayo) provided temporal cortex (TCX) samples from 147 individuals. The mean age at death was 82.6 (SD = 8.0).

The Mount Sinai Brain Bank study (herein MSBB) provided samples from the same individuals across multiple brain regions. The mean age at death was 83.3 (SD = 7.4). 134 individuals had bulk-tissue RNAseq data sampled from Frontal Pole (FP), Brodmann area 10; 112 from Inferior Frontal Gyrus (IFG), Brodmann area 44; 104 from Parahippocampal Gyrus (PHG), Brodmann area 36; and 117 from Superior Temporal Gyrus (STG), Brodmann area 22.

Tissue Preparation and Bulk Tissue RNA Sequencing

Details pertaining to the handling and processing of post-mortem samples in ROS/MAP (De Jager et al., 2018), Mayo (Allen et al., 2016), and MSBB (Wang M. et al., 2018) have been previously published (Wan et al., 2020). RNA sequencing procedures differed between studies:

For ROS/MAP, RNA sequencing on DLPFC tissue was carried out in 13 batches within three distinct library preparation and sequencing pipelines (see Supplementary Methods). Sequencing was carried out using the Illumina HiSeq (pipeline #1: 50M 101 bp paired end reads) and NovaSeq6000 (pipeline #2: 30M 100 bp paired end; pipeline #3: 40–50M 150 bp paired end reads). Full details on RNA extraction and sequencing are available on the Synapse AMP-AD Knowledge Portal (syn3219045).

For Mayo, sequencing was carried out on the Illumina HiSeq 2000 (101 bp paired end reads). Details available on the AMP-AD Knowledge Portal (syn5550404).

For MSBB, sequencing was carried out on the Illumina HiSeq 2500 (100 bp single end reads). Details available on the AMP-AD Knowledge Portal (syn3159438).

Processing of Bulk Tissue RNA Sequencing Datasets

Bulk-tissue based RNA-seq read counts from all three studies underwent uniform quality control (QC) and filtering protocols, using the same approach as described previously (Felsky et al., 2020). Briefly, genes with a median expected count less than or equal to 15 were removed and multidimensional scaling was performed. Subjects were deemed outliers and removed if they differed from the sample median of any of the first 5 latent components by more than 3 interquartile ranges. Gene counts were log2 transformed with an offset of 0.5, to coerce any log2(expected count) value differing from the sample median by 3 interquartile ranges to its nearest most extreme value within that range. After sample- and gene-level filtering, the log2(expected counts) were transformed back to the expected count scale. Trimmed mean of m-values (TMM) normalization (using edgeR calcNormFactors) and mean-variance derived observational-level weights were calculated. Variance related to technical factors, including sequencing batch, percent of mapped bases, percent usable bases, RNA integrity number (RIN), and post-mortem interval, were removed using the removeBatchEffect function in limma (Ritchie et al., 2015).

Consensus Definition of Alzheimer’s Disease for Mega-Analysis

We applied a harmonized definition of AD case/control diagnosis as defined previously in the same studies (Wan et al., 2020). This definition combined clinical and post-mortem neuropathological data for categorization, where controls were defined as individuals with a low burden of plaques and tangles, as well as no evidence of cognitive impairment (if available). To define AD case status, Braak stage, global cognitive status at time of death, and CERAD (Consortium to Establish a Registry for Alzheimer’s Disease) scores were used in ROS/MAP, with Clinical Dementia Rating (CDR) scores being used instead of global cognitive status in MSBB. For the Mayo dataset, only neuropathological criteria were used to define case/control status, with details previously published (Allen et al., 2018). In total, 704 individuals across the three studies met the established AD case or control criteria.

Cognitive and Neuropathological Measures in Religious Orders Study/Rush Memory and Aging Project

All subjects in ROS/MAP were administered 17 cognitive tests annually spanning five cognitive domains. Raw scores for tests within each domain were z-scored (using the mean and standard deviation of the entire ROS/MAP combined study population at baseline) and averaged to form the composites. The list of individual cognitive tasks and their corresponding domains has been published (Felsky et al., 2019). Prior to autopsy, the average post-mortem interval was 9.3 h (SD = 8.1). All brains were examined by a board-certified neuropathologist blinded to clinical data. We analyzed 11 neuropathologies measured brain-wide: total amyloid β peptides, neuritic and diffuse plaques, paired helical filament tau protein, neurofibrillary tangles, Braak stage (tau), gross cerebral infarcts, cerebral atherosclerosis, degree of alpha-synucleinopathy, transactive response DNA binding protein 43 kDa (TDP43) proteinopathy, and hippocampal sclerosis. Detailed descriptions of all neuropathological variables have been previously published (Felsky et al., 2019).

Single-Nucleus RNA Sequencing Datasets

In total, we used expression data from four human cortical single-nucleus RNA sequencing (snRNAseq) datasets for this study. First, we used an ultra high-depth SMART-seq based snRNAseq dataset from the human neocortex provided by the Allen Institute for Brain Sciences (AIBS) (Hodge et al., 2019) to define our reference cell type taxonomy and derive cell type specific marker genes (see Supplementary Methods). We used all nuclei sampled from the cingulate gyrus (5,939 nuclei) and medial temporal cortex (15,519), as these correspond most closely with the bulk expression samples described above. Given that nuclei from non-neuronal cell types were relatively undersampled in this dataset, we supplemented this dataset with 2,620 nuclei corresponding to non-neurons sampled from other cortical regions, including visual, auditory, somatosensory and motor cortex (502, 742, 595, and 781 nuclei, respectively). We further used three AD case/control snRNAseq datasets collected from subjects sampled from the ROS/MAP cohort (Mathys et al., 2019; Cain et al., 2020; Zhou et al., 2020). Cells from each of the three ROS/MAP snRNAseq datasets were bioinformatically mapped onto the AIBS snRNAseq dataset (see Supplementary Methods).

Estimation of Relative Cell Type Proportions From Bulk RNA Sequencing Samples

Relative cell type proportions were estimated with the MarkerGeneProfile (MGP) R package, as described previously (Mancarci et al., 2017; Toker et al., 2018), using our derived cell type-specific marker genes with default parameters. The output of the mgpEstimate function was taken as the relative cell-type proportion estimates (rCTPs), providing an indirect measure of cell type abundance in each sample. To ensure consistency in rCTP definitions across individual bulk datasets, rCTPs were estimated using only cell type-specific marker genes passing QC in each of the six bulk-tissue datasets. rCTPs were converted to standardized z-scores within each dataset prior to downstream analysis.

Estimation of Single-Nucleus RNA Sequencing-Derived Cell Type Proportions

Cell type proportions from snRNAseq datasets (snCTPs) were directly estimated from snRNAseq datasets by counting nuclei annotated to each cell type and normalizing by the total count of all QC-passing nuclei per individual subject. We note that such calculations were only performed on nuclei passing quality control and also met our mapping criteria to our reference cell type taxonomy. Direct comparisons between bulk and snRNAseq derived cell type proportions for subjects from the ROS/MAP cohort were performed by identifying subjects in common between both datasets and correlating rCTPs with snCTPs values across subjects.

Statistical Analyses

Mega-Analysis of Bulk RNA Sequencing Cell Type Proportions With Alzheimer’s Disease

The lme4 package in R (Bates et al., 2015) was used to perform a set of mega-analyses, one per cell type, across all bulk RNAseq datasets. Linear mixed effect models were fitted as follows, for each cell type (i), including a random effect of sample to account for correlation structure between brain samples taken from multiple regions of the same individuals in the MSBB study:

Adjustment of two-sided p-values to account for multiple cell types was performed. The Bonferroni method and false discovery-rate (FDR) method (Benjamini and Hochberg, 1995) were applied separately to identify highly stringent and more relaxed thresholds for statistical significance. Corrected p-values are labeled specifically within results (i.e., p, p).

Analysis of Single-Nucleus Cell Type Proportions in Alzheimer’s Disease and Controls

For snCTPs, AD association analyses were performed as for rCTPs, with an additional covariate added for post-mortem interval (PMI; for bulk analyses, variation in gene expression due to PMI was removed during the preprocessing phase) and using a linear model due to the limited overlap in individuals sampled between snRNAseq datasets.

Association of Bulk Tissue Relative Cell-Type Proportions With Neuropathology, Cognitive Decline, and Residual Cognition in Religious Orders Study/Rush Memory and Aging Project

In ROS/MAP we performed detailed analyses associating each rCTP to measures of individual brain pathologies, global cognitive decline, and cognitive status proximal to death. We specifically utilized the full set of individuals in ROS/MAP with bulk expression and other measures available, as opposed to the more limited set of individuals in our cross-cohort mega-analysis of AD case/control criteria described above. “Residual cognition” was calculated per individual as in White et al. (2017); the measure corresponds to residuals of a linear model regressing global cognition proximal to death on observed neuropathological variables and demographic factors. Associations between rCTPs for each cell type and cognitive and neuropathological phenotypes were assessed using linear models covarying for age at death, sex, and PMI. For models of cognitive outcomes, we covaried for sex, educational attainment, and age at time of baseline study assessment. Correction for multiple testing across all cell types and outcomes (19 × 14 = 266 tests) was performed using the FDR method. To estimate variance explained (R2) by rCTPs over and above covariates and measured neuropathologies, we built a series of nested linear models and compared them using likelihood ratio tests. To improve the generalizability of our estimates, models were bootstrapped (100 iterations) using the 0.632 + method (Efron and Tibshirani, 1997).

Causal Mediation Modeling of Intra-Telencephalic and Somatostatin Relative Cell-Type Proportions, Alzheimer’s Disease Neuropathology, and Cognitive Performance in Religious Orders Study/Rush Memory and Aging Project

The R mediation package (v4.5.0) was used for causal mediation modeling. To model pathological burden, we used an established measure of global post-mortem AD neuropathology: the mean of brain-wide diffuse plaques, neuritic plaques, and neurofibrillary tangles. Four models were tested with SST and IT rCTPs as either predictor or mediator, AD pathology the other role and global cognitive performance at last study visit was always the outcome measure. To estimate confidence intervals for average indirect, direct, and total effects, 1000 Monte Carlo draws were used for non-parametric bootstrapping.

Results

Deriving Cell Type-Enriched Marker Genes for the Major Neuron Subclasses of the Human Neocortex

To build a high-quality foundation for investigating cell subtype proportions–focusing on subclasses of neocortical neurons–in AD, we first established representative marker gene sets using ultra-high depth single-nucleus data RNA sequencing (snRNAseq) data from multiple regions of human cortex collected by the Allen Institute for Brain Sciences (Hodge et al., 2019). In these datasets, nuclei were annotated to an established cell type reference taxonomy where transcriptomically defined cell clusters are linked to orthologous multi-modal taxonomies established in other species and neocortical regions (Tasic et al., 2018; Gouwens et al., 2020; Berg et al., 2021). This annotation enables the inference of additional aspects of cellular identity for these human snRNAseq-based cell clusters that include axonal projection patterns and cell morphology (Hodge et al., 2019).

We used these snRNAseq reference data to derive cell type-specific “marker genes” (illustrated in Supplementary Table 1), focusing our analyses primarily on the subclass cell type resolution. This taxonomic grouping serves as an intermediate resolution (e.g., somatostatin-expressing GABAergic interneurons) between more coarse-grained (e.g., inhibitory neurons) and fine-grained cell taxonomic divisions (e.g., Martinotti neurons). A key benefit of these markers is their specificity and good cross-dataset replicability, including in snRNAseq datasets specific to aging and AD samples (Supplementary Figures 1, 2).

Bulk Tissue Analysis Implicates Fewer Inhibitory and Excitatory Neurons in Alzheimer’s Disease, With Most Robust Associations With Somatostatin Interneurons and Intra-Telencephalic Pyramidal Cells

We first sought to understand how the abundance of specific cell types are different in brains of individuals with a pathologically confirmed AD diagnosis compared to controls. We estimated the relative cell type proportions (rCTPs) of each post-mortem bulk tissue RNAseq sample across all six bulk expression datasets separately using the Marker Gene Profile (MGP) method (Mancarci et al., 2017; Toker et al., 2018) and our novel cell type-enriched marker gene sets described above. We then performed mega-analysis for AD case/control status with rCTPs across each of the major cell subclasses and all six datasets using a linear mixed effects model. In aggregate, we found lower rCTPs in most GABAergic subclasses, mostly fewer but some greater rCTPs among glutamatergic subclasses, and higher rCTPs for most non-neuronal subclasses (Figure 1A).

FIGURE 1

Differences in relative cell type proportions of neuronal and non-neuronal subclasses in Alzheimer’s disease. (A) Consensus associations of Alzheimer’s disease (“AD”) vs. control (“C”) status and cell type-specific relative cell type proportions (rCTPs) across six bulk RNAseq datasets. Y-axis shows standardized beta coefficients estimated using a mixed effects model to pool associations across datasets. Positive (negative) standardized beta coefficients indicate an increase (decrease) in the cell type-specific rCTP in AD. Error bars indicate one standard deviation. Asterisks (dots) above each cell type indicate two-sided p < 0.05 (or less stringent FDR < 0.1). (B) Comparisons between rCTPs between controls and AD cases in each of the six bulk gene expression datasets, ROS/MAP, sampling the dorsolateral prefrontal cortex (DLPFC), MSSB, sampling the Frontal Pole (FP), Inferior Frontal Gyrus (IFG), Parahippocampal Gyrus (PHG), and Superior Temporal Gyrus (STG), and the Mayo cohort, sampling the Temporal Cortex (TCX). Y-axis shows estimates of rCTPs for somatostatin (SST) interneurons from individual post-mortem samples (each dot reflects one individual), after covarying for age and sex. Numbers show p-values from a statistical model comparing residualized rCTPs between controls and AD cases, uncorrected for multiple comparisons across datasets and cell types. Subjects with outlier values of rCTPs not shown. (C) Heatmaps showing AD case/control associations for marker genes for SST inhibitory cells. White dots indicate specific associations where FDR < 0.1. (D,E) Same as panels (B,C) for intra-telencephalic-projecting (IT) excitatory pyramidal cells (IT cells).

Specifically, at a stringent Bonferroni correction threshold of p < 0.05, our analysis identified lower rCTPs for SST (β = −0.48, p = 8.98 × 10–9) and PVALB (β = −0.28, p = 0.0072) GABAergic interneurons, as well as IT (β = −0.45, p = 4.32 × 10–7), L4 IT (β = −0.24, p = 0.023), L5 6 NP (β = −0.23, p = 0.039), and L6 CT (β = −0.25, p = 0.023) glutamatergic neurons in AD. At the same threshold, we observed greater rCTPs for L5 6 IT Car3 glutamatergic neurons (β = 0.023, p = 0.034) and VLMC cells (β = 0.24, p = 0.041). At a less stringent threshold (FDR < 0.1), we also observed lower rCTPs for lysosomal associated membrane protein family member 5 (LAMP5) and vasoactive intestinal peptide-expressing (VIP) GABAergic interneurons, and greater rCTPs for L6b glutamatergic cells and most non-neuronal cells, except endothelial cells. One important caveat of these analyses is the focus on relative proportions, which are not absolute cell counts (Mancarci et al., 2017; Toker et al., 2018); therefore, potentially paradoxical reported differences in some rCTPs here, such as greater L5 6 IT Car3 and L6b glutamatergic cell rCTPs in AD, may not necessarily indicate that these cell types are actually increasing in their absolute cell counts. These results are consistent with prior observations that AD is characterized by relatively fewer neuronal cells and corresponding relatively more non-neuronal cells (Patrick et al., 2020).

The strongest AD-associated cell type in mega-analysis was the somatostatin (SST) interneuron (β = −0.48, p = 8.98 × 10–9); in each individual dataset, SST rCTPs were lower in AD cases relative to controls (Figure 1B), though the differences were not significant in all regions. Our findings mirror those of Cain et al. (2020) highlighting SST interneurons as particularly associated with AD phenotypes among ROS/MAP participants and further generalizes this finding to additional studies and brain regions. Importantly, in addition to the SST gene mRNA transcript itself, we also observed lower mRNA expression of other SST interneuron marker genes, including GRIK1 and COL25A1 across most brain regions (Figure 1C). Moreover, the SST rCTP signal is robust, albeit attenuated, to the removal of the SST mRNA as a marker gene (β = −0.39, p = 6.02 × 10–6), suggesting the relevance of fewer SST-expressing neurons rather than a lower expression of the SST gene specifically.

Among excitatory cell types, rCTPs for intratelencephalic- projecting (IT) pyramidal cells showed the greatest difference between AD and controls (β = −0.45, p = 4.32 × 10–7). Like SST rCTPs, proportionally fewer IT neurons were observed in AD cases relative to controls in each of the six bulk expression datasets (Figure 1D), albeit not significantly in all regions. The IT cell subclass includes both superficial layer pyramidal cells, such as CUX2-positive cells, as well as more infragranular cells, including RORB- and THEMIS-positive neurons (Hodge et al., 2019). As expected, we observed lower expression in many of the individual IT cell marker genes across each of the datasets in AD, including LINC00507 and LINC01202 (Figure 1E).

Single-Nucleus Analysis Suggests Selective Vulnerability of Specific Inhibitory Subclasses in Alzheimer’s Disease, Including Lysosomal Associated Membrane Protein Family Member 5 and Somatostatin Interneurons, but Not Intra-Telencephalic-Projecting Pyramidal Cells

To complement the indirect estimates of rCTPs from bulk tissue, we accessed three AD case/control snRNAseq datasets sampled from subsets of participants from the ROS/MAP cohort (Table 1; Mathys et al., 2019; Cain et al., 2020; Zhou et al., 2020). We first harmonized cell type annotations from these snRNAseq datasets by mapping cells to the same Allen Institute-based human cortical cell type reference taxonomy used in our rCTP analyses, finding good qualitative agreement between the annotated cell type identities provided within the original publications and those following QC and re-mapping (Supplementary Figure 3). We then estimated single-nucleus CTPs (snCTPs) per post-mortem sample by counting nuclei annotated to each cell type, as a percentage of the total nuclei sampled in each subject. Correlations between cell type-specific snCTPs and bulk-tissue derived rCTPs in subjects with overlapping data types were modest, but overall positive (Supplementary Figure 4), in agreement with published deconvolution efforts validated using immuno-histochemistry (Patrick et al., 2020).

In line with our bulk tissue rCTP analysis, a mega-analysis across the three snRNAseq datasets indicated that AD samples showed lower snCTPs in most inhibitory subclasses, both higher and lower snCTPs among excitatory subclasses, and greater snCTPs for most non-neuronal subclasses (Figure 2A). At Bonferroni-corrected p < 0.05, we found only LAMP5 GABAergic interneurons to be lower in AD compared to controls (β = −0.95, p = 0.011). At a less stringent threshold (FDR < 0.1), we also observed lower snCTPs for PAX6 (β = −0.62, FDR = 0.093) and SST (β = −0.74, FDR = 0.093) interneurons (Figure 2B). We note that we did not see strong evidence for lower expression of SST mRNA among SST-annotated nuclei in AD after covarying for SST cell proportion differences (β = −0.43, p = 0.47; Supplementary Figure 5), providing additional evidence for SST cell-specific vulnerability highlighted by bulk tissue analysis. In contrast to our bulk tissue results, we did not find any effects for IT pyramidal cells (Figure 2C). To assess the overall consistency between our bulk tissue rCTP and single-nucleus approaches, we correlated standardized effect coefficients for each cell type between analyses (Figure 2D). Effects were strongly correlated (Spearman rho = 0.65), illustrating broad concordance between methods, with some exceptions.

FIGURE 2

Differences in single-nucleus derived cell type proportions (snCTPs) of neuronal and non-neuronal subclasses in Alzheimer’s disease. (A) Consensus associations of AD status and snCTPs across three AD snRNAseq case/control datasets. Y-axis shows standardized beta coefficients estimated using a mixed effects model to pool associations across datasets. Positive (negative) standardized beta coefficients indicate an increase (decrease) in the cell type-specific snCTP in AD. Error bars indicate one standard deviation. Asterisks (dots) above each cell type indicate two-sided p < 0.05 (or less stringent FDR < 0.1). (B) Comparisons between snCTPs between controls and AD cases in each of three snRNAseq datasets. Y-axis dots show snCTPs for somatostatin (SST) interneurons (as a percentage of all nuclei sampled) from individual post-mortem samples. Numbers show p-values from t-test (uncorrected for multiple comparisons across datasets and cell types) from a statistical model comparing snCTPs between controls and AD cases. Subjects with outlier values of rCTPs not shown. (C) Same as panel (B) for intratelencephalic-projecting (IT) excitatory pyramidal cells. (D) Consistency of AD standardized effect sizes between bulk rCTPs and snCTPs based on single-nucleus analyses. X-axis shows point estimates of standardized beta coefficients of AD effects on rCTPs in the ROS/MAP cohort (as in Figure 1A) and y-axis is the same as point estimates shown in A. Diagonal line denotes the unity line. Inset correlation value denotes overall Spearman correlation (rho) between rCTP and snCTP estimated effects.

Somatostatin Interneurons and Intra-Telencephalic-Projecting Pyramidal Cells Specifically Are Correlated With Alzheimer’s Disease Neuropathologies and Residual Cognition

Having identified SST interneurons and possibly IT-projecting pyramidal cells as especially diminished in AD cases vs. controls, we explored the associations of rCTPs with specific aging-related neuropathologies and rates of longitudinal cognitive decline. We utilized a larger set of individuals from ROS/MAP with available data (889 subjects), as opposed to only those meeting the consensus AD case/control criteria. After FDR correction, we observed a striking pattern of association whereby only SST rCTPs were negatively associated with each tau and beta-amyloid-related neuropathology–most strongly with neuritic plaques (p = 3.1 × 10–4)–and positively associated with rates of cognitive decline (p = 3.9 × 10–6) and cognition measured proximal to death (p = 5.7 × 10–5) (Figure 3). In contrast, IT neurons were also associated with both cognitive measures (decline: p = 8.3 × 10–5; proximal to death: p = 1.2 × 10–7) but were not associated with canonical AD-related neuropathologies. However, an association with TDP-43 proteinopathy was observed (p = 0.015). At a relaxed FDR < 0.1 threshold, several other neuronal and non-neuronal associations were observed for individual pathologies (Figure 3), though none as strong as those for SST and IT neurons.

FIGURE 3

Associations between cell type specific relative proportions and neuropathologies, cognition, and residual cognition. Inset values denote the FDR statistics of specific associations, where FDR < 0.1. Note that while pathology scores are coded such that greater levels of pathology indicate worsening brain health, global cognition scores are coded such that higher scores indicate better cognition and less dementia. Std Beta, standardized beta coefficients; PHF, paired helical filaments; AA, amyloid angiopathy.

Finally, we sought to determine if cell type-cognition associations were either independent or a reflection of accumulating brain pathology. Therefore, we calculated a measure of residual cognition for all individuals, as described previously (White et al., 2017), which represents global cognitive performance (proximal to death) after accounting for variability due to neuropathology and demographics (see Section “Materials and Methods”). After correction, IT rCTPs were the only cell type significantly associated with residual cognition (p = 1.2 × 10–5), though we note SST neurons showed a marginal association as well (p = 0.069). To quantify the additional variance in cognition explained by IT rCTPs over and above measured neuropathology, we first established a baseline model of cognition, where demographic and neuropathological variables alone explained 40.3% of total variance (R2 = 0.403). Adding IT rCTPs to this model increased the variance explained by 1.9% (likelihood ratio test p = 7.1 × 10–8, R2 = 0.422). By contrast, SST rCTPs increased the explained variance to a much lesser extent (additional 0.53%; likelihood ratio test p = 0.0049). These findings were supported by mediation analyses, including apolipoprotein E ε4 AD risk variant as a co-variate, which found bidirectional mediation of SST and AD pathology on cognitive performance proximal to death (p < 0.0001), but no mediation of the relationship between pathology and cognition by IT neurons (p = 0.31), or of IT neurons and cognition by pathology (p = 0.35) (results in Supplementary Table 2).

Discussion

Our analysis leveraged three aging and AD studies with multi-region post-mortem bulk gene expression to determine which neocortical cell subtypes are most strongly associated with AD. Based on marker gene expression specific to cellular subclasses, we observed lower relative proportions of most neuronal types and greater relative proportions of most non-neuronal types. In particular, our analyses highlighted fewer somatostatin-expressing (SST) interneurons and intra-telencephalic projecting (IT) pyramidal cells in AD that were replicated across studies and neocortical regions. Cellular proportions directly derived from three additional AD case/control single-nucleus RNAseq datasets provided partial corroboration of our bulk-tissue based results, suggesting that such cellular changes are likely the result of cellular loss as opposed to a coordinated global change in cellular identity. The results of our analyses support previous literature implicating the loss of SST interneurons in AD and further indicate that the preservation of IT pyramidal cells may contribute to cognitive resilience despite the presence of AD neuropathologies.

Our conclusion implicating SST interneurons, a key subpopulation of cortical GABAergic cells that provide synaptic inhibition to pyramidal cells by targeting their distal dendrites (Urban-Ciecko and Barth, 2016), are consistent with a broad literature on the role of SST in neurological illness, recently reviewed (Song et al., 2021). The association of SST with AD is decades-old, beginning with seminal findings reporting reduced levels of somatostatin immunoreactivity in AD brain (Davies et al., 1980), and more recently with cross-study differential expression analyses finding ubiquitous reductions of SST RNA in AD brain tissue (with the exception of the cerebellum[1]). However, it remains unknown if this association is driven by a selective loss of SST mRNA or losses of populations of SST-expressing neurons. Taken as a whole, our bulk and single-nucleus based findings support the latter conclusion, though certainly do not provide definitive evidence for it. In agreement, one study found that SST interneurons were uniquely lost in AD among neuronal subtypes in prefrontal cortex ROS/MAP samples (Cain et al., 2020). While the precise role of selective SST interneuron vulnerability in AD remains to be understood, a recent publication pointed to a role for a potential direct biochemical interaction between the SST neuropeptide and amyloid beta (Solarski et al., 2018).

In addition, we observed negative associations between four neuronal subtypes (GABAergic: LAMP5, VIP, and PVALB; glutamatergic: IT) and TDP-43 neuropathology, with the most consistent effects among GABAergic neurons. We are cautious to not over-interpret this result given its non-specificity and relatively weak statistical signal, but note that GABA-specific involvement in TDP-43 neurodegeneration has been observed in preclinical (Tsuiji et al., 2017) and human contexts (Lin et al., 2021). In addition, recent work has shown that gene expression modules enriched for GABAergic neurons (module 18) are strongly regulated by the TMEM106B genetic locus, a known risk factor for TDP-43 neurodegenerative conditions (Yang et al., 2020).

We also observed strong associations between IT pyramidal cells and AD, and, intriguingly, this was the only cell type significantly associated with residual cognition. IT cells are defined by their cortico-cortical and cortico-striatal projections (Tasic et al., 2018) and encompass supragranular pyramidal cells, such as CUX2-positive cells, and infragranular cells, including RORB- and THEMIS-positive pyramidal cells (Hodge et al., 2019). Immunohistochemical studies corroborate these results in part, suggesting that SMI32-immunoreactive neurons, labeling large pyramidal neurons in Layers 3 and 5, are selectively lost in the frontotemporal cortex in AD (Hof et al., 1990; Bussière et al., 2003; van de Nes et al., 2008). More recently, evidence from snRNAseq studies of the entorhinal cortex and superior frontal gyrus have implicated RORB-expressing excitatory neurons as selectively vulnerable in AD (Leng et al., 2021). As one caveat, we note that we saw some, but overall limited replicability of decreased IT CTPs in our bulk- compared to our single-nucleus analyses.

Our study has several key limitations and considerations. First, the backbone of our study is a neocortical cell type taxonomy derived from deep transcriptomic sequencing of single-nuclei from normotypic individuals (Hodge et al., 2019); it remains unclear how comparable these transcriptional profiles are to those in the elderly and in individuals with AD. In addition, the fact that our marker gene sets were required to pass quality control in each dataset means that some biological cell-specific signals may have been missed, though the observed consistency in effects between our individual study samples (Figures 1B–E) is encouraging. Second, the conclusion of our study rely on the accurate estimation of relative differences in cell-type proportions between individuals. Such estimates are highly dependent on the choice and quality of the constituent marker genes that serve as representatives of each cell type (Mancarci et al., 2017; Toker et al., 2018) as well as the particular choice of method for cellular deconvolution (Jew et al., 2020; Patrick et al., 2020; Park et al., 2021). However, the relatively high degree of consistency between our bulk- and single-nucleus based results help partially mitigate this concern. Third, the focus of this work is the study of robust changes in cell-type proportions in AD, but does not tackle the question of within-cell type transcriptional regulation (Mathys et al., 2019; Wang and Li, 2021). Lastly, all of the results presented here are associational; further studies are needed to determine how and when cell type-specific loss occurs relative to the emergence of AD pathologies and cognitive decline.

Conclusion

Overall, our study provides a comprehensive consensus overview of the vulnerability of neocortical neuronal subpopulations in AD. Our results demonstrate that losses of SST interneurons and IT pyramidal cell populations are those most strongly associated with AD. In addition, IT pyramidal cells are uniquely associated with residual cognition, suggesting that efforts to preserve or maintain this key neuronal subpopulation might mitigate cognitive decline in the face of AD neuropathologies. Our hope is that these results will inform future studies to further disentangle the causal progression of AD neuropathological burden, cellular loss, and cognitive decline.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The Synapse AMP-AD Knowledge Portal (https://adknowledgeportal.synapse.org/ and doi: 10.7303/syn2580853). ROSMAP resources can be requested at https://www.radc.rush.edu.

Ethics Statement

The studies involving human participants were reviewed and approved by the Rush University, Mayo Clinic, Mount Sinai, and JJ Peters VA Medical Center Institutional Review Boards. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MC was responsible for data processing, statistical analysis, and manuscript writing and editing. YC contributed to statistical analyses, data visualization, and manuscript writing. DF and ST were responsible for data access, ensuring data quality control, study design, and manuscript writing and editing. VM, YW, PD, DB, and JS were responsible for aspects of data collection, collaborative input on study design, and manuscript editing. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.