Bulk RNA-seq expression profiles of genes and transcripts were retrieved from the Religious Orders Study/Memory and Aging Project (ROSMAP) (Raj et al., 2018). A total of 454 brain samples, including 145 non-demented control samples and 309 disease samples. The disease samples were further divided into AsymAD (n = 140) and AD samples (n = 169). AsymAD was considered as the early stage of AD (Jack et al., 2018). AD samples were from sporadic AD cohorts (Haure-Mirande et al., 2019). Postmortem neuropathological evaluation of neuritic plaque distribution was performed according to the CERAD score which is a semiquantitative measure of neuritic plaques, and the severity of neurofibrillary tangle (NFT) pathology was assessed with the Braak staging system. Other neuropathological diagnoses were made in accordance with established criteria and guidelines. Sample classification harmonization was performed according to the previous study (Johnson et al., 2020). Particularly, samples with CERAD score 3–4 and Braak 0–3 without dementia at last evaluation were defined as non-demented control; samples with CERAD score 1–3 and Braak 3–6 without dementia at last evaluation were defined as AsymAD; samples with CERAD 1–2 and Braak 3–6 with dementia at last evaluation were defined as AD. Dementia was defined as MMSE < 24. These datasets were downloaded from the AMP-AD Knowledge Portal.¹ The expression profiles were normalized to control different sequencing strategies and sequencing depths in values of transcript per million mapped reads (TPM). Both gene and transcript expression profiles were filtered by removing low-expression genes or transcripts. Specifically, transcripts or genes with expression level >0.1 TPM in no less than 5% of all samples were kept for further analysis.

The differentially expressed transcript (DET) analysis was performed between the AsymAD and non-demented control, AD and non-demented control samples, respectively. The Wilcoxon’s rank-sum test was employed in the difference evaluation. A transcript with fold change >1.5 or < 0.67 and FDR value <0.05 was considered as DET in AsymAD or AD samples. Protein-coding genes with DETs were adopted to perform functional enrichment analysis by using the clusterProfiler R package (version 4.1.4) (Wu et al., 2021). Gene lists of biological processes were derived from the Molecular Signature Database (MSigDB) (Liberzon et al., 2015).

For each transcript, a Pearson correlation coefficient was computed with the corresponding host gene by using normalized expression values across all samples. Host genes are those genes that the corresponding transcripts are transcribed from, and transcript-host gene pairs were extracted from the GENCODE annotation (Frankish et al., 2019). The expression correlation between a transcript and host gene was considered positive if their Pearson correlation coefficient was higher than 0.3, negative if their correlation coefficient was lower than -0.3. The threshold of p-values was set at 0.05. Otherwise, the expression of a transcript was considered as irrelevant with its host gene.

The alternative splicing events (ASEs) were identified by using the SUPPA2 (version 2.3) software (Trincado et al., 2018). Briefly, ASEs were inferred from the comparison of transcript structures and expression levels in the same genes. The value of percent spliced in (PSI) was calculated for each ASE in each sample (Li et al., 2019). ASEs were further processed to generate high-confidence events by retaining events with a PSI value greater than 0.1 in no less than 5% of samples. To detect differentially spliced ASEs, we used the Wilcoxon singed-rank test to compare the PSI values between non-demented control and AsymAD or AD samples. ASEs with a median PSI > 0.1 and a corrected p-value < 0.05 were considered differentially spliced.

Isoform switching events of individual genes were identified in each sample. Briefly, summed expression of each gene was calculated from all transcripts (isoforms) that belong to corresponding host genes. The isoform fraction (IF) of each isoform was calculated as follows:

where IF_i represents the IF of isoform i, T_i is the expression level of isoform i, E_j denotes the expression level of isoforms in gene j, n indicates the number of isoforms in gene j.

The significance of the change in isoform usages between AD and non-demented control, or AsymAD and non-demented control group was evaluated by Student’s t-test. The difference in IF of each isoform was calculated as follows:

where dIF_i is the difference in IF of isoform i, IF_i denotes the IF value of isoform i in group k, n is the sample number of group k, IF_j denotes the IF value of isoform i in group j, m is the sample number of group j.

A gene that has at least one isoform with | dIF| > 0.1 and FDR < 0.05 was considered to have isoform switching events in AsymAD or AD samples.

