Ontology Specific Alternative Splicing Changes in Alzheimer's Disease.

Alternative splicing (AS) is a common phenomenon and correlates with aging and aging-related disorders including Alzheimer's disease (AD). We aimed to systematically characterize AS changes in the cerebral cortex of 9-month-old APP/PS1 mice. The GSE132177 dataset was downloaded from GEO and ENA databases, aligned to the GRCm39 reference genome from ENSEMBL via STAR. Alternative 3'SS (A3SS), alternative 5'SS (A5SS), skipped exon (SE), retained intron (RI), and mutually exclusive exons (MXE) AS events were evaluated using rMATS, rmats2sashimiplot, and maser. Differential genes or transcripts were analyzed using the limma R package. Gene ontology analysis was performed with the clusterProfiler R package. A total of 60,705 raw counts of AS were identified, and 113 significant AS events were finally selected in accordance with the selection criteria: 1) average coverage >10 and 2) delta percent spliced in (ΔPSI) >0.1. SE was the most abundant AS event (61.95%), and RI was the second most abundant AS type (13.27%), followed by A3SS (12.39%), thereafter A5SS and MXE comprised of 12.39%. Interestingly, genes that experienced SE were enriched in histone acetyltransferase (HAT) complex, while genes spliced by RI were enriched in autophagy and those which experienced A3SS were enriched in methyltransferase activity revealed by GO analysis. In conclusion, we revealed ontology specific AS changes in AD. Our analysis provides novel pathological mechanisms of AD.

Introduction

More than 50,000 genes and 140,000 transcripts for humans have been documented in the Ensemble database, indicating an average of approximately three transcripts or isoforms for each gene caused by alternative splicing (AS). AS is mediated by the spliceosome, a large macromolecular complex that coordinates and catalyzes splicing reactions (Wilkinson et al., 2020). Introns are defined by 5′ splice site (5′SS), branch point (BP), and 3′ splice site (3′SS). The small nuclear ribonucleoproteins (snRNP) U1 and U2 recognize 5′SS and BP sequence, respectively, and the intron is released in lariat form and degraded after two transesterification reactions. Finally, the adjacent two exons are joined by a phosphodiester bond (Leung et al., 2011). Generally, five basic modes of AS are classified: alternative 3′SS (A3SS), alternative 5′SS (A5SS), skipped exon (SE), retained intron (RI), and mutually exclusive exons (MXE) (Sammeth et al., 2008). AS process orchestrates the temporal and tissue-specific expression of numerous genes during development and aging (Baralle and Giudice, 2017; Bhadra et al., 2020).

Alzheimer’s disease (AD) is the top cause of dementia and is characterized by a progressive decline of cognitive functions including memory, language, and visuospatial skills (Querfurth and LaFerla, 2010). The neuropathology of AD includes extracellular β-amyloid deposition and neurofibrillary tangles formed by hyperphosphorylated tau (Ferreira and Klein, 2011). Age, environmental and genetic risk factors also contribute to AD pathogenesis. Genetic risk factors involve amyloid precursor protein, presenilin 1 (PSEN1), sortilin-related receptor (SORL1), and triggering receptor expressed on myeloid cells 2 (TREM2) (Kulkarni et al., 2021).

Accumulating transcriptomic studies revealed AS disruption in AD (Yang et al., 2021). Variants within the fourth exon of PSEN1 cause exon 4 skipping and produce a loss of function PSEN1, leading to an early onset of AD (Tysoe et al., 1998). Monti G et al observed an inclusion of a novel exon in SORL1 gene that encoded a truncated protein. Moreover, this novel transcript was mainly located in neuronal dendrites and was decreased in AD patients (Monti et al., 2021). Han S et al observed that the second exon of TREM2 was more frequently skipped in individuals having at least one low-frequency TREM2 variant, leading to an enrichment of immune-related functional pathways revealed by gene ontology (GO) analysis of differently expressed genes (DEGs) (Han et al., 2021). SE and RI events have been widely observed in AD (Li et al., 2021; Yang et al., 2021)).

In this study, we systematically characterized AS changes in the cerebral cortex of 9-month APP/PS1 mice. We revealed ontology specific AS changes in AD.

Materials and Methods

Dataset and Reference Genome

