Abstract
Alzheimer’s disease (AD) is the leading cause of dementia in aging individuals. However pathophysiological processes involved in the brain are still poorly understood. Among numerous strategies, a comprehensive overview of gene expression alterations in the diseased brain has been proposed to help for a better understanding of the disease processes. In this work, we probed the differential expression of genes in different brain regions of healthy and AD adult subjects using data from three large studies: Mayo Clinic; Mount Sinai Brain Bank (MSBB) and ROSMAP. Using a combination of differential expression of gene and isoform switch analyses we provide a detailed landscape of gene expression alterations in the temporal and frontal lobes, harboring brain areas affected at early and late stages of the AD pathology, respectively. Next, we took advantage of an indirect approach to assign the complex gene expression changes revealed in bulk RNAseq to individual cell types of the adult brain. This strategy allowed us to identify cell type/subtype specific isoform switches in AD brains previously overlooked. Among these alterations, we show isoform switches in the AD causal gene APP (Amyloid Beta Precursor Protein) and the risk gene BIN1 (Bridging Integrator 1), which could have important functional consequences in neuronal cells. Altogether, our work proposes a novel integrative strategy to analyze RNAseq data in AD based on both gene/transcript expression and regional/cell-type specificities.
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 processes (Masters et al., 2015). 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 individuals (Allen et al., 2016; De Jager et al., 2018; Wang et al., 2018). 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) (Raj et al., 2018) or single-cells RNA sequencing (scRNAseq) (Grubman et al., 2019; Mathys et al., 2019). However, these studies have focused on samples obtained from different brain regions, namely the dorsolateral prefrontal (Mathys et al., 2019; Raj et al., 2018) and entorhinal cortices (Grubman et al., 2019), 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) (Braak & Braak, 1991). 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 cerebellum (Xu et al., 2019).
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 transcripts (Yi et al., 2018). Recently, new approaches using transcripts-level analysis, so called differential transcript usage (DTU), enables identification of alternative splicing and isoform switches with prediction of functional consequences (Anders et al., 2012; Vitting-Seerup & Sandelin, 2019). 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 the Mayo’s clinic study, both the temporal cortex and cerebellum were used to obtain bulk RNAseq (Allen et al., 2016). In the Religious Orders Study (ROS) and Memory and Aging Project (MAP), henceforth called ROSMAP, the dorsolateral prefrontal cortex was used (De Jager et al., 2018). Finally, in the Mount Sinai/JJ Peters VA Medical Center Brain Bank (MSBB), 4 different Brodmann areas of the brain were studied: areas 22 and 36 from the temporal lobe, areas 10 and 44 in the frontal lobe (Wang et al., 2018). 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 AD (Allen et al, 2016; Wang et al, 2018; De Jager et al, 2018). Considering the regional progression of AD pathology (Braak and Braak, 1991), 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 (Figure 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 3,348 (1244 down- and 2104 up-regulated genes) and 2,172 (1170 down and 999 up-regulated genes in BM22 and BM36; 3 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 (Figure 2A-B; 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) (Figure 2A-B; 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 up-regulated genes with FC>1.3; Canchi et al., 2019). Among DEGs detected in the FL, 62.5% were also detected in the TL (Figure 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 (Figure 1) and differences are observed even in samples obtained from the same donors (compare MSBB_TL and MSBB_FL in Figure 2). Thus, changes in gene expression are much more prominent in brain areas affected at early stages of AD pathology.
Datasets obtained from 3 consortia (Mayo, MSBB and ROSMAP) were grouped according to the brain region sampled in frontal lobe (FL) or temporal lobe (TL). Next, RNAseq data was pseudo-aligned using Kallisto and analyzed using the packages from R (version 3.6) DESeq2, IsoformSwitchAnalyzeR (ISAR) and gene-set enrichment analysis (GSEA). Assignment of differentially expressed genes or isoform switches to specific cell types was performed indirectly using scRNAseq signatures obtained from the Allen Brain Atlas transcriptomic data.
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).
To select genes consistently altered in AD brains, considering the several sources of measurement variations in RNAseq experiments (Van den Berge et al., 2017), 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 diseases (Abuznait & Kaddoumi, 2012); 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 accumulate (Benoit et al., 2013); REST (RE1 Silencing Transcription Factor), which regulates neural circuit activity during aging (Zullo et al., 2019); GAD1 and 2 (Glutamate Decarboxylase 1 and 2), SLC32A1 (Solute Carrier Family 32 GABA Vesicular Transporter, Member 1), CALB1 (Calbindin 1), PVALB (Parvalbumin), SST (Somatostatin) and VIP (Vasoactive Intestinal Peptide), all expressed in GABAergic neurons and involved in cognitive decline in AD and other neurological diseases (Prévot & Sibille, 2020). 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 disease (Beckmann et al., 2020).
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 analysis (Subramanian et al., 2005). 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 (Figure 2C; 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. (2019). 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 splicing (Vitting-Seerup & Sandelin, 2017; Yi et al., 2018). 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 2,509 and 1,843 genes with differential transcript usage (gDTU) in the temporal lobe of AD brains studied in the Mayo and MSBB datasets, respectively (Figure 3A-B; 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, Figure 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 (Figure 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 (Figure 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 types (Mathys et al., 2019; Grubman et al., 2019) 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.
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.
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.
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 AD (Kamenetz et al., 2003); NSG1 (Neuronal Vesicle Trafficking Associated 1), which has been implicated in the regulation of AMPA receptors (AMPAR) and APP trafficking, thus affecting synaptic transmission, plasticity and Aβ production (Alberi et al., 2005; Norstrom et al., 2010); RELN (Reelin) that plays important role in synaptic transmission and has been associated with AD (Yu et al., 2016); GABRA1 (Gamma-Aminobutyric Acid Type A Receptor Subunit Alpha1), 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 brain (Limon et al., 2012).
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 (Figure 5; 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 (Figure 5B). Other common splicing events observed were alternative 3’ or 5’ splice sites (A3 and A5, respectively), multiple exon skipping (MES) and intron retention (IR) (Figure 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 brain (Raj et al., 2018).
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.
Alternative splicing events may have diverse functional consequences for protein expression, such as shifting the frequency of transcripts containing introns (non-coding) or mRNA stability (nonsense mediated decay) or leading to gain/loss of protein domains, intrinsically disordered regions or signaling peptides (Vitting-Seerup & Sandelin, 2019). Quantification of these consequences revealed some interesting differences between Mayo and MSBB BM36 (Figure 5C), the two datasets with 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 (Figura 5B). These differences could be at least partly explained by the larger number of gDTUs detected in the Mayo compared to MSBB TL (Figure 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; Wan et al., 2020).
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 (Figure 6; Supplementary Table 5). For this purpose, we subdivided samples in 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 (Figure 6). We observed that most gDTUs were detected only while comparing high with either low or mid Braak stages (Figure 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.
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.
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 (Figure 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 network (Raj & Blencowe, 2015), such as RBFOX 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) (Figure 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 (Figure 6F). Changes in the expression of those genes were hardly observed in FL (low number of gDTUs – Figure 3), but were frequent in TL samples (high number of gDTUs – Figure 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 (Figure 6F). This could help to explain the dissimilar consequences of alternative splicing events observed in those datasets (Figure 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 (Figure 7; Supplementary Figure 1). We found that, out of the 1135 genes with altered expression, i.e. gDTU + DEG, in the TLI (Figures 2 and 3), 839 were found in at least one cell-type using as cut-off 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 (Figure 7A; Supplementary Figures 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 (Figure 7B; 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 (Figure 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 (Figure 7C). Comparison of the cell-type gene expression signatures identified in our work with previous studies using scRNAseq in AD (Mathys et al., 2019; Grubman et al., 2019) showed a similar degree of overlap (Supplementary Figures 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.
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 Figure 1). Out of 839 single-cell TLI genes (scTLI), 281 were expressed in a unique cell-type, 249 in 2 to 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.
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 AD (Kunkle et al., 2019; Lambert et al., 2013). However, the functional variants and their target genes remain mostly elusive (Dourlen et al., 2019). To contribute with the identification of target genes, we first evaluate the expression of 176 genes located within the 45 loci associated with AD AD risk (Supplementary table 10; Dourlen et al., 2019) 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 (Figure 8A; 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.
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.
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 (Figure 8A), we observed that 54 were also DEGs/gDTUs in at least one of the bulk RNAseq datasets analyzed. Among those genes, 2 were exclusively identified in the FL (Figure 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) (Figure 8A-B). Eight genes were altered in both datasets (Figure 6B, 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 (Figure 8C and 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 (Figure 8C). This pattern could lead to a decrease of the neuronal specific BIN1 isoform 1 expression (Zhou et al., 2014), given that retained introns are non-coding 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 Figure 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 (Figure 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 deposition (Tanzi et al., 1989; Johnson et al., 1989; Kitaguchi al., 1988). 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
A comprehensive knowledge of gene expression alterations associated with the onset and progression of human diseases is a key step towards the understanding of their cellular and molecular mechanisms (Lee & Young, 2013). 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 data (Grubman et al, 2019; Mathys et al, 2019). Finally, we show that a subset of AD causal/risk factors such as APP or BIN1 are 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 patients (Allen et al, 2016; Wang et al, 2018; De Jager et al, 2018), a systematic evaluation of the gene expression changes in the AD brain, as well as comparisons of these changes in distinct brain region, was missing. To the best of our knowledge, only one study aimed at comparing gene expression levels in different AD brain regions (Haroutunian et al., 2009), but this work was based on microarray data which has a limited gene coverage. We here 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 pathogenesis (Braak and Braak, 1991). Conversely, few changes are present in biological samples derived from the frontal lobe, where cells are affected only at advanced stages of the 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) (Xu et al., 2019). They can also help to explain the low number of DEGs identified in scRNAseq data obtained from the frontal lobe (Mathys et al., 2019) compared to a similar study in the entorhinal cortex (Grubman et al., 2019).
In order to minimize the variability in RNAseq experiments (Van den Berge et al., 2017) 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 pathology, such as ABCA1, ABCA2, CALB1, C1R, C1S, GAD1/2, PVALB, REST, SLC32A1, SST, VGF and VIP (Abuznait & Kaddoumi, 2012; Beckmann et al., 2020; Benoit et al., 2013; Prévot & Sibille, 2020; Zullo et al., 2019). 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 reported (Canchi et al., 2019). However, this study analyzed only the ROSMAP dataset and considered genes with FDR<0.05 as significant, regardless of the fold change, identifying 1,722 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 advised (Subramanian et al., 2005).
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 non-coding sequences (Costa-Silva et al., 2017), are insensitive to the dynamics of gene expression that could, for example, lead to isoform switches with important functional consequences (Vitting-Seerup and Sandelin, 2017). 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, GABRA1 (Alberi et al., 2005; Kamenetz et al., 2003; Limon et al., 2012; Norstrom et al., 2010; Yu et al., 2016). 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-signaling (Canter et al., 2016). Similar gene set enrichment has been described using the analysis of co-expression modules in bulk RNAseq (Milind et al., 2020; Wan et al., 2020) or DEG analysis DEGs in unique cell types in scRNAseq (Mathys et al., 2018; Grubman et al., 2018). 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 studies (Mathys et al., 2019; Grubman et al., 2019).
Notably, we show that several DEGs/gDTUs associated with AD pathogenesis, such as NSG1, CALB1, RELN (Alberi et al., 2005; Norstrom et al., 2010; Yu et al., 2015; Odero et al., 2010) are exclusively assigned to GABAergic neurons. These genes may be particularly relevant for AD pathogenesis, given the central role of GABAergic neurons for generation of oscillatory rhythms, network synchrony, and memory in different animal models of AD (Verret et al., 2012; Etter et al., 2019). 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 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 plasticity (Mehr et al., 2020). 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 brains (Canter et al., 2016; Dörrbaum et al., 2020).
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 previous study using a different strategy to detect isoform switch (Raj et al., 2018).
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 down regulated in the MSBB BM36 and could contribute to the altered rate of exon skipping observed in this region (Alam et al., 2014; Raj & Blencowe, 2015). Noteworthy, reduced expression of RBFOX1 has been associated with an increased inclusion of exon 7 in the APP gene, leading to an enhanced expression of APP isoforms 770 and 751 containing the KPI domain (Alam et al., 2014). A similar switch in the APP isoforms has also been associated with somatic gene recombination in AD (M. H. Lee et al., 2018), 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 APP-KPI expression levels, amyloid plaque deposition and AD pathology progression (Tanzi et al., 1988, Johnson et al. 1990; Kitaguchi et al., 1988), 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 AD (Lambert et al., 2013; Kunkle et al., 2019). 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 1 (Zhou et al., 2014). 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 is reduced. This observation is in line with decreased BIN1 isoform 1 protein expression in the AD brain compared with controls (our own results; Glennon et al., 2013). This would be also in agreement with the observation that an overexpression of the BIN1 isoform 1 may be protective in a model of Tauopathy (Sartori et al., 2019).
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. Firstly, we observe that the percentage of up and down regulated 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 individuals (Mathys et al., 2019; Grubman et al., 2019), ruling out significant changes in cellular composition of AD brains. Lastly, the 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 susceptibility (Raj et al., 2018) 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 (Mayo: Allen et al, 2016; MSBB: Wang et al, 2018; ROSMAP: De Jager et al, 2018). 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 were 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/42 (Strotzer, 2009), 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. Co-variates such as “Post-mortem 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.
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 pseudo aligner pipeline
Using human GRCh38 cDNA release 94 (ftp://ftp.ensembl.org/pub/release-94) as reference, we built an index to align all our fastq files. Next, we used pseudoaligner Kallisto (Bray et al., 2016; 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 library (Love et al., 2014, 2019). To facilitate kallisto output import, transcript-level estimated counts, length and abundance was extracted using tximport function (Soneson et al., 2016). As described by Michael Love group, transcript-level differential gene expression enhances analysis resolution (Love et al., 2019). 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 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 IsoformSwitchAnalyzeR (Vitting-Seerup & Sandelin, 2019). 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 cut off = 10, isoform expression cut off = 3, differential isoform fraction (dIF) cut off = 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; (Anders et al., 2012)). For that reason, we chose isoformSwitchAnalysisPart1 function using DRIMSeq (Soneson et al., 2016) to test differential transcript usage. Using part1 fasta files, all external analysis was performed and used as input to soformSwitchAnalysisPart2 function. We used CPC2, Pfam, SignalIP and Netsurfp2 as indicated in the pipeline. Next, we performed a confirmation stage using stageR (Van den Berge et al., 2017) 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 seurat, 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 was 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 gprofiler2 was used. Using gost function, correction_method=“fdr’’ and significant=TRUE. To minimize the enrichment of gene ontologies based on small set of genes, we used three conditions for significance assessment: false discovery rate (FDR)<0.01; intersection size (intersection between gene set vs. 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 for 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 the Supplementary Table 10. Briefly, AD risk loci were selected from previous work using Genome-wide Association Studies and Whole Exome Sequencing (Lambert et al., 2013; Kunkle et al., 2019). 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 interest (Dourlen et al., 2019). 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 3 non-pathology (Age: 80.33 ± 3.78 years; Braak: 2.66 ± 1.15; PMI: 37.33 ± 22.50 hours) and 6 AD patients (Age: 79.57 ± 6.70 years; Braak: 6; PMI: 26.57 ± 13.40 hours) were lysed with RIPA buffer and sonicated at 100% during 10 seconds 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 48h after transfection and frozen for further analyses.
Protein quantification was performed using the BCA protein assay (Thermo Scientific). 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; 1h at RT) to block non-specific binding sites during 30min 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) for 1h at RT on 20 RPM. The membranes were washed 3 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. 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 Fiji (Schindelin et al., 2012).
Data Availability
All data referred to in the manuscript is publicly available
https://doi.org/10.7303/syn3388564
Data availability
MSBB RNAseq (De Jager et al, 2018) - https://doi.org/10.7303/syn3388564
ROSMAP RNAseq (Wang et al, 2018) - https://doi.org/10.7303/syn3157743
MAYO RNAseq (Allen et al, 2018) - https://doi.org/10.7303/syn5550404
Code availability
Codes used for analyses are available, differentially expressed genes and differential transcript usage in individual datasets can be assessed using https://diegomscoelho.github.io/AD-IsoformSwitch/
Author Contributions
M.R.C. conceived of the presented idea and supervised the experiments. D.M-C. performed all bioinformatic analyses and prepared the figures and tables of the manuscript. L.I.C.C. contributed with bioinformatic analyses. A.R.M.F. performed western blot analysis in human samples. D.M-C., J-C. L. and M.R.C. wrote the manuscript. All authors revised and approved the final version of the manuscript.
Competing Interests statement
The authors declare no competing financial and non-financial interests.
Acknowledgments
We thank Charles Duyckaerts (Neuroceb, GIE Neuro-CEB BB-0033-00011) for providing the brain samples. We would also like to thank Dr. Tarciso Velho, Pierre Dourlen and Fabien Delahaye for reading and suggestions on the manuscript. The results published here are in whole or in part based on data obtained from the AMP-AD Knowledge Portal (https://adknowledgeportal.synapse.org/). Study data were provided by the following sources: The Mayo Clinic Alzheimer’s Disease Genetic Studies, led by Dr. Nilufer Taner and Dr. Steven G. Younkin, Mayo Clinic, Jacksonville, FL using samples from the Mayo Clinic Study of Aging, the Mayo Clinic Alzheimer’s Disease Research Center, and the Mayo Clinic Brain Bank. Data collection was supported through funding by NIA grants P50 AG016574, R01 AG032990, U01 AG046139, R01 AG018023, U01 AG006576, U01 AG006786, R01 AG025711, R01 AG017216, R01 AG003949, NINDS grant R01 NS080820, CurePSP Foundation, and support from Mayo Foundation. Study data includes samples collected through the Sun Health Research Institute Brain and Body Donation Program of Sun City, Arizona. The Brain and Body Donation Program is supported by the National Institute of Neurological Disorders and Stroke (U24 NS072026 National Brain and Tissue Resource for Parkinson’s Disease and Related Disorders), the National Institute on Aging (P30 AG19610 Arizona Alzheimer’s Disease Core Center), the Arizona Department of Health Services (contract 211002, Arizona Alzheimer’s Research Center), the Arizona Biomedical Research Commission (contracts 4001, 0011, 05-901 and 1001 to the Arizona Parkinson’s Disease Consortium) and the Michael J. Fox Foundation for Parkinson’s Research. Study data were also provided by the Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago. Data collection was supported through funding by NIA grants P30AG10161 (ROS), R01AG15819 (ROSMAP; genomics and RNAseq), R01AG17917 (MAP), the Illinois Department of Public Health (ROSMAP), and the Translational Genomics Research Institute (genomic). For ROSMAP, additional phenotypic data can be requested at www.radc.rush.edu. And part based on data generated from postmortem brain tissue collected through the Mount Sinai VA Medical Center Brain Bank and were provided by Dr. Eric Schadt from Mount Sinai School of Medicine.
This work was co-funded by the European Union under the European Regional Development Fund (ERDF) and by the Hauts de France Regional Council (contrat n°18006176), the MEL (contract_2016_ESR_05), and the French State (contract n°2018-3-CTRL_IPL_Phase2) to MRC. This work was also funded by the Lille Métropole Communauté Urbaine and the French government’s LABEX DISTALZ program (Development of innovative strategies for a transdisciplinary approach to Alzheimer’s disease) to JCL. DMC and LICC are supported by Coordenaç ão de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) scholarships.
Footnotes
In this revised version, we have included new analysis describing the splicing events responsible for the isoform switches observed in the AD brain, as well as the consequences of these splicing events for protein expression. We have also added detailed clinico, pathological and demographic information about the samples used and improved the quality of figures.
References