A list of RBPs were retrieved from the Encyclopedia of RNA Interactomes (ENCORI) database (Li et al., 2014). ENCORI database also provided RBP binding regions that were generated from the integration of large-scale high-throughput cross-linking immunoprecipitation (CLIP) experiments. Spearman correlations were calculated between each RBP and differentially spliced ASE. In the correlation calculation, we only considered RBP-ASE pairs that RBPs had binding evidence of CLIP signals within regions +/-300 bp around the corresponding splicing sites. RBPs that showed no less than 10% expression changes between non-demented control samples and AsymAD or AD samples were considered in the correlation calculation. RBP-ASE pairs that had a Spearman correlation >0.5 or <-0.5 and p-value < 0.05 were considered significantly correlated.

Statistical analysis and data visualization in this study were performed by using the R software (R Foundation for Statistical Computing, Vienna, Austria).² Unless specifically stated, all tests were two-sided, and p or FDR values <0.05 were considered statistically significant.

To characterize a high-resolution dysregulation in AD transcriptome, we performed transcript-level investigation in 454 AD-related samples. In terms of median value, 77,258, 75,253, and 75,355 transcripts were detected in the non-demented control, AsymAD, and AD samples, respectively (Figure 1A). The vast majority (n = 125,426, 84.5%) of these expressed transcripts were derived from protein-coding genes, while approximately 10 percent (n = 13,924, 9.4%) were from long non-coding RNA (lncRNA) genes and 5.3% of transcripts were from pseudogenes (Figure 1B). Most genes expressed multiple transcripts, especially protein coding genes, wherein over 80% of protein coding genes expressed more than five transcripts (Figure 1C). More than 60% of lncRNA genes expressed no less than two transcripts, while about one third pseudogenes expressed at least two transcripts. We next identified DETs in AsymAD and AD samples, respectively. A total of 1,039 upregulated transcripts and 1,545 downregulated transcripts were identified in AD samples compared to non-demented control samples (Figure 1D and Supplementary Table 1). In AsymAD samples, 135 upregulated and 144 downregulated transcripts were identified (Supplementary Table 2). The smaller number of DETs in AsymAD indicated that AsymAD might have lower level of molecular pathological changes compared to AD. Some dysregulated transcripts were shared between AsymAD and AD, while some were specific in AsymAD or AD samples. Specifically, 1,470 transcripts were downregulated in AD but showed no change in AsymAD samples, and 975 transcripts showed no change in AsymAD but upregulated in AD samples (Figure 1E). There were 98 transcripts that were downregulated and 89 transcripts that were upregulated in AsymAD but exhibited no change in AD samples. Among these dysregulated transcripts, 43 were upregulated, whereas 41 were downregulated in both AsymAD and AD samples. Further enrichment analysis revealed that upregulated DETs in AsymAD samples were enriched in metabolism-related processes, such as “Regulation of steroid biosynthetic process,” “Regulation of lipid biosynthetic process,” and “Regulation of steroid metabolic process,” while downregulated DETs showed enrichment of localization and transport processes, such as “Protein localization to nucleus” and “Nuclear export” (Figure 1F). Upregulated DETs in AD samples showed enrichment in synaptic functions, such as “Regulation of trans-synaptic signaling,” “Modulation of chemical synaptic transmission,” and “Synapse organization,” whereas downregulated DETs were enriched in cellular response processes, such as “Cellular response to abiotic stimulus,” “Cellular response to environmental stimulus,” and “Cellular response to zinc ion” (Figure 1G). In addition, we compared the AsymAD and AD transcriptomes and identified 95 upregulated and 37 downregulated transcripts in AD samples (Supplementary Figure 1A, Supplementary Table 3). The upregulated transcripts were enriched in “Positive regulation of cell-matrix adhesion,” “Cell cycle arrest,” and “Iron-ion transport,” while downregulated transcripts were enriched in “Negative regulation of cell-matrix adhesion,” “Negative regulation of cysteine-type endopeptidase activity involved in apoptotic process,” and “Negative regulation of blood vessel endothelial cell migration” (Supplementary Figure 1B). In summary, we characterized AsymAD and AD transcriptome and its dysregulation at transcript resolution.