GSE132177 was deposited in GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132177). The raw fastq files were downloaded from the ENA database (https://www.ebi.ac.uk/ena/browser/view/PRJNA546262). This dataset performed RNA sequencing of the cerebral cortex from 9-month APP/PS1 and control mice. The reference genome was obtained from the ENSEMBL database. The gtf (http://ftp.ensembl.org/pub/release-105/gtf/mus_musculus/Mus_musculus.GRCm39.105.gtf.gz), DNA fasta (http://ftp.ensembl.org/pub/release-105/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz), and cDNA fasta (http://ftp.ensembl.org/pub/release-105/fasta/mus_musculus/cdna/Mus_musculus.GRCm39.cdna.all.fa.gz) files were downloaded for reading alignment.

Reads Alignment

The raw fastq files were downloaded from the ENA database using Axel (v2.17.5), and the md5 numbers were checked using md5sum (v8.30). Quality control was performed using FastQC (v0.11.9) and multiqc (v1.11). Adaptors and low-quality bases were removed using trim-galore (v0.6.7) with the following parameters: trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3. The clean data were aligned to the reference genome using STAR (v 2.7.10a) (Dobin et al., 2013).

Gene/Transcript Expression Matrix

The aligned bam files after STAR alignment were further used to generate an expression matrix on gene level with featureCounts (v2.0.1) (Liao et al., 2014). Clean fastq files from trim-galore were further used to generate an expression matrix on transcript level with salmon (v1.0.4) (Patro et al., 2017).

Differential Expressed Gene/Transcript Analysis

The gene expression matrix and transcript expression matrix were used for differentially expressed gene/transcript analysis using limma R package (v3.50.0) (Ritchie et al., 2015). The cutoff p value was 0.5, and the cutoff logFC value was 1.5.

Alternative Splicing Analysis

The aligned and sorted bam files after STAR alignment was used for AS analysis using rMATS (v4.1.2) and rmats2sashimiplot (v2.0.4) (Shen et al., 2014). A3SS, A5SS, SE, RI, and MXE events were evaluated. Significant AS events were identified with the following criteria (Wilkinson et al., 2020): average coverage >10 (Leung et al., 2011); delta percent spliced in (ΔPSI) >0.1. The AS event statistics were performed using maser (v1.12.1) and rtracklayer (v1.54.0) R packages, and the results were plotted using echarts4r (v0.4.3), maser, and ggplot2 (v3.3.5) R packages.

Gene Ontology Enrichment Analysis

Gene ontology (GO) enrichment analysis was used to perform enrichment analysis of differentially spliced gene sets with clusterProfiler (v4.2.2) (Wu et al., 2021), and GOplot (v1.0.2) R packages (Walter et al., 2015). Significant GO terms were collected if the q value < 0.05.

Results

Global Summary of Alternative Splicing Events

A total of 60,705 AS events were initially detected in the GSE132177 dataset (Figure 1A). Eventually, 113 significant AS events were collected according to our strict criteria (average coverage >10, ΔPSI >0.1) (Figure 1B; Table 1). The principal component analysis (PCA) of these AS events could clearly classify AD from control mice (Figure 1C). Increased PSI was found in A3SS, A5SS, and RI in the AD group, although ΔPSI showed similar changes globally (Figure 1D,E). The top three events in each type of AS pattern were included in Table 2, and all 113 significant AS events were shown in Supplementary Table S1.

FIGURE 1

Global summary of alternative splicing events. (A) Summary of the raw counts in each AS event. (B) The number and ratio of significant AS events. (C) The PCA plot of the PSI values. AD and control groups were clearly classified. (D) The percent of ΔPSI in AS events. (E) The distribution of ΔPSI in AS events.

TABLE 1

Event summary.

Event type	Total events junction counts	Average coverage >10	FDR<0.05, deltaPSI >0.1	Ratio in WT (%)	Ratio in AD (%)
A3SS	4474	2932	14	10.71	14.04
A5SS	2605	1586	11	7.14	12.28
SE	44018	34816	70	66.07	57.89
RI	3377	2046	15	12.50	14.04
MXE	6231	5189	3	3.57	1.75
Total	60705	46569	113	100	100

A3SS, alternative 3′ splice site; A5SS, alternative 5′ splice site; SE, skipped exon; RI, retained intron; MXE, mutually exclusive exons.

TABLE 2

Top three events in each type of AS pattern.

Id	p Value	FDR	ΔPSI	Gene	Type
5780	1.64E-11	7.36E-08	0.125	Rsrp1	A3SS
5943	1.95E-07	0.00029	0.276	Tmem138	A3SS
1546	5.30E-06	0.002963	0.264	Prdm16	A3SS
2080	8.34E-06	0.00345	0.299	Dab2ip	A5SS
3531	1.31E-05	0.004267	−0.145	Zfp652	A5SS
2196	2.22E-05	0.005772	−0.105	Gria3	A5SS
19005	1.23E-09	8.99E-06	0.151	Dnajc6	SE
25241	1.65E-08	4.55E-05	−0.297	1500004A13Rik	SE
31257	2.26E-08	5.85E-05	0.108	Borcs8	SE
3901	0	0	0.229	Rsrp1	RI
3018	1.07E-05	0.006047	0.278	Farsa	RI
916	1.70E-05	0.007192	0.269	Pced1a	RI
1649	4.44E-16	2.77E-12	−0.523	Tmem234	MXE
1647	4.66E-12	1.45E-08	−0.53	Tmem234	MXE
3622	9.61E-06	0.00998	0.115	Eps15l1	MXE

SE Alternative Splicing Events

SE events were the most common AS events, accounting for 61.95% of all significant AS events (Figure 1B). GO analysis of genes with significant SE patterns were enriched in histone acetyltransferase (HAT) complex including Inhibitor of growth family member 3 (Ing3) and Lysine acetyltransferase 5 (Kat5) (Figure 2A,B). Both Ing3 and Kat5 tended to harbor longer exons in AD (Figure 2C,D).

FIGURE 2

SE alternative splicing events. (A) Significant GO terms involving genes with significant SE events. (B) cnetplot showed the relation of involved genes and GO terms. Both Ing3 and Kat5 colored in red showed increased PSI level. (C,D) The detailed sashimi plots for Ing3 and Kat5. Red arrow represented increased PSI.

RI Alternative Splicing Events

RI is usually associated with decreased protein translation. The PSI referred to the inclusion level of the retained intron. A total of 15 (13.27%) significant RI events were identified. GO analysis revealed enrichment of mitochondrion autophagy, involving Optineurin (Optn) and Beclin 1 (Becn1) (Figure 3A,B). Optn showed an increased intron inclusion level (Figure 3C), while Becn1 showed decreased retained intron (Figure 3D). Becn1 plays a central role in autophagy (Wei et al., 2008).

FIGURE 3

RI alternative splicing events. (A) Significant GO terms involving genes with significant RI events. (B) cnetplot showed the relation of involved genes and GO terms. Becn1 colored in green indicated decreased RI PSI levels, while Optn colored in red represented increased PSI level. (C,D) The detailed sashimi plots for Optn and Becn1. The red arrow represented increased PSI, while green arrow indicated decreased PSI.

A3SS Alternative Splicing Events

A3SS occurred due to alternative acceptor sites, expressing shorter or longer exons. In A3SS, the PSI indicates the inclusion level for the longer exon. A total of 14 significant A3SS events were identified (Figure 1B and Supplementary Table S1). GO analysis was performed for A3SS related genes. Interestingly, genes with significant A3SS pattern were enriched in methyltransferase activity, including TRNA Methyltransferase 1 (Trmt1), and PR domain containing 16 (Prdm16) (Figure 4A,B). The sashimi plot showed decreased PSI for Trmt1, indicating a shorter exon of Trmt1 in AD (Figure 4C). Prdm16 showed increased PSI and a longer exon of Prdm16 in AD (Figure 4D). Trmt is a tRNA-modifying enzyme, which acts as a di-methyltransferase using S-adenosyl methionine as a methyl donor to modify the 26th guanine residue of the tRNA. Trmt1 was associated with neurological dysfunction (Jonkhout et al., 2021). Prdm16 displayed histone methyltransferase activity and functioned as a transcriptional regulator (Nishikata et al., 2003).

FIGURE 4

A3SS alternative splicing events. (A) Significant GO terms involving genes with significant A3SS events. (B) cnetplot showed the relation of involved genes and GO terms. Trmt1 colored in green indicated decreased A3SS PSI level, while Prdm16 colored in red represented increased PSI level. (C,D) The detailed sashimi plots for Trmt1 and Prdm16. The red arrow represented increased PSI, while the green arrow indicated decreased PSI.

Differentially Expressed Genes and Transcript

Although aberrant mRNAs arising from AS would be degraded via the nonsense-mediated mRNA decay (NMD) pathway. However, in some cases, these mRNAs may bypass NMD and be translated, especially for RI events (Forrest et al., 2004). We also examined the differentially expressed genes and transcripts in AD. These samples were well correlated with the group (Figure 5A). 404 up- and 370 down-regulated genes were identified (Figure 5B). The expression of genes with significant AS events did not show a significant difference between AD and control on gene level (Figure 5C). On transcript level, 992 up- and 908 down-regulated transcripts were identified. Several genes with significant AS events harbored transcripts with opposite changes, including Trmt1, Ing3, Rsrp1, and Tmem138 (Figure 5D; Table 3). Given Tmem138 for instance, the A3SS event resulted in a shorter transcript (ENSMUST00000235990) which underwent NMD, and a longer transcript (ENSMUST00000025568) with protein coding ability (Figure 5E). In the A3SS event, Tmem138 showed increased PSI for the longer transcript (Figure 5F). The gene level expression showed no significant difference in Tmem138 (Figure 5G). However, the longer ENSMUST00000025568 was increased while the shorter ENSMUST00000235990 was decreased in AD (Figure 5H,I), consistent with the A3SS result (Figure 5F). These results suggest that AS may alter the gene expression on transcript level rather than gene level.

FIGURE 5

Differentially expressed genes and transcripts. (A) Sample correlation analysis using gene expression matrix. (B) Heatmap of dysregulated genes between AD and control groups. (C) Volcano plot of dysregulated genes. The black circle represented genes with significant AS event, and these genes did not show significant change on gene level. (D) Volcano plot of dysregulated transcripts. The labeled genes harbored different transcripts with opposite changes. (E) Genomic information of Tmem138. The red arrow indicated the longer transcript with coding ability, and the green arrow revealed the shorter transcript which may undergo NMD. (F) The detailed sashimi plots for Tmem138, and the red arrow represented increased PSI. (G–I) Relative expression of Tmem138 in gene and transcript levels, respectively. The longer ENSMUST00000025568 was increased while the shorter ENSMUST00000235990 was decreased in AD.

TABLE 3

Dysregulated transcripts involved in significant AS genes.

Id	Gene	Transcript	logFC	p Value	Change
18692	Borcs8	ENSMUST00000123760	1.020448	0.009799	UP
15339	Borcs8	ENSMUST00000110139	1.051684	0.029002	UP
29162	Btbd8	ENSMUST00000152474	−0.79053	0.029368	DOWN
39178	Chfr	ENSMUST00000197968	0.826185	0.028575	UP
15538	Dcun1d2	ENSMUST00000110839	−1.07889	0.00516	DOWN
11548	Ddhd1	ENSMUST00000087320	−0.60375	0.002007	DOWN
14296	Dnajc6	ENSMUST00000106930	0.618619	0.00106	UP
46691	Ermard	ENSMUST00000226599	0.94063	0.016823	UP
17011	Ing3	ENSMUST00000115389	0.899232	0.003086	UP
23256	Ing3	ENSMUST00000136200	−0.63228	0.017412	DOWN
33408	Mdga1	ENSMUST00000168044	1.104151	0.034109	UP
19211	Optn	ENSMUST00000125203	−0.62723	0.033085	DOWN
37077	Pdzd2	ENSMUST00000187398	−1.40437	0.043823	DOWN
11483	Pus1	ENSMUST00000086643	1.62418	0.031663	UP
10747	Rsrp1	ENSMUST00000078084	−0.72837	0.000284	DOWN
26558	Rsrp1	ENSMUST00000145364	0.988984	0.03247	UP
40879	Sema6c	ENSMUST00000204709	−1.30559	0.023373	DOWN
34638	Siglech	ENSMUST00000173835	0.826854	0.047645	UP
6304	Sik2	ENSMUST00000041375	−0.69521	0.03334	DOWN
31551	Tle3	ENSMUST00000160882	0.981684	0.038223	UP
31491	Tle3	ENSMUST00000160724	0.727575	0.046688	UP
48977	Tmem138	ENSMUST00000235990	−1.62863	0.000732	DOWN
2774	Tmem138	ENSMUST00000025568	1.354949	0.005409	UP
21420	Tor1aip2	ENSMUST00000131359	1.211702	0.023502	UP
35047	Trmt1	ENSMUST00000175980	1.207114	0.006327	UP
35484	Trmt1	ENSMUST00000177531	−1.04671	0.016455	DOWN
27869	Zfp652	ENSMUST00000148945	0.819442	0.014392	UP

Discussion

AS is a widespread phenomenon, generating multiple isoforms of a single gene. Here, we investigated the main five AS patterns in the cerebral cortex of 9-month APP/PS1 mice. SE was the most abundant AS event (61.95%), and RI was the second most abundant AS type (13.27%), followed by A3SS (12.39%), then A5SS and MXE comprised 12.39%.

SE is reported to be the most common and widely investigated AS event due to shifting of the open reading frame or loss of functional domains/sites, leading to numerous diseases and considered therapeutic targets, and a variety of SE events have been found in AD human and transgenic mice (Yang et al., 2021). Apolipoprotein E receptor 2 generated an isoform lacking exon 19 though SE and was associated with impairments of spatial learning and long-term memory storage in AD (Forrest et al., 2004; Hinrich et al., 2016). Antisense oligonucleotide targeting intronic splicing silencer, which increased exon 19 inclusion, improved the mentioned cognitive defects in AD mice (Hinrich et al., 2016), providing novel therapeutic strategies for AD. We observed that HAT complex-related genes were spliced via the SE pathway. 16 proteins localized to the h4 HAT complex including Ing3 and Kat5(Rouillard et al., 2016). Both Ing3 and Kat5 showed increased exon inclusion of the skipped exons, leading to the production of longer transcripts. Ing3 activates p53 trans-activated promoters and interacts with TP53 to suppress cell growth and induce apoptosis (Wu et al., 2020). Kat5 is a major histone acetyltransferase that plays important roles in the regulation of DNA repair, autophagy, proteasome-dependent protein turnover, and learning and memory. Disruption of Kat5 mediated histone acetylation was an early common event in neurodegenerative disorders including AD (Beaver et al., 2020; Li and Rasmussen, 2020).

RI is the least understood mode of AS and ranked the second most common AS event in this study. Most retained introns were degraded by the NMD pathway due to the premature termination codons in the newly formed transcript. However, some mRNAs may bypass NMD and be translated (Beaver et al., 2020; Li and Rasmussen, 2020). Li HD et al performed integrative functional genomic analysis of RI in AD and observed remarkable correlations of RI transcripts within innate immune genes. They further identified splicing-related genes which may regulate RI events in AD (Beaver et al., 2020; Li and Rasmussen, 2020). Here, we found that genes spliced by RI were enriched in autophagy. Becn1 is the first described mammalian autophagy gene and plays a central role in autophagy. Mitochondria selective autophagy plays a critical role in various biological processes including the elimination of the damaged mitochondria. Cheng B et al reported that a short Becn1 isoform, which resulted from exon skipping of exon 10–11, was indispensable for mitochondria-selective autophagy, while the full-length Becn1 was essential for nonselective macroautophagic induction (Cheng et al., 2015). Whether novel Becn1 isoforms derived from RI contribute to autophagy regulation deserves further investigation.

On the other hand, A3SS and A5SS are relatively poorly characterized, although they were found to be associated with several diseases due to aberrant splicing. During A3SS and A5SS events, exons are flanked on one constitutive splice site and are flanked by two competing splice sites on the opposite side, leading to an included or excluded alternate region within the transcript. Here we found that A3SS was enriched in methyltransferase activity. Prdm16 is a transcription factor that contains an N-terminal PR domain. Discriminative AS events of Prdm16 transcripts were noted. Exon16 skipped Prdm16 showed a stronger effect than full-length isoform on promoting the transcriptional activity of the PGC-1α promoter (Chi and Lin, 2018). These results support the diverse functions for differentially spliced isoforms.

AS may result in the expression changes of corresponding transcripts or genes. We identified 404 up- and 370 down-regulated genes. On the other hand, we identified 992 up- and 908 down-regulated transcripts. Moreover, genes with significant AS events did not show a significant change of expression on gene level (Figure 5C), however, some of these genes harbored different transcripts with opposite changes on transcript level (Figure 5D). More importantly, we found that Tmem138 showed increased A3SS of ENSMUST00000025568 transcript, which was positively correlated with the expression of this transcript. The decreased transcript of Tmem138 (ENSMUST00000235990) showed reduced transcript expression, while Tmem138 showed no significant change in gene level (Figure 5E–I). These data suggest that AS may cause expression changes on the transcript level. We have also explored the potential splicing factors (SFs) involved in AS (Sveen et al., 2011), however, we did not find any SFs with significant change on gene level (Supplementary Figure S1).

In conclusion, we revealed ontology-specific AS changes in AD. SE genes were enriched in HAT complex, RI genes were enriched in autophagy, and A3SS genes were enriched in methyltransferase activity. Our analysis provides novel pathological mechanisms of AD.

Abbreviations

AS, alternative splicing; 5′SS, 5′ splice site; BP, branch point; 3′SS, 3′ splice site; snRNP, small nuclear ribonucleoprotein; A3SS, alternative 3′SS; A5SS, alternative 5′SS; SE, skipped exon; RI, retained intron; MXE, mutually exclusive exons; AD, Alzheimer’s disease; PSEN1, presenilin 1; SORL1, sortilin-related receptor; TREM2, and triggering receptor expressed on myeloid cells 2; PSI, percent spliced in; GO, gene ontology; DEGs, differently expressed genes; PCA, principal component analysis; HAT, histone acetyltransferase; Ing3, Inhibitor of growth family member 3; Kat5, Lysine acetyltransferase 5 (Kat5); Optn, Optineurin; Becn1, Beclin 1; Trmt1, TRNA Methyltransferase 1; Prdm16, PR domain.