Introduction

Changes in gene expression characterize a multitude of human diseases and have been successfully used to predict molecular and cellular mechanisms associated with pathological processes1. Alzheimer’s disease (AD) is the most prevalent type of dementia and causes a progressive cognitive decline, for which there is no effective treatment or cure. Although expression analyses in brain diseases are generally limited by tissue availability, RNA sequencing (RNAseq) data have been generated from postmortem brain samples of healthy and AD individuals2,3,4. However, a comprehensive description of the gene expression alterations in the AD brain remains elusive.

Recent work has begun to address this important gap in the study of AD pathology using bulk brain tissue RNA sequencing (RNAseq)5 or single-cell RNA sequencing (scRNAseq)6,7. However, these studies have focused on samples obtained from different brain regions, namely the dorsolateral prefrontal5,7 and entorhinal cortices6, which could lead to important discrepancies in the results. In fact, AD pathology shows a progressive impact on different brain regions, characterized at early stages by the presence of TAU protein inclusions in the locus coeruleus, the transentorhinal and entorhinal regions (stages I and II). This is followed by the presence of TAU inclusions in the hippocampal formation and some parts of the neocortex (stages III and IV), followed by large parts of the neocortex (stages V and VI)8. This temporal progression of AD pathology could differently impact gene expression in those brain areas. Accordingly, a recent study has shown that changes in protein expression are much more prominent in areas affected at early and intermediate stages, such as the hippocampus, entorhinal cortex, and cingulate cortex in the temporal lobe, compared to other brain regions affected at later stages of AD pathology, such as sensory cortex, motor cortex, and cerebellum9.

Another important aspect to consider is the descriptive relevance of gene expression analysis based solely on the identification of differentially expressed genes (DEG), which fails to detect dynamics in the expression of multiple related transcripts10. Recently, new approaches using transcripts-level analysis, so-called differential transcript usage (DTU), enables identification of alternative splicing and isoform switches with the prediction of functional consequences11,12. Therefore, important gene expression modifications in the AD brain could occur at the transcript level and be overlooked in classical DEG analyses.

Here, we took advantage of three available RNAseq datasets, generated using samples from different brain regions, to systematically probe gene expression changes (DEG and DTU) in AD. In Mayo’s clinic study, both the temporal cortex and cerebellum were used to obtain bulk RNAseq2. In the Religious Orders Study (ROS) and Memory and Aging Project (MAP), henceforth called ROSMAP, the dorsolateral prefrontal cortex was used3. Finally, in the Mount Sinai/JJ Peters VA Medical Center Brain Bank (MSBB), four different Brodmann areas of the brain were studied: areas 22 and 36 from the temporal lobe, areas 10 and 44 in the frontal lobe4. We also added another level of complexity using an indirect approach to assign DEGs and gDTUs to unique cell types in order to identify AD gene expression signatures for neural cells, microglia, and endothelial cells. Finally, we linked these alterations with AD causal and risk genes, identifying novel isoform switches in BIN1 and APP genes of potential functional consequences for pathology progression.

Results

Regional gene expression alterations in the AD brain correlates with pathological progression

Several consortia have generated RNAseq data from brains of individuals with a clinical and/or pathological diagnostic of AD2,3,4. Considering the regional progression of AD pathology8, we set out to identify and compare differentially expressed genes (DEG) in the temporal lobe (TL), encompassing brain regions affect at early stages of the AD such as the hippocampus and entorhinal cortex, and in the frontal lobe (FL), affect at more advanced stages of the pathology (Fig. 1). Comparisons between control and AD individuals were performed independently for each dataset and only genes with fold change >1.3 and FDR > 0.01 were considered as DEGs. We found 3348 (1244 down- and 2104 upregulated genes) and 2172 (1170 down and 999 upregulated genes in BM22 and BM36; three genes regulated in opposite directions in these two areas) DEGs in the TL of AD individuals compared to their respective controls in the MSBB_TL and Mayo datasets, respectively (Fig. 2A, B and Supplementary Table 1). Of those DEGs, 734 genes (145 down and 520 up) were commonly regulated in both Mayo and MSBB_TL (88.4% of genes altered in the same direction; 15.33% of overlap; P = 8.56 × 10−59, hypergeometric test). In contrast, only 327 (113 down and 214 up) and 209 (97 down and 112 up) DEGs were detected in the MSBB_FL and ROSMAP, respectively. Of those, 31 genes (18 down and 13 up) were found in both datasets (7.34% of overlap; P = 1.67 × 10−14, hypergeometric test) (Fig. 2A, B and Supplementary Table 1). This small number of DEGs in the FL is in agreement with previous data obtained from the DLPFC (106 down- and 158 upregulated genes with FC > 1.3)13. Among DEGs detected in the FL, 62.5% were also detected in the TL (Fig. 2B), suggesting that similar molecular changes occur in these brain areas, but at different stages of the disease progression. The differences in the number of DEGs detected in the FL and TL can neither be attributed to lack of statistical power nor potential biases due to tissue processing, since the number of samples in the FL is larger than in the TL groups (Fig. 1) and differences are observed even in samples obtained from the same donors (compare MSBB_TL and MSBB_FL in Fig. 2). Thus, changes in gene expression are much more prominent in brain areas affected at the early stages of AD pathology.

Fig. 1: Schematic summary of the methodology.
figure 1