To investigate whether most transcripts showed consistent expression patterns with their host genes, we next calculated expression correlations between transcripts and the corresponding host genes (see Methods). Approximately half of all transcripts showed irrelevant expression correlation with their host genes (Figure 2A). As expected, the number of irrelevant transcripts increased with larger number of transcripts in single genes (Figure 2B). These observations highlighted the importance of transcript-level investigation of AD transcriptome. In addition, the expression fractions of positively correlated transcripts decreased with the increasing numbers of transcripts in single genes. More specifically, the median expression fraction of positively correlated transcripts was 0.63 in genes that had two transcripts (Figure 2C), while the median values of expression fractions were all below 0.5 in genes that expressed more than two transcripts. For example, the gene FUCA2 expressed four transcripts, i.e., ENST00000002165.5, ENST00000438118.2, ENST00000451668.1, and ENST00000367585.1 (Figure 2D). Transcript ENST00000002165.5 showed high expression correlation with host gene FUCA2, while neither of the other three transcripts exhibited any expression correlation. Conclusively, our results highlighted that transcript-level analysis revealed a higher resolution of AD transcriptome than those in gene level.

Alternative splicing is supposed to make a major contribution to transcriptional variations in human diseases. We next identified ASEs in all samples, including seven different ASE types, i.e., alternative 3’ splice site (A3), alternative 5’ splice site (A5), alternative first exon (AF), alternative last exon (AL), mutually exclusive exons (MX), retained intron (RI), and skipping exon (SE). As expected, we identified the most ASEs in protein coding genes (Figure 3A). AF was the most frequent (n = 44,960, 42.2%) ASE type in protein coding genes, followed by SE (n = 24,088, 22.6%) (Figure 3B). In long non-coding RNA (lncRNA) genes, AF (n = 1,598, 29.5%) and AL (n = 1,371, 25.3%) were the most frequent ASE types. In pseudogenes, SE was the ASE type with the largest number. We next examined whether sequences affected by ASEs had intact codons. In all ASE types, over half sequences did not have intact codons, indicating that most ASEs might cause frame shift (Figure 3C). We also examined these in separate sample groups and found that the number distribution of different ASE types (Supplementary Figure 2), percentages of different ASE types (Supplementary Figure 3), and intact codons (Supplementary Figure 4) showed no notable differences in non-demented control, AsymAD, and AD sample groups. Comparing with non-demented control samples, AD samples showed more differentially spliced ASEs than AsymAD (Figure 3D). In total, 4,056 and 1,200 differentially spliced ASEs were identified in AD and AsymAD samples, respectively (Supplementary Tables 4, 5). AF and SE events had the largest number of differentially spliced events in both AsymAD and AD samples. Differentially spliced genes (DSGs) in AsymAD were enriched in “Adipogenesis,” “HEME metabolism,” and “MYC targets,” while DSGs in AD showed enrichment in “Myogenesis,” “Mitotic spindle,” and “Fatty acid metabolism” (Figure 3E). These results revealed the differences of alternative splicing variations between AsymAD and AD.

To further investigate the transcript variations in AsymAD and AD, we examined transcripts that showed greater than 20% changes compared to the non-demented control samples. A large number of transcripts showed consistent expression changes (>20% higher or lower) in both AsymAD and AD samples (Figure 4A). There were 5,229 upregulated and 4,843 downregulated transcripts that showed expression changes in only AD samples, whereas 4,843 transcripts were downregulated and 2,584 transcripts were upregulated in only AsymAD samples. In addition, 1,340 transcripts showed downregulated expression in AsymAD but upregulated expression in AD, while 581 transcripts were upregulated in AsymAD but downregulated in AD. Next, we identified isoform switching events in AsymAD and AD samples, respectively (see Methods). A total of 163 and 119 transcripts showed increased usage, while 124 and 103 transcripts exhibited decreased usage in AsymAD and AD, respectively (Figures 4B, C and Supplementary Table 6). We performed enrichment analysis of dysregulated genes that showed isoform switching events. Those genes that had isoform switching events in AsymAD samples were enriched in biological processes of visual system, such as “Eye development,” “Visual system development,” and “Retinal rod cell development” (Supplementary Figure 5A). In AD samples, genes that had isoform switching events were enriched in regulation of synaptic functions, such as “Modulation of chemical synaptic transmission,” “Regulation of trans-synaptic signaling,” and “Regulation of postsynaptic membrane potential” (Supplementary Figure 5B). For example, in the OPN4 gene, transcript ENST00000372071.2 had one more exon than transcript ENST00000241891.5 (Figure 4D). The OPN4 gene, a photosensitive protein, not only participates in the non-visual imaging system of animals (Tian and Copenhagen, 2003), but also in the regulation of circadian rhythm (Van Gelder, 2001). Circadian rhythm dysfunction and sleep disturbances are associated with aging and neurodegenerative diseases, including mild cognitive impairment (MCI) and AD. The OPN4 gene has been reported to be associated with AD (Bacalini et al., 2022), which showed no expression changes between AsymAD and non-demented control samples. Transcript ENST00000241891.5 showed significant upregulation, while transcript ENST00000372071.2 was significantly downregulated in AsymAD samples. Additionally, among transcripts expressed from gene OPN4, the usage of transcript ENST00000241891.5 was significantly increased, while the usage of transcript ENST00000372071.2 decreased in AsymAD. In the APOA2 gene, transcript ENST00000367990.3 expressed four exons, whereas transcript ENST00000463812.1 had three exons (Figure 4E). Transcript ENST00000367990.3 was upregulated, whereas transcript ENST00000463812.1 showed downregulation in AD samples. Compared to non-demented control samples, gene APOA2 expressed higher proportion of transcript ENST00000367990.3 in AD. The regulation and metabolism of lipid have been demonstrated to play important roles in the disease progression of AD (Yin, 2023). The APOA2 gene belongs to the apoprotein gene family, wherein APOE has been demonstrated to play roles in AD. The switch between transcript ENST00000367990.3 and ENST00000463812.1 indicates the potential regulatory role of the APOA2 gene in the progression of AD. In addition, we examined all isoform switching events in previously reported AD-related risk genes (Keren-Shaul et al., 2017; Grubman et al., 2019; Leng et al., 2021; Mertens et al., 2021; Bellenguez et al., 2022), and found 34 risk genes have at least one isoform switching event. For example, the APP gene showed switch between transcript ENST00000348990.5 and ENST00000354192.3 (Supplementary Figure 6), which was reported to have AD-related risk loci in GWAS analysis (Bellenguez et al., 2022). In summary, these results revealed the transcriptional reprograming in AsymAD and AD.

To further examine the regulatory role of RBPs, we built RBP-ASE regulatory networks in AsymAD and AD, respectively (see Methods). In the RBP-ASE network of AsymAD, four upregulated (QKI, RBM47, RBM6, and PAPD5) and three downregulated (KHDRBS2, ALYREF, and FTO) RBPs modulated 5,190 ASE events (Figure 5A and Supplementary Table 7). RBM6 regulated the largest number of ASEs, while FTO regulated the most pathways (Figure 5B). For example, gene ALG2 was predicted to be regulated by five RBPs (FTO, KHDRBS2, ALYREF, RBM6, and PAPD5) in AsymAD, which expressed three transcripts, i.e., ENST00000476832.1, ENST00000319033.6, and ENST00000238477.5 (Figure 5C). Transcript ENST00000319033.6 showed upregulation and increased usage, while transcript ENST00000476832.1 exhibited downregulation and decreased usage in AsymAD. In the RBP-ASE regulatory network of AD, 15 upregulated (LARP7, DKC1, QKI, RBM47, TROVE2, SLTM, MOV10, ACIN1, PTBP1, HNRNPA2B1, FBL, RBM6, ZFP36, MSI2, and VIM) and 10 downregulated (GNL3, KHDRBS2, ZNF184, LIN28B, LSM11, ALYREF, RBFOX2, DDX54, KHDRBS1, and YWHAG) RBPs modulated 9,918 ASE events (Figure 5D and Supplementary Table 8). KHDRBS2 regulated the largest number of ASEs in AD, whereas HNRNPA2B1 regulated the most pathways that were enriched by differentially spliced ASEs (Figure 5E). For instance, multiple RBPs were predicted to regulate gene HTR1E in AD (Figure 5F). In gene HTR1E, transcript ENST00000369584.1 expressed one more exon than transcript ENST00000305344.4. The expression level and usage of transcript ENST00000369584.1 significantly increased, while those of transcript ENST00000305344.4 decreased in AD.