A Datasets obtained from three consortia (Mayo, MSBB, and ROSMAP) were grouped according to the brain region sampled in the frontal lobe (FL) or temporal lobe (TL). Next, RNAseq data were pseudo-aligned using Kallisto. Clinico-pathological classifications were included as metada. B scRNAseq data from the midle temporal gyrus (MTG, Allen Brain Atlas) were analyzed using the R package SEURAT. C Gene expression analyses were performed using the R packages DESeq2, IsoformSwitchAnalyzeR (ISAR), and gene set enrichment analysis (GSEA). Assignment of differentially expressed genes or isoform switches to specific cell types/subtypes was performed indirectly using scRNAseq signatures obtained from the MTG (B).

Full size image
Fig. 2: Gene expression alterations are more prominent in the temporal than the frontal lobe of AD patients.
figure 2

A Volcano plots showing differentially expressed genes (DEG, red dots; FC > 1.3 and FDR < 0.01) in the frontal lobe (ROSMAP and MSBB FL - BM10 and BM44) and temporal lobe (Mayo and MSBB TL - BM22 and BM36). B Upset plot showing the total number of DEGs identified in each dataset (horizontal bars) and the number of DEGs exclusive of one dataset (first four vertical bars) or shared by different datasets (other vertical bars). Black dots below vertical bars indicate datasets quantified. Venn diagram illustrates the same results in colors and circle sizes. C Gene ontology terms enriched for DEGs identified in the TL or FL intersections (TLI and FLI, respectively).

Full size image

To select genes consistently altered in AD brains, considering the several sources of measurement variations in RNAseq experiments14, we decided to focus only on DEGs replicated in at least two independent datasets obtained from related brain areas. This resulted in a set of 734 DEGs detected in both Mayo and MSBB TL (temporal lobe intersection—TLI), and 31 DEGs shared between ROSMAP and MSBB FL (frontal lobe intersection—FLI) (Supplementary Table 2). Among TLI DEGs, we observed ABCA1 and 2 (ATP-binding cassette subfamily A member 1 and 2), primarily involved in the maintenance of normal brain homeostasis and associated with AD and other neurological diseases15; Complement C1R and C1S, involved in the immune/inflammatory response and previously shown to be upregulated in the brain of a 3 × Tg mouse model of AD when Aβ plaques start to accumulate16; RE1 silencing transcription factor (REST), which regulates neural circuit activity during aging17; glutamate decarboxylase 1 and 2 (GAD1 and 2), solute carrier family 32 GABA vesicular transporter, member 1 (SLC32A1), calbindin 1 (CALB1), parvalbumin (PVALB), somatostatin (SST), and vasoactive intestinal peptide (VIP), all expressed in GABAergic neurons and involved in cognitive decline in AD and other neurological diseases18. Among the few DEGs common to TLI and FLI, we observed a significant downregulation of the neurosecretory protein VGF (VGF nerve growth factor inducible), recently suggested as a key regulator of Alzheimer’s disease19.

Next, we used gene set enrichment analyses (GSEA) to assess the functional profile of the DEGs identified in our analysis. Again, we used only genes commonly altered in two datasets (TLI or FLI) to avoid inaccurate results associated with the use of large gene sets in functional analysis20. We found that TLI DEGs were significantly enriched for terms (GO:BP, GO:CC and KEGG) associated with generic biological processes, such as cell-signaling pathways and cell-cell signaling, whereas the small number of DEGs in the FLI were not significantly enriched for any term (Fig. 2C and Supplementary Table 3). The limited number of significant gene set enrichment observed in our analysis after inputting DEGs is in disagreement with results reported by Canchi et al.13. This discrepancy can likely be explained by the use of stringent criteria to detect TLI DEGs in our study (only genes detected in at least two independent datasets with FC > 1.3 and FDR < 0.01), which significantly reduce the number of genes used in the GSEA.

Differential transcript usage analysis reveals novel genes associated with AD pathology

Gene-level expression analysis lacks the sensitivity to detect possible changes at the transcript-level caused, for example, by alterations in alternative splicing10,21. To overcome this limitation, we used differential transcript usage (DTU) analysis to identify additional alterations of gene expression in the AD brains compared to controls. We observed 2509 and 1843 genes with differential transcript usage (gDTU) in the temporal lobe of AD brains studied in the Mayo and MSBB datasets, respectively (Fig. 3A, B and Supplementary Table 1). Similar to what we observed for DEGs, a much smaller number of gDTUs were detected in the frontal lobe, both in ROSMAP and MSBB studies (59 and 855 genes with transcripts altered, respectively). We found 435 gDTUS in TLI (11.1% of overlap; P = 6.16 × 10−25, hypergeometric test) and 13 gDTUs in FLI (1.47% of overlap; P = 2.56 × 10−3, hypergeometric test) (Supplementary Table 2). In TLI, most gDTUs did not overlap with DEGs (TL—34 gDTUs that are DEGs out of 435 gDTUs, Fig. 4A), whereas in FLI, we found no overlap at all. Consistent with this small overlap, GSEA using only DEGs, only gDTUs or both showed complementary results (Fig. 4B). GSEA using gDTUs (alone or in combination with DEGs) showed significant enrichment for vesicle-mediated transport and other synapse-related terms, which were not observed while inputting only DEGs (Figs. 3C and 4B; Supplementary Table 3). The functional enrichment annotation using both DEGs and gDTUs is in agreement with previous studies using scRNAseq to identify gene expression alterations in unique cell types6,7 and clearly improves the annotation observed using only DEGs, suggesting that the use of DTU analysis could contribute to unraveling gene expression alterations overlooked in the classical DEG analysis.

Fig. 3: Differential transcript usage analysis identifies gene expression alterations in AD associated with synapse transmission.
figure 3

A Volcano plots showing genes with differential transcript usage (gDTU, yellow dots; Differential isoform fraction (dIF) >0.05 and FDR < 0.05) in the frontal lobe (ROSMAP and MSBB_FL - BM10 and BM44) and temporal lobe (Mayo and MSBB_TL - BM22 and BM36). B Upset plot showing the total number of gDTUs identified in each dataset (horizontal bars) and the number of gDTUs exclusive of one dataset (first four vertical bars) or shared by different datasets (other vertical bars). Black dots below vertical bars indicate datasets quantified. Venn diagram illustrates the same results in colors and circle sizes. C Synapse-related terms enriched for gDTUs in the TLI are not observed in the FLI.

Full size image
Fig. 4: Differential transcript usage analysis in AD brains reveals gene expression alterations overlooked in DEG analysis.
figure 4

A Venn diagram showing DEGs and gDTUs identified in the TLI. B Comparison of GO and KEGG terms enriched for DEG, gDTU, or DEG + gDTU identified in the TLI.

Full size image

Among genes with isoform switches enriched in synaptic-related terms, we observed the AD causal gene APP, previously associated with regulation of synapse transmission and long-term plasticity in AD22; neuronal vesicle trafficking associated 1 (NSG1), which has been implicated in the regulation of AMPA receptors (AMPAR) and APP trafficking, thus affecting synaptic transmission, plasticity, and Aβ production23,24; RELN (Reelin) that plays important role in synaptic transmission and has been associated with AD25; gamma-aminobutyric acid type A receptor subunit alpha 1 (GABRA1), which encodes for a subunit of the main ionotropic GABA receptor in the brain and has previously been shown to be downregulated in the AD brain26.

Alternative splice events in AD brains and functional consequences

To identify the causes subjacent to gene isoform switches in the AD brain, we quantified the frequency of splicing events associated with the isoform switches detected in AD compared to control brains (Fig. 5 and Supplementary Table 4). We found that alternative transcription start site (ATSS), alternative transcription termination site (ATTS), and exon skipping (ES) were the most frequent splicing events in AD brains (Fig. 5B). Other common splicing events observed were alternative 3′ or 5′ splice sites (A3 and A5, respectively), multiple exons skipping (MES) and intron retention (IR) (Fig. 5B). These observations suggest that changes in alternative splicing could be implicated in AD pathogenesis, corroborating previous analyses in the ROSMAP cohort using intronic usage ratios to identify abnormal splicing events in the AD brain5.

Fig. 5: Alternative splicing mechanisms associated with isoform switches and consequences for protein expression.
figure 5

A Schematic showing different splicing events that can lead to gene isoform switches. B Quantification of the number of isoforms showing more or less splicing events in AD compared to controls for each dataset. C Quantification of the number of isoforms showing (i) gain or loss of coding potential, domains/signal peptides identified, intrinsically disordered regions (IDR), intron retention, open-reading frame (ORF) sequencing similarity; (ii) switch (simultaneous gain and loss) of domains identified or IDR; (iii) sensitive or insensitive to nonsense-mediated decay (NMD); and (iv) longer or shorter ORF sequencing similarity.

Full size image

Alternative splicing events may have diverse functional consequences for protein expression, such as shifting the frequency of transcripts containing introns (noncoding) or mRNA stability (nonsense-mediated decay) or leading to gain/loss of protein domains, intrinsically disordered regions, or signaling peptides12. Quantification of these consequences revealed some interesting differences between Mayo and MSBB BM36 (Fig. 5C), the two datasets with the largest numbers of gDTUs. Whereas in the Mayo dataset, a high number of isoforms showed loss of coding potential and protein domains, in the MSBB BM36 isoforms showed an even distribution of loss and gain of coding potential or protein domains (Fig. 5B). These differences could be at least partly explained by the larger number of gDTUs detected in the Mayo compared to MSBB TL (Fig. 3) and are likely related to the different median read depth of these datasets (Mayo—12.58 billion bases; MSBB BM22—3.23 billion bases; MSBB BM36—3.56 billion bases)27.

Differential expression of genes involved in alternative splicing correlates with isoform switches during disease progression

To evaluate whether the emergence of gDTUs could be correlated with AD pathology hallmarks, we quantified the total of gDTUs observed at different disease stages in the MSBB dataset using the Braak classification (Fig. 6 and Supplementary Table 5). For this purpose, we subdivided samples into three groups: low Braak (0, 1, and 2)— 196 samples (clinical diagnosis: 15 AD and 181 controls); mid-Braak (3 and 4)—133 samples (clinical diagnosis: 58 AD and 75 controls); and high Braak (5 and 6)—308 samples (clinical diagnosis: 305 AD and 3 controls). Next, we evaluated the number of gDTUs when comparing individuals at these different stages (Fig. 6). We observed that most gDTUs were detected only while comparing high with either low or mid-Braak stages (Fig. 6A–D). This pattern was observed both in the FL (BM10 and BM44) and TL (BM22 and BM36), suggesting that gene isoform switches positively correlate with AD pathology progression.

Fig. 6: Coincidence between altered expression of splicing-related genes and gDTUs in advanced pathologic stages of AD.
figure 6

A–D Upset plots showing the total number of gDTU identified in the comparison between different Braak stages (low vs. high, low vs. mid, and mid vs. high) in BM10 (A), BM44 (B), BM22 (C), or BM36 (D). Horizontal bars show the total number of gDTUs identified in each comparison (low vs. high, low vs. mid, and mid vs. high), whereas vertical bars indicate the gDTUs exclusive or common to different comparisons. Black dots below vertical bars indicate stages analyzed. E, F Differential expression of genes associated with splicing/spliceosome after comparison of different Braak stages (E) or AD vs controls in different datasets (F). Red and blue squares indicate, respectively, up- and downregulated genes. Gene symbols highlighted in red indicate genes belonging to the neuronal splicing machinery.

Full size image

Next, we set out to evaluate alterations in the expression of genes encoding for proteins of the splicing machinery between the same Braak stages. We found that among 441 genes related to “splicing” or “spliceosome” terms (Supplementary Table 6), 79 were DEGs at high compared to low or mid-Braak stages (Fig. 6E). In contrast, we could not detect any DEG in the comparison of mid vs low Braak stages. Among DEGs detected in the comparison between high and low/mid-Braak stages, we observed that several genes specifically associated with the neuronal splicing regulatory network28, such as RBFOX1 and 2 (RNA binding Fox-1 homolog 1 and 2), ELAVL2 (ELAV like RNA binding protein 2), MBNL3 (muscleblind like splicing regulator 3), PTBP1 (polypyrimidine tract binding protein 1), and NOVA2 (NOVA alternative splicing regulator 2) (Fig. 6E, highlighted in red). A similar correlation between pathological burden and differential expression of the same 441 splice-related genes was observed in the comparison between all AD versus control subjects of the different datasets (Fig. 6F). Changes in the expression of those genes were hardly observed in FL (low number of gDTUs—Fig. 3), but were frequent in TL samples (high number of gDTUs—Fig. 3), albeit to a lesser extent than that observed in the comparison between different Braak stages (likely due to the effects of combining low, mid and high Braak stages in the AD group). Remarkably, the majority of the splicing-related genes with altered expression in the Mayo dataset was not reproduced in the MSBB BM36 dataset, and vice versa (Fig. 6F). This could help to explain the dissimilar consequences of alternative splicing events observed in those datasets (Fig. 5C) and suggest that a myriad of proteins could be involved in altered splicing in the AD brains.

Differential gene expression in separate cell types of the human brain

Considering the cellular diversity in the brain, we took an indirect approach to sort DEGs and gDTUs according to individual cell types. To that, we used scRNAseq data obtained from the adult human brain to identify cell types expressing the genes altered in our DEG/gDTU analysis (Fig. 7 and Supplementary Fig. 1). We found that, out of the 1135 genes with altered expression, i.e., gDTU + DEG, in the TLI (Figs. 2 and 3), 839 were found in at least one cell-type using as cutoff the expression in more than 10% of cells assigned for a specific cell-type (Supplementary Table 7). From these, 239 were identified in unique cell-types/subtypes, 396 in multiple (2–4 cell-types), and 211 in all cell-types analyzed (Fig. 7A, Supplementary Figs. 1 and 2, Supplementary Table 7). Confirming the efficacy of our strategy, GO analyses using cell-type-specific genes revealed that DEGs/gDTUs in the TLI of AD patients were significantly enriched for biological processes associated with inflammation in microglial cells, whereas those associated with cell adhesion were enriched in endothelial cells (Fig. 7B and Supplementary Table 8). Similarly, DEGs/gDTUs identified in neuronal cells were enriched for GO terms such as synaptic signaling, synaptic plasticity, and synapse vesicle cycle (Fig. 7C). Notably, these enrichments were more significant in GABAergic neurons, which could suggest a more pronounced pathological burden on these cells compared to glutamatergic neurons (Fig. 7C). Comparison of the cell-type gene expression signatures identified in our work with previous studies using scRNAseq in AD6,7 showed a similar degree of overlap (Supplementary Figs. 3 and 4; Supplementary Table 9), further supporting the effectiveness of our strategy to assign gene expression alterations to unique cell types in the AD brain.

Fig. 7: Cell-type expression pattern for genes altered in AD brains.
figure 7

A Schematic representation showing our strategy to assign DEGs and gDTUs identified in the TLI to specific cell types of the adult human brain (see also Supplementary Fig. 1). Out of 839 single-cell TLI genes (scTLI), 281 were expressed in a unique cell-type, 249 in 2–4 cell-types, and 77 in all cell-types/subtypes analyzed. B Gene ontology terms enriched for scTLI DEGs, gDTUs, or both per cell type. C Selected GO terms associated with synaptic transmission.

Full size image

DEG/gDTU analyses identify cell-type-specific alterations in AD risk/causal genes

Genomic association studies have revealed about 45 loci containing variants related to an increased or decreased probability of developing AD29,30. However, the functional variants and their target genes remain mostly elusive31. To contribute to the identification of target genes, we first evaluate the expression of 176 genes located within the 45 loci associated with the AD risk (Supplementary Table 10)31 and 3 causal AD genes—PSEN1, PSEN2, and APP—in individual cell types of the adult human brain. We found that 116 out of the 179 AD risk/causal genes were expressed by at least one of the major cell types identified in the brain (Fig. 8A and Supplementary Table 11). Subsets of these genes were exclusively expressed either in microglial cells (14 out of 116), neurons (12), astrocytes (2), oligodendrocytes (6), or endothelial cells (6), suggesting cell-type specific roles for these AD risk/causal genes.

Fig. 8: Expression of AD risk/causal genes is mostly altered in the TL of patients.
figure 8

A Heatmap showing the expression of predicted AD risk/causal genes in different cell types of the adult human brain. DEGs and gDTUs in at least one dataset are highlighted in red. B Venn diagram showing the number of AD risk/causal DEGs or gDTUs identified in the different datasets analyzed. The intersection between Mayo and MSBB TL is highlighted in yellow, and genes identified are shown in the yellow box. C Representation of the 6 most significant BIN1 isoforms altered (left) and quantification of the differential isoform fraction (dIF) in AD brains compared to controls (right). Main protein domains are indicated with different colors. D Similar representation for APP. * dIF >0.05 and FDR > 0.01.

Full size image

Next, we set out to evaluate the differential expression or transcript usage for these genes. Out of the 116 AD risk/causal genes expressed by brain cell types (Fig. 8A), we observed that 54 were also DEGs/gDTUs in at least one of the bulk RNAseq datasets analyzed. Among those genes, two were exclusively identified in the FL (Fig. 8B). We, therefore, decided to focus on the 52 AD risk/causal genes identified in the temporal lobe for further analyses. In this region, we identified 27 and 17 DEGs/gDTUs in the MSBB_TL and Mayo datasets, respectively, including some well-characterized AD risk genes, such as ADAM10 (ADAM metallopeptidase domain 10), BIN1, CLU (Clusterin), and TREM2 (triggering receptor expressed on myeloid cells 2), and the causal AD genes APP, PSEN1 and 2 (presenilin 1 and 2) (Fig. 8A, B). Eight genes were altered in both datasets (Fig. 8B, yellow box; 15,38% of overlap) and were selected for further analysis of isoform switch. Using ISAR to identify the isoforms altered in the AD brains compared to controls, we observed some patterns of isoform switch that could have important functional relevance (Fig. 8C, D). For instance, while BIN1 transcripts ENST00000316724.9 (NP_647593.1—isoform 1) and ENST00000409400.1 (NP_647600.1—isoform 9) were downregulated, transcripts ENST00000393040.7 (NP_647598.1—isoform 6) and ENST00000462958.5, ENST0000046611.5 and ENST00000484253.1 (intron retention) were upregulated (Fig. 8C). This pattern could lead to a decrease of the neuronal-specific BIN1 isoform 1 expression32, given that retained introns are noncoding sequences. Using western blotting analysis, we confirmed this decrease of BIN1 isoform 1 protein in the frontal cortex and hippocampus of AD brain samples compared to controls (Supplementary Fig. 5).

We also observed isoform switches in the AD causal gene APP with possible functional consequences in neuronal cells. While two APP isoforms were downregulated (ENST00000348990 and ENST00000354192), the isoforms ENST00000346798 and ENST00000357903 were upregulated in Mayo and MSBB datasets (Fig. 8D). Noteworthy, significantly downregulated APP isoforms lack exon 7, which contains the Kunitz protease inhibitor (KPI) domain. KPI is one of the main serine protease inhibitors and increased KPI( + )APP mRNA and protein expression levels have been described in AD brains and are associated with increased amyloid-beta deposition33,34,35. At the exception of ENST00000354192, the other transcripts are mostly expressed in neurons (Marques-Coelho and Costa, unpublished data), indicating that these cells may have a selective increase in the expression of KPI(+)APP and, consequently, enhanced production of Aβ1–42.

Discussion

Comprehensive knowledge of gene expression alterations associated with the onset and progression of human diseases is a key step toward the understanding of their cellular and molecular mechanisms36. In this work, we provide a novel framework to identify cell-type-specific gene expression alterations in AD using patient-derived bulk RNAseq. Comparing RNA sequencing data obtained from distinct brain regions of control and AD patients, we show that changes in gene expression are more significant in the temporal than frontal lobe. We also show that a large number of genes present isoform switches without changes in the global expression levels. As a consequence, these genes are overlooked in classical differential expression analysis but can be detected through differential transcript usage analysis. Gene isoform switches are mostly evident at late stages of the pathology and correlate with altered expression of genes encoding for splicing-related proteins. Using an indirect approach to assign genes to unique cell types, we are also able to map DEGs/gDTUs to unique cell populations of the adult brain, and our results are comparable to previously published scRNAseq data6,7. Finally, we show that a subset of AD causal/risk factors such as APP or BIN1 is differentially expressed in the AD brain. Altogether, our work provides a comprehensive description of regional and cell-type-specific gene expression changes in the AD brain and suggests that alternative splicing could be an important mechanism for pathological progression.

Despite the availability of RNAseq datasets generated from healthy subjects and AD patients2,3,4, a systematic evaluation of the gene expression changes in the AD brain, as well as comparisons of these changes in distinct brain regions, was missing. To the best of our knowledge, only one study aimed at comparing gene expression levels in different AD brain regions37, but this work was based on microarray data which has limited gene coverage. We show, using bulk tissue RNAseq data, that alterations in gene expression are highly prominent in biological samples obtained from the temporal lobe, which harbors the first brain regions affected in the AD pathogenesis8. Conversely, few changes are present in biological samples derived from the frontal lobe, where cells are affected only at advanced stages of AD. These observations are in line with recent data showing that changes in protein expression levels in AD brains are much more prominent in the temporal lobe (hippocampus, entorhinal cortex, and cingulate gyrus) than in the frontal lobe (motor cortex)9. They can also help to explain the low number of DEGs identified in scRNAseq data obtained from the frontal lobe7 compared to a similar study in the entorhinal cortex6.

In order to minimize the variability in RNAseq experiments36, we here focused on DEGs (genes with FC > 1.3 and FDR < 0.01 in AD versus control) detected independently in at least two datasets containing samples of similar brain regions (TLI or FLI). These stringent criteria limited the number of DEGs used in subsequent analyses, but still allowed the uncovering of several genes previously associated with AD pathologies, such as ABCA1, ABCA2, CALB1, C1R, C1S, GAD1/2, PVALB, REST, SLC32A1, SST, VGF, and VIP15,16,17,18,19 The reduced number of DEGs in FLI and TLI likely explains our failure to detect functional annotations associated with synaptic transmission and immune response in GSEA, as previously reported13. However, this study analyzed only the ROSMAP dataset and considered genes with FDR < 0.05 as significant, regardless of the fold change, identifying 1722 DEGs in AD versus control brains. Besides the questionable meaning of DEGs with very small fold changes, the use of such a large set of genes for GSEA can artificially increase the number of significantly enriched functional annotations and is not advised20.

Nevertheless, our failure to detect key functional annotations associated with AD pathology while inputting TLI DEGs is puzzling and could suggest that DEG analysis fails to detect relevant alterations in gene expression in the AD brain. Indeed, classical DEG analysis using DESeq or edgeR, which rank all gene transcripts, including noncoding sequences38, are insensitive to the dynamics of gene expression that could, for example, lead to isoform switches with important functional consequences21. Therefore, important gene expression alterations could occur at the level of transcripts, without significant changes in the global expression of genes. According to this possibility, we provide convincing evidence that a high number of genes in the AD brain show isoform switches (DTU) but are not detected by DEG analysis, including several genes associated with the regulation of synapse transmission, such as APP, NSG1, RELN, GABRA122,23,24,25,26. Moreover, gDTUs identified in two independent datasets (TLI), alone or in combination with TLI DEGs, were enriched for key biological processes involved in AD pathogenesis, such as synaptic communication, immune response, inflammation, endocytosis, and cell-signaling39. Similar gene set enrichment has been described using the analysis of co-expression modules in bulk RNAseq27,40 or DEG analysis DEGs in unique cell types in scRNAseq6,7. This could suggest that the combination of DEG and DTU to analyze bulk RNAseq is comparable to scRNAseq regarding the sensitivity to detect gene expression alteration in AD brains. In agreement with this possibility, we were able to assign DEGs and gDTUs to unique cell types and confirm the similarities among cell-type-specific functional annotations observed in our work compared to previous scRNAseq studies6,7.

Notably, we show that several DEGs/gDTUs associated with AD pathogenesis, such as NSG1, CALB1, RELN23,24,25,41 are exclusively assigned to GABAergic neurons. These genes may be particularly relevant for AD pathogenesis, given the central role of GABAergic neurons for the generation of oscillatory rhythms, network synchrony, and memory in different animal models of AD42. Isoform switches in the APP gene could particularly affect GABAergic neurons, which express high levels of that gene, contributing to AD pathogenesis. According to this possibility, conditional knockout of APP/APLP2 only in GABAergic forebrain neurons using DlxCre mice leads to cognitive deficits in hippocampus-dependent spatial learning and memory tasks, associated with altered neuronal morphology and synaptic plasticity43. It is tempting to speculate that GABAergic neurons could be particularly vulnerable in AD, contributing to the increased neuronal activity and synapse downscaling observed in AD brains39,44.

The high number of gDTUs observed in AD brains compared to controls can likely be explained by altered expression of genes encoding for proteins of the splicing machinery, affecting alternative splicing. According to this interpretation, we show that a high number of isoform switches is associated with alternative transcription start site, alternative transcription termination site, exon skipping, alternative 3′ or 5′ splice sites, multiple exon skipping and intron retention. Moreover, we show that several genes encoding for proteins of the splicing machinery have their expression altered in AD brains, especially those showing a high degree of pathology (Braak >5). Also in agreement with the regional differences in gene expression described above, alterations in the splicing machinery are more prominent in the TL than in the FL, which could help to explain the low number of gDTUs in the latter brain region identified in our work and in the previous study using a different strategy to detect isoform switch5.

Particularly interesting, several genes encoding for proteins involved in the control of alternative splicing in neurons are differently expressed in the TL of AD brains. For instance, RBFOX1 and 2 are downregulated in the MSBB BM36 and could contribute to the altered rate of exon skipping observed in this region28,45. Noteworthy, reduced expression of RBFOX1 has been associated with increased inclusion of exon 7 in the APP gene, leading to an enhanced expression of APP isoforms 770 and 751 containing the KPI domain45. A similar switch in the APP isoforms has also been associated with somatic gene recombination in AD46, indicating that increased ratios of APP isoforms containing the KPI domain could be detrimental to neurons. Considering these findings and the well-established associations between KPI(+)APP expression levels, amyloid plaque deposition, and AD pathology progression33,34,35, it is tempting to speculate that controlling APP isoform switches by manipulating RBFOX family proteins could be a potential therapeutic strategy to hamper disease progression.

Altered exon skipping could also help to explain the isoform switch observed for BIN1, which is a major risk factor for AD29,30. BIN1 comprises a N-BAR domain involved in membrane curvature sensing, an SH3 domain that binds to proline-rich motifs, and a clathrin-binding domain (CLAP) specific of the neuronal isoform 132. We show that the transcript encoding for this latter isoform is significantly reduced in the temporal lobe, suggesting that expression of BIN1 isoform 1 in neurons could be reduced. This observation is in line with decreased BIN1 isoform 1 protein expression in the AD brain compared with controls (our own results)47. This would be also in agreement with the observation that overexpression of the BIN1 isoform 1 may be protective in a model of Tauopathy48.

Although we cannot formally rule out that a stage-dependent increase in the number of DEGs and gDTUs could be due to the loss of neuronal cells in brain regions affected by the pathology, several lines of evidence indicate that this is not the most parsimonious explanation for the data described here. First, we observe that the percentage of up- and downregulated genes in GABAergic and glutamatergic neurons are close to 50%, ruling out the possibility that changes in cell numbers could explain these changes. Secondly, previous scRNAseq studies in AD observed a consistent fraction of cell types isolated across control and AD individuals6,7, ruling out significant changes in the cellular composition of AD brains. Lastly, a large number of genes with total expression levels unchanged but presenting isoform switches in the AD brains may likely presuppose a steady cellular composition of the tissue.

Altogether, our work proposes a novel strategy to analyze bulk RNAseq data and identify meaningful gene expression alterations in the diseased brain. It also corroborates previous work implicating alternative splicing in AD susceptibility5 and suggests that isoform switches in the gene BIN1 are involved in the reduced expression of the main neuronal BIN1 isoform 1 in AD brains.

Methods

Bulk RNAseq data from human control and disease banks

RNAseq datasets obtained from different brain regions were used (Mayo2; MSBB3; ROSMAP4). Datasets were downloaded from AMP-AD Knowledge Portal (https://www.synapse.org) following all terms and conditions to use the data. The brain area analyzed and the number of individuals per condition was the following: Mayo - Temporal cortex, which neuroanatomically subdivides into the inferior, middle, and superior temporal gyri (STG), and cytoarchitectonically can be subdivided into Brodmann areas (BM, instead of BA) 20/21/22/41/4249, N = 160 subjects (82 AD and 78 controls); MSBB - BM22, which is part of the Wernicke’s area in the STG, N = 159 subjects (98 AD and 61 controls); MSBB BM36, corresponding to the lateral perirhinal cortex, N = 154 subjects (88 AD and 64 controls); MSBB BM10, corresponding to the anterior prefrontal cortex, N = 176 subjects (105 AD and 71 controls); MSBB BM44, corresponding to the inferior frontal gyrus, N = 153 subjects (90 AD and 63 controls); and ROSMAP - Dorsolateral prefrontal cortex (DLPFC), containing BM46 and part of BM9, N = 423 subjects (222 AD and 201 controls). Unless stated otherwise, data obtained from different analyses were grouped in “temporal lobe” (TL) - Mayo, MSBB BM22 and MSBB BM26; or “frontal lobe” (FL) - ROSMAP, MSBB BM10 and MSBB BM44.

Metadata obtained from each study was used to classify patients into Control and Alzheimer’s disease groups (Supplementary Table 12). Briefly, for the MSBB dataset, we used patients with Neuropathology Category (NP.1) labeled as “Control” and “definitive Alzheimer”. For the Mayo dataset, we used the “Diagnosis” column of the metadata, selecting only “AD” and “Control” patients. For the ROSMAP dataset, we also used the column “Diagnosis” of the metadata, selecting only “Control” (value = 1) and “Alzheimer with no other conditions” (value = 4). In all those datasets, subjects marked as “AD” showed Braak stage values higher than 4. In the MSBB dataset, CDR scores of AD patients were consistently higher than 2. In the Mayo and ROSMAP datasets, all AD patients had also a definitive diagnosis according to NINCDS criteria. Covariates such as “Postmortem interval (PMI)”, “RNA integrity number (RIN)”, “Age of death”, and “Sex” were balanced among the different groups (Table 1; Chi-square test, p > 0.05). We used RIN and PMI as covariates in our model to control for possible “batch effects” (linear regression). For detailed information of all individual samples used in this study, please refer to Supplementary Table 13.

Table 1 Summary of clinical, demographic, and technical variables of samples analyzed from different datasets.
Full size table

For single-cell RNA sequencing (scRNAseq) analyses, we used processed data obtained from the middle temporal gyrus (MTG), available at the Allen Brain Atlas consortium (https://celltypes.brain-map.org/rnaseq).

Realignment of human reads into single pseudoaligner pipeline

Using human GRCh38 cDNA release 94 (ftp://ftp.ensembl.org/pub/release-94) as a reference, we built an index to align all our fastq files. Next, we used pseudoaligner Kallisto50 (version 0.43.1) with our pre-built index to align fastq files.

Differential gene expression analyses

Differentially expressed genes (DEGs) were identified using differential gene expression at transcript-level using DESeq2 R library51,52. To facilitate kallisto output import, transcript-level estimated counts, length, and abundance were extracted using tximport function53. As described by Michael Love group, transcript-level differential gene expression enhances analysis resolution52. Using DESeqDataSetFromTximport, a DESeq2 object was created and filtered using rows with sum of all counts bigger than 10. Next, DESeq function was used with default values. Using the results function, we selected all genes with a false discovery rate (FDR) < 0.01 and fold change (FC) > 1.3. We also used RIN and PMI as covariates (linear regression).

Differential transcript usage (DTU) analysis was performed using the R library IsoformSwitchAnalyzeR12. Following pipeline instructions, kallisto abundance tables were imported using importIsoformExpression and importRdata functions to create a switchAnalyzeRlist object. Same cDNA release used in kallisto alignment and correspondent annotation (ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.chr_patch_hapl_scaff.gtf.gz) were applied as input. We filtered data using a gene expression cutoff = 10, isoform expression cutoff = 3, differential isoform fraction (dIF) cutoff = 0.05 and removed single isoform genes. Although DEXSeq is recommended to test differential isoform usage, it does not work efficiently for large datasets (more than 100 samples)11. For that reason, we chose isoformSwitchAnalysisPart1 function using DRIMSeq53 to test differential transcript usage. Using part1 fasta files, all external analysis was performed and used as input to isoformSwitchAnalysisPart2 function. We used CPC2, Pfam, SignalIP and Netsurfp2 as indicated in the pipeline. Next, we performed a confirmation stage using stageR14 to generate isoforms overall false discovery rate (OFDR). We selected all isoforms with OFDR < 0.01 and dIF >0.05. RIN and PMI metadata were used as covariates (linear regression).

Statistical significance of the intersections among different datasets was calculated using the hypergeometric test (phyper).

Splicing events and event consequences

We used extractSplicingSummary and extractConsequenceSummary functions to quantify gain/loss of predicted splicing events (such as exon skipping and intron retention); and gain/loss (also sensitive/insensitive, shorter/longer and switch) of predicted functional consequences (such as coding potential and domain identified), respectively.

Single-cell RNAseq

Using R library seurat54, we created a seurat object (CreateSeuratObject), normalized data (NormalizeData), found variable genes (FindVariableFeatures), and rescaled data using a linear model (ScaleData, use.umi = F). After that, we generated 50 PC’s (RunPCA) but only used 35 of them based on the PC’s visualization distribution (ElbowPlot). Since Allen data were already annotated, we only used tSNE (RunTSNE) to facilitate visualization. A group classified as “None” by Allen metadata were removed from our analysis. This strategy generated 7 main different cell types: Astrocytes, Endothelial cells, Glutamatergic Neurons, GABAergic Neurons, Microglia, Oligodendrocytes and oligodendrocyte precursor cells (OPCs). To assign genes to specific cell types, we used the AverageExpression function. Using pct.exp bigger than 0.1, we created a list of genes that were expressed by each cell type.

Gene set enrichment analysis (GSEA)

For gene ontology analysis, R library gprofiler255 was used. Using gost function, correction_method = “fdr” and significant=TRUE. To minimize the enrichment of gene ontologies based on a small set of genes, we used three conditions for significance assessment: false discovery rate (FDR) < 0.01; intersection size (intersection between gene set vs. a number of genes in a term) >3; and precision (intersection size divided by gene set) >0.03. We used Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC), Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome (REAC), WikiPathways (WP), TRANSFAC (TF), miRTarBase (MIRNA), Human Protein Atlas (HPA), CORUM (CORUM), Human phenotype ontology (HP) as sources. For improved visualization, we plotted results only for GO:BP, GO:CC and KEGG and show only FDR related to terms reaching all criteria of significance.

Selection of splicing-associated genes

To select splicing-related genes, we searched for terms containing the words “splicing” or “spliceosome” in gProfiler bank (https://biit.cs.ut.ee/gprofiler/gost). Taking only GO and WP datasets, 25 terms and 441 genes related to those terms were selected (Supplementary Table 6).

Selection of AD risk/causal genes

The complete list of AD risk/causal genes used in this study is described in Supplementary Table 10. Briefly, AD risk loci were selected from previous work using genome-wide association studies and whole exome sequencing29,30. AD risk genes within these loci were determined based on regional association plots, assuming that the functional risk variants are located in the vicinity of the SNP producing the top signal and taking into account the linkage disequilibrium patterns and the recombination peaks within the loci of interest31. Early-onset AD causal genes used in this study are APP, PSEN1, and PSEN2.

Western blotting

Frozen brain samples obtained from the frontal cortex (FCx) and hippocampus (hip) of three non-pathology (age: 80.33 ± 3.78 years; Braak: 2.66 ± 1.15; PMI: 37.33 ± 22.50 h) and six AD patients (age: 79.57 ± 6.70 years; Braak: 6; PMI: 26.57 ± 13.40 h) were lysed with RIPA buffer and sonicated at 100% during 10 s before use for the western blotting analyses. The controls for BIN1 isoforms 1 (Iso1) and 9 (Iso9) were obtained using HEK cells transiently transfected with 1 µg/ml DNA solution containing plasmids encoding for BIN1 isoforms mixed with the transfection reagent FuGENE HD (Promega) at the ratio 1:3. Cells were lysed using RIPA buffer 48 h after transfection and frozen for further analyses.

Protein quantification was performed using the BCA protein assay (Thermo Scientific). In total, 10–20 μg of total protein from extracts were separated in SDS–polyacrylamide gels 4–12% (NuPAGE Bis-Tris, Thermo Scientific) and transferred to nitrocellulose membranes (Bio-Rad). Next, membranes were incubated in milk (5% in Tris-buffered saline with 0.1% Tween-20, TTBS) to block non-specific binding sites during 1 h at RT, followed by several washes with TTBS. Immunoblotting was carried out with primary antibodies anti-BIN1 (Abcam, ab182562), anti-β-ACTIN (Sigma-Aldrich, A1978), and anti-GAPDH (Millipore, AB2302) overnight at 4 °C on 20 RPM. The membranes were washed three times in TTBS, followed by incubation with HRP-conjugated secondary antibodies (Jackson, anti-mouse 115-035-003 and anti-rabbit 111-035-003; Thermo Scientific, anti-chicken A16054) overnight at 4 °C on 20 RPM agitation. The membranes were washed three times in TTBS, and the immunoreactivity was revealed using the ECL chemiluminescence system (SuperSignal, Thermo Scientific) and imaged using the Amersham Imager 600 (GE Life Sciences). Optical densities of bands were quantified using “Gel Analyzer” plugin in Fiji56.

All western blot experiments were performed in compliance with relevant guidelines provided by the Neuro-CEB - Biological Resources Platform and the protocols were approved by the National Institute of Health and Medical Research (INSERM), at the Institut Pasteur de Lille, University of Lille.

Reporting summary

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