Cross-platform transcriptional profiling identifies common and distinct molecular pathologies in Lewy body diseases.

Parkinson's disease (PD), Parkinson's disease with dementia (PDD) and dementia with Lewy bodies (DLB) are three clinically, genetically and neuropathologically overlapping neurodegenerative diseases collectively known as the Lewy body diseases (LBDs). A variety of molecular mechanisms have been implicated in PD pathogenesis, but the mechanisms underlying PDD and DLB remain largely unknown, a knowledge gap that presents an impediment to the discovery of disease-modifying therapies. Transcriptomic profiling can contribute to addressing this gap, but remains limited in the LBDs. Here, we applied paired bulk-tissue and single-nucleus RNA-sequencing to anterior cingulate cortex samples derived from 28 individuals, including healthy controls, PD, PDD and DLB cases (n = 7 per group), to transcriptomically profile the LBDs. Using this approach, we (i) found transcriptional alterations in multiple cell types across the LBDs; (ii) discovered evidence for widespread dysregulation of RNA splicing, particularly in PDD and DLB; (iii) identified potential splicing factors, with links to other dementia-related neurodegenerative diseases, coordinating this dysregulation; and (iv) identified transcriptomic commonalities and distinctions between the LBDs that inform understanding of the relationships between these three clinical disorders. Together, these findings have important implications for the design of RNA-targeted therapies for these diseases and highlight a potential molecular "window" of therapeutic opportunity between the initial onset of PD and subsequent development of Lewy body dementia.

Introduction

The Lewy body diseases (LBDs) comprise three neurodegenerative diseases, which are characterised by accumulation of Lewy bodies (α-synuclein-containing aggregates) in neurons and neuronal processes [55, 84]. These disorders, which include Parkinson’s disease (PD), Parkinson’s disease with dementia (PDD) and dementia with Lewy bodies (DLB), have a prevalence in the general population aged ≥ 65 years of 2–3% [87], 0.3–0.5% [3] and 1–2% [55], respectively. Together, PDD and DLB are collectively known as the Lewy body dementias and they are second only to Alzheimer’s disease (AD) in prevalence among people with dementia [4]. All three LBDs are associated with disability and reduced quality of life; DLB is associated with earlier mortality and a higher cost of care compared with AD [18, 83, 107]. With no disease-modifying therapies available for any of the LBDs, these diseases present a major unmet clinical need [95].

While a variety of mechanisms, including mitochondrial and lysosomal dysfunction, oxidative stress, α-synuclein misfolding and neuroinflammation, have been implicated in PD pathogenesis [17, 87], less is known about the mechanisms underlying PDD and DLB. Elucidating these mechanisms could provide a biological basis for the clinical distinction between PDD and DLB, which remains controversial in the field [14, 55, 56, 89, 112]. Clinically, PDD and DLB are arbitrarily separated by the diagnostic "1-year rule": if dementia is diagnosed before or within 1 year of the onset of parkinsonism, it is considered to represent DLB, whereas PDD is defined by dementia first presenting more than 1 year after the onset of parkinsonism [38, 74]. Thus, PDD and DLB are clinically distinguished based only on the relative timing of motor and cognitive impairments, despite sharing many symptoms (e.g. dementia, depression, parkinsonism, REM sleep behaviour disorder and visual hallucinations). Arguably, two of the core clinical features of DLB, fluctuating cognition and visual hallucinations, are more prevalent in DLB compared with PD/PDD [39, 81], suggesting two separate disorders. However, the overlap of these core clinical features could also be evidence that the disorders are on a spectrum of disease, where DLB represents a more severe form of PDD.

Neuropathologically, all three LBDs are classed as synucleinopathies, but at the end stage of disease they often present with concomitant pathologies, such as tau neurofibrillary tangles and amyloid-β [44, 92, 99]. PD and PDD are thought to be purer synucleinopathies, whereas over 90% of DLB cases have some, often substantial, AD pathology [38, 52, 71, 74, 92, 99]. While some neuropathological differences have been reported between the Lewy body dementias and PD (e.g. tau and amyloid-β pathology at a more advanced stage in the Lewy body dementias [99]), these differences do not permit confident distinction between the LBDs when no clinical diagnosis is present. Genetically, the differences between PDD and DLB are not well-characterised, although APOE, GBA and SNCA mutations have been implicated in both [2, 112]. More is known about the genetic risk factors contributing to PD and DLB, which share some risk loci (GBA, TMEM175 and SNCA) and pathways (lysosomal and endocytic pathways) [21, 29, 50, 77, 93]. However, there is also evidence that association signals at SNCA may be distinct in PD and DLB (i.e. located at the 3’ and 5’ end of the SNCA gene, respectively) [21, 29, 48, 50], and while risk pathways are shared, PD genetic risk factors only explain a small portion of DLB phenotypic variance [29, 49].

Identifying therapeutic targets that could modify the development of PDD or DLB requires an understanding of the cellular and molecular features of these diseases. Transcriptomic profiling, through RNA-sequencing of patient-derived tissue, would aid in the identification of such targets, but remains limited in all three LBDs. Of all transcriptomic studies of PD and Lewy body dementia highlighted in two recent systematic reviews (33 and 31 gene expression studies in brain, respectively [17, 30]), only 5 used RNA-sequencing. Furthermore, among transcriptomic studies of the three LBDs, few have addressed possible alternative splicing or the confounding of bulk-tissue transcriptomic profiling by differences in cellular composition.

Here, we pair bulk-tissue and single-nucleus RNA-sequencing to gain a comprehensive view of cell-type-specific transcriptional changes in the LBDs. This combined approach is used, because, while single-nucleus RNA-sequencing can address confounding by cellular composition, providing previously unattainable insight into cell-type-specific transcriptomic pathology [60, 61], compared with bulk-tissue RNA-sequencing it has little ability to resolve transcriptomic diversity via splicing. This limitation arises due to the trade-off that exists between choosing a single-nucleus RNA-sequencing protocol that has high throughput but only sequences 3′ ends of transcripts versus a protocol whose library construction permits sequencing full-length transcripts but has reduced throughput [27]. Using this combined sequencing approach, we found transcriptional changes in multiple cortical cell types across the LBDs, with more differentially expressed genes and pathways identified in PDD and DLB than in PD. We also observed widespread alternative splicing, particularly in PDD and DLB, with evidence suggesting that specific splicing factors play a role in orchestrating the disease-related splicing changes. Collectively, these results identify common and distinct molecular pathology in the LBDs across several cell types and provide insight into the extent to which the LBDs represent discrete diseases with unique pathogenic processes.

Results

Paired single-nucleus and bulk-tissue RNA sequencing of anterior cingulate cortex in individuals with Lewy body disease

We applied single-nucleus and bulk-tissue RNA-sequencing to adjacent anterior cingulate cortex tissue sections from 28 individuals, including non-neurological control individuals and individuals with Lewy body disease (Fig. 1). The latter were split into three disease groups, consisting of PD, PDD and DLB, based on clinical assessments of retrospectively reviewed case records (n = 7 per group). We sampled from the anterior cingulate cortex, as it is one of the first cortical areas to be affected by α-synuclein pathology [6, 105] and a region where Lewy body densities correlate with cognitive impairment in PD [59]. Although selected individuals were matched, where possible, for demographic and pathologic factors, there were significant differences in the proportions of sexes between the groups in keeping with previous literature describing a male bias in DLB [78] (proportion female: control = 1/7, PD = 5/7, PDD = 2/7, DLB = 0/7; p value = 0.020; Chi-squared test; Supplementary Fig. 1, Supplementary Table 1). Disease duration also differed significantly between groups, with DLB cases having a shorter duration of disease before death, reflecting the fact that PDD cases have PD motor symptoms for several years before development of dementia (median disease duration in years: PD = 12, PDD = 11, DLB = 6; p value = 0.0099; Kruskal–Wallis rank sum test; Supplementary Fig. 1, Supplementary Table 1). Using this sample set, we report a total of 205,948 droplet-based single-nucleus and 24 bulk-tissue transcriptomic profiles, with an average of 1,398 genes per nucleus and 27,802 genes per bulk-tissue sample detected, respectively (Supplementary Fig. 2, Supplementary Fig. 3, Supplementary Table 1).

Fig. 1

Overview of approach. In this study, anterior cingulate cortex was sampled from a cohort of 28 individuals divided equally between four groups: non-neurological controls; Parkinson’s disease without cognitive impairment (PD); Parkinson’s disease with dementia (PDD); and dementia with Lewy bodies (DLB) (Supplementary Fig. 1, Supplementary Table 1). For each individual, a frozen tissue block derived from the anterior cingulate was sectioned (sectioned area indicated with green shaded box), with adjacent sections used for single-nucleus or bulk-tissue RNA-sequencing (Supplementary Fig. 2, Supplementary Fig. 3, Supplementary Table 1). Following data pre-processing, single-nucleus RNA-sequencing data was used to generate cell-type-specific differential gene expression profiles and to deconvolute bulk-tissue RNA-sequencing data. Bulk-tissue RNA-sequencing was used in differential gene expression and splicing analyses, with cell-type proportions included as model covariates in both analyses. Results from single-nucleus RNA-sequencing and bulk-tissue RNA-sequencing were used in downstream gene set enrichment analyses to identify disease-relevant pathways. Furthermore, common risk variants for Alzheimer’s disease (AD), PD risk and PD age of onset (PD AOO) were mapped to cell-type-specific expression profiles and cell-type-specific differential expression. The image of the human brain displays a coronal section cut at the level of nucleus accumbens. RNA-seq RNA-sequencing, UMI unique molecular identifier

Increased proportions of microglia and vascular cells across Lewy body diseases

Quality control, clustering and classification of major cell types in the anterior cingulate cortex was first performed on nuclear RNA from each of the 28 individuals, after which we used the Conos framework to generate a joint graph of nuclei across all individuals [13]. Clusters were assigned to 7 broad cell types by significant overlap (Fisher’s exact test, p value < 2.2 × 10–16) with a merged list of marker genes derived from two human single-cell datasets (Supplementary Fig. 2) [61, 109]. In total, we identified 75,826 excitatory neurons, 26,467 inhibitory neurons, 46,662 oligodendrocytes, 25,726 astrocytes, 13,788 microglia, 12,497 oligodendrocyte precursors (OPCs), and 4532 vascular cells (which represented a merge of endothelial cells and pericytes), with each cell type consistently identified across all individuals in each disease group (Fig. 2a, Supplementary Fig. 4a, b).

Fig. 2

Cellular diversity of the anterior cingulate cortex across disease states. a Joint graph of all nuclei derived from all individuals visualised using UMAP embedding. Nuclei are coloured by cell type. b Cell-type proportions derived from Scaden deconvolution (available in Supplementary Table 2). Cell-type proportions (upper panel) are grouped by cell type and disease status and displayed relative to the median of controls (within a cell type). Significant differences in cell-type proportions between disease groups (lower panel) were determined using the Wilcoxon rank sum test, with FDR correction for multiple testing. Non-significant results (FDR > 0.1) were coloured white; **FDR < 0.05; *FDR ≤ 0.1. OPC oligodendrocyte precursor cell, UMAP uniform manifold approximation and projection

Next, we sought to identify significant changes in the proportions of these major cell types across all disease groups. Although single-nucleus RNA-sequencing shows less sampling bias than single-cell sequencing [10], its suitability for estimation of cell-type proportions remains in question [33]. Thus, we used Scaden [76], a deep-learning-based deconvolution algorithm that can train on artificial bulk-tissue RNA-sequencing samples simulated from tissue-matched single-nucleus RNA-sequencing data, to estimate cell-type proportions across disease groups. Importantly, Scaden permitted pairing of our single-nucleus and bulk-tissue transcriptomic profiles and modelling of inter-subject variability. We observed a low overall correlation between single-nucleus-estimated and Scaden-predicted cell-type proportions (Spearman’s ρ = 0.25, p value = 0.0009), although per-cell-type correlations were higher for some cell types (highest in microglia, Spearman’s ρ = 0.79, p value = 8.2 × 10–6; Supplementary Fig. 4c).

Using Scaden predictions, we identified a significantly increased proportion of microglia in all disease groups compared with the control group, and a significantly increased proportion of OPCs and vascular cells in DLB cases compared with controls (Fig. 2b, FDR-corrected p < 0.05, Wilcoxon rank sum test). In addition, we observed a nominally significant increase in vascular proportions in PDD and PD cases compared with controls (FDR-corrected p < 0.1, Fig. 2b). By applying Scaden to a second, larger independent PD case–control bulk-tissue RNA-sequencing dataset [37], we were able to replicate the observed increase in microglial and vascular proportions in PD cases compared with controls (FDR-corrected p < 0.05, Supplementary Fig. 5).

Differential gene expression analysis highlights transcriptional alterations in multiple cell types and differentiates Lewy body dementias from PD

Differential gene expression analyses were separately performed with bulk-tissue and single-nucleus RNA-sequencing data to characterise molecular changes across the disease groups (“Materials and methods”). Following correction for changes in Scaden-predicted cell-type proportions in bulk-tissue gene expression, only 60 genes (53 unique genes) were found differentially expressed (DE) across the six pairwise comparisons (FDR < 0.05, Supplementary Table 3). Despite the low number of bulk-tissue DE genes identified, we noted that gene expression adjusted for cell type and experimental covariates resulted in much clearer clustering of samples by disease group (as determined through visual inspection) compared with uncorrected gene expression and gene expression adjusted for experimental covariates alone (Supplementary Fig. 6a–c). Notably, separation of disease groups was primarily observed on the same axis of variation (i.e. the first principal component, PC1), suggesting that (i) the genes contributing most to variation between groups are similar across disease groups, and thus PD, PDD and DLB may represent a neuropathological continuum and (ii) that there are gene expression changes between disease groups that are independent of differences in cell-type proportions (Supplementary Fig. 6a–c). Using pathway enrichment, we found that the top 100 genes contributing to PC1 were associated with immune-related GO terms (e.g. peptide antigen binding and MHC protein complex), as well as terms relating to endocytic vesicles and unfolded protein binding (Supplementary Fig. 6d, Supplementary Table 4).

Consistent with the view that gene expression changes exist between disease groups independent of differences in cell-type proportions, using single-nucleus RNA-sequencing data, 9,242 unique genes were found DE across cell-type-specific pairwise comparisons (all six pairwise comparisons, |log2(fold change)|> log2(1.5), FDR < 0.05, Supplementary Table 5). Focusing only on comparisons with the control group, these analyses highlighted three main themes.

First, differential gene expression was widespread and involved glia and neurons. While we found that DE genes were detected across all three case–control comparisons and across all major cell types, the largest numbers of DE genes were observed in excitatory neurons, followed by oligodendrocytes (Fig. 3a). In fact, across case–control comparisons, the number of DE genes identified in oligodendrocytes exceeded that in inhibitory neurons by a factor of up to 11.4-fold (depending on the case–control comparison; Fig. 3a). Comparison of the Lewy body diseases to each other yielded similar results; that is, transcriptional alterations across all major cell types, but with the largest number of DE genes observed in excitatory neurons, followed by oligodendrocytes (Supplementary Fig. 7).

Fig. 3

Cell-type-specific gene expression changes and pathway enrichments across disease states. a Number of differentially expressed (DE) genes across each cell type in pairwise comparisons of disease groups to the control group (|log2(fold change)|> log2(1.5), FDR < 0.05). The intensity of the grey colour is proportional to the number of DE genes. b Binary plot indicating with bars whether a gene (column) is down-regulated (upper panel) or up-regulated (lower panel) in a given cell type (rows). Number of DE genes in each comparison indicated on the x-axis. c Reduced gene ontology (GO) terms associated with cell-type-specific down- and up-regulated DE genes identified across pairwise comparisons of disease groups with the control group. Due to the magnitude of pathway enrichments, original GO term enrichments (referred to as “child terms”) were grouped using semantic similarity. The number of enriched child GO terms assigned to each reduced parent term across all cell types and comparisons in the panel is indicated in parentheses on the y-axis. Reduced GO terms were ordered on the y-axis by the number of cell types and comparisons in which the term was found enriched. The fill of each tile indicates the − log10(FDR) of the most significant child term associated with the parent term within that comparison/cell type. Non-significant results (FDR > 0.05) were coloured white. Results for pairwise comparisons between disease groups are displayed in Supplementary Fig. 7. All cell-type-specific DE genes and pathway enrichments are available in Supplementary Table 5 and Supplementary Table 6, respectively. DEG differentially expressed gene, OPC oligodendrocyte precursor cell

Second, DE genes were commonly specific to a cell type. Indeed, of the 1131, 2535 and 4816 down-regulated DE genes identified across comparisons of PD, PDD and DLB with control, 79%, 66% and 67%, respectively, were DE in only one cell type (Fig. 3b). Among up-regulated DE genes, these percentages ranged from 74 to 76% across the three case–control comparisons.

Third, the Lewy body dementias, as distinct from PD, were characterised by the predominant down-regulation of gene expression relative to control in most cell types; the only exception were inhibitory neurons in PDD, where the number of up-regulated DE genes exceeded the number of down-regulated DE genes (Fig. 3a, b). Furthermore, the transcriptomic profile of the two Lewy body dementias was very similar, with 303 down-regulated and 87 up-regulated DE genes identified in a comparison of DLB with PDD (Supplementary Fig. 7). In contrast, comparisons of the two Lewy body dementias with PD identified > 2000 down-regulated and > 1000 up-regulated DE genes, suggesting that while there are transcriptional commonalities between PDD and DLB, PD is transcriptionally distinct from the Lewy body dementias in the anterior cingulate cortex.

Pathway enrichment was used to explore the biological implications of cell-type-specific differential gene expression. Focusing on case–control comparisons, we found that down- and up-regulated DE gene sets were enriched for 306 and 272 GO terms, respectively (each pathway was only counted once, even if it appeared across > 1 case–control comparison). Using measures of semantic similarity to cluster GO terms, and thus reduce pathway redundancy, we identified 29 down-regulated and 27 up-regulated GO terms (Fig. 3c, Supplementary Table 6). Despite the high proportion of cell-type-specific DE genes, we identified GO terms that were perturbed across multiple cell types in a given case–control comparison. For example, in comparisons of PD with control, terms related to glutamatergic synapses, the mitochondrial inner membrane, and post-translational protein modification were enriched across ≥ 5 cell types. These commonalities in GO term enrichment were a feature of both down- and up-regulated DE gene sets but were more apparent among (i) down-regulated DE gene sets and (ii) comparisons of PDD and DLB with control, with pathway perturbations affecting a median of 3–5 cell types, as compared with 1–3 in comparisons of PD with control (Supplementary Fig. 8a). Furthermore, we noted that consistent with the high number of DE genes detected for excitatory neurons, a high number of enriched pathways were observed in this cell type across all case–control comparisons, particularly in PDD and DLB (Supplementary Fig. 8b). This observation was even more pronounced in comparisons of the Lewy body dementias with PD, where the number of enriched pathways identified in excitatory neurons was almost twofold higher than the second most-affected cell type. Overall, this analysis served to highlight disproportionately large transcriptional differences in PDD and DLB, as compared with PD, particularly in excitatory neurons and, to a lesser extent, oligodendrocytes.

Genes and pathways genetically associated with PD implicate physiological variability of SNCA expression in selective vulnerability of neurons

Many of the GO terms enriched among down- and up-regulated genes, such as receptor-mediated endocytosis, have been previously implicated in PD. With this in mind, we narrowed our focus to the cell-type-specific expression of genes and pathways genetically associated with PD pathogenesis [12, 16].

PD-associated genes were derived from a recent review of mutations that have been reported to cause PD, including well-known examples, such as SNCA [16]. Of the 21 genes considered, 13 were DE in at least one major cell type and one case–control comparison (Fig. 4a). For example, excitatory neurons, inhibitory neurons, astrocytes and oligodendrocytes all showed significant up-regulation of SNCA in PD cases when compared with controls (fold change: 0.64–1.30; FDR: 2.6 × 10–7–7.2 × 10–157, Fig. 4a).

Fig. 4

Cell-type-specific alterations of PD-associated genes and pathways. a Differential expression of PD-associated genes (associated by mutations reported to cause PD) across cell types and pairwise comparisons of disease groups with the control group. Fill of the tile indicates the log2(fold change), with non-significant results (FDR > 0.05) coloured grey. b UMAP plot of excitatory and inhibitory neurons (upper panel, 102,293 nuclei), with SNCA expression levels (lower panel). c Ridgeline plot of distribution of SNCA expression levels in excitatory neurons across disease groups. Distributions have been split into 3 cumulative quantiles, highlighting, where 0–50%, 50–90% and 90–100% of the nuclei in each disease group lie. d Number of enriched pathways (FDR < 0.05) identified using cell-type-specific down- and up-regulated DE genes from each pairwise comparison together with 46 PD-associated pathways (associated in a large-scale polygenic risk score-based assessment of 2199 gene sets). DEG differentially expressed gene, GO gene ontology, OPC oligodendrocyte precursor cell, UMAP uniform manifold approximation and projection. PD-associated genes and pathways were derived from references [12, 16], respectively

There is robust genetic evidence linking increased SNCA dosage to PD pathogenesis, including (i) duplication and triplication events in the SNCA gene that underlie autosomal dominant forms of PD [26, 97] and (ii) the association of PD risk loci with increased SNCA expression [66, 100]. In view of this evidence, we further explored SNCA expression, finding that, while SNCA expression was up-regulated in PD in all four cell types with a similar fold change (Fig. 4a), SNCA expression in control individuals was highly variable across cell types (Supplementary Fig. 9). This variability in control SNCA expression extended to (i) the proportion of nuclei expressing SNCA, with 61% of excitatory neurons expressing SNCA, as compared with < 22% across all other cell types and (ii) the range of observed SNCA expression, which was wider in excitatory neurons compared with all other cell types (Supplementary Fig. 9). These differences in cell-type-specific SNCA expression were particularly apparent between inhibitory and excitatory neurons, irrespective of disease group, with a higher proportion of excitatory neurons expressing SNCA (Fig. 4b, Supplementary Fig. 9). Furthermore, these differences were visible in a cell type across disease groups. Indeed, SNCA expression in excitatory neurons from the Lewy body dementias, as compared with the control group, was marked by (i) a decrease in the proportion of SNCA-expressing nuclei in PDD and (ii) a shift in the expression range of the top 10% highest expressing nuclei to lower levels of SNCA expression (Fig. 4c). This was not, however, the case for PD, which maintained a similar distribution of SNCA expression to the control group, with a slight shift in the expression range of the top 10% highest expressing nuclei to higher levels of SNCA expression. The absence of a population of cells expressing higher levels of SNCA suggests that variability in SNCA expression within control ranges may contribute to the selective vulnerability of subpopulations of excitatory neurons to Lewy body pathology.

PD-associated pathways were leveraged from a recent study identifying 46 pathways implicated in PD through pathway-specific polygenic risk score and rare variant burden analyses [12]. Based on case–control comparisons, we found that pathways that have been genetically associated with PD causation (such as terms related to synaptic transmission and vesicle-mediated transport) were dysregulated in all major cell types, with the exception of vascular cells, wherein only 3 pathways were implicated (Fig. 4d, Supplementary Fig. 10, Supplementary Table 7). We noted that the number of dysregulated pathways tended to increase with increasing clinical disease severity (i.e. PD < PDD < DLB) in excitatory neurons and glia, but not inhibitory neurons and vascular cells, supporting the notion of a disease spectrum. In general, fewer pathways were dysregulated in inhibitory neurons, with 12 of 46 pathways dysregulated in at least one case–control comparison, as compared with excitatory neurons, astrocytes and oligodendrocytes (23–27 of 46 pathways).

Differentially expressed genes in glia enrich for heritability of PD age of onset and risk

To identify cell types through which common genetic variants associated with PD risk and dementia may be acting, we used Hi–C-coupled Multi-marker Analysis of GenoMic Annotation (H-MAGMA) [96] and stratified LD score regression (sLDSC) [41]. As age of PD onset is correlated with clinical progression [34, 58, 85], and there is a significant negative genetic correlation between the GWAS for PD age of onset (AOO) and PD risk [15], we included both GWASs in our analysis. Furthermore, given the potential cooccurrence of Alzheimer’s disease (AD) pathology in the Lewy body dementias, we used a recent late-onset AD GWAS [54].

Genetic association analyses with H-MAGMA and sLDSC were run with two sets of annotations: (i) the top 10% most cell-type-specific genes from each disease group and (ii) cell-type-specific DE genes (|log2(fold change)|> log2(1.5), FDR < 0.05). The latter were tested on the basis that DE genes better capture gene expression signatures representative of a given disease state. Using the top 10% most cell-type-specific genes, we observed a significant association between AD genetic risk and genes highly expressed in microglia derived from control, PD and PDD groups (control, FDRLDSC = 0.038; PD, FDRLDSC = 0.019; PDD, FDRLDSC = 0.035; Fig. 5a; Supplementary Table 8), replicating previous literature [5, 22, 54]. Furthermore, we observed a significant association between genetic determinants of PD age of onset and genes highly expressed in OPCs derived from the DLB group (FDRHMAGMA = 0.022) and PD genetic risk and genes highly expressed in oligodendrocytes (a cell type of increasing interest to the PD field [5, 22]) derived from the control group (FDRHMAGMA = 0.013).

Fig. 5

Genetic associations with top 10% most cell-type-specific genes and cell-type-specific differentially expressed genes. Genetic associations using a top 10% most cell-type-specific genes in each disease group and b cell-type-specific differentially expressed genes in disease comparisons with controls. Two methods were used to identify associations: Hi–C-coupled MAGMA (H-MAGMA) and stratified LD score regression (sLDSC). The heatmap is coloured by degree of significance with both or either method, with * and ** indicating nominal significance (unadjusted p value < 0.05) or significance (FDR-corrected p value < 0.05; corrected for number of cell types tested). Results available in Supplementary Table 8. AD Alzheimer’s disease, OPC oligodendrocyte precursor cell

Using cell-type-specific DE genes, we identified a significant association between genetic determinants of PD age of onset and genes found DE in astrocytes and OPCs from comparisons of PD with control (astrocytes, FDRLDSC = 0.0085; OPCs, FDRLDSC = 0.0085; Fig. 5b). Splitting differentially expressed genes by their direction of effect showed that this signal was driven by up-regulated genes (Supplementary Fig. 11). In addition, we identified a nominal association using both methods between PD genetic risk and genes found DE in oligodendrocytes from comparisons of PD with control (PHMAGMA = 0.011, PLDSC = 0.041; Fig. 5b), which was driven by up-regulated genes (FDRHMAGMA = 0.013, PLDSC = 0.044; Supplementary Fig. 11). Finally, we noted that genes up-regulated in excitatory neurons from comparisons of PDD with control were significantly associated with PD genetic risk (FDRLDSC = 0.040; Supplementary Fig. 11).

Differential splicing distinguishes PDD from DLB and highlights the role of specific RNA-binding proteins

Given the limitations of single-nucleus RNA-sequencing in the detection of splicing, we applied Leafcutter to our bulk-tissue RNA-sequencing to assess differential splicing (DS) [65]. Leafcutter captures changes in local splicing events through construction of intron clusters, wherein overlapping introns are connected by the splice junction(s) they share. We identified a total of 4656 DS intron clusters in 3751 genes (FDR < 0.05, |∆PSI|≥ 0.1; Supplementary Table 9) across all pairwise comparisons, with the highest number identified in comparisons of DLB with control or PD (Supplementary Fig. 12a). Notably, between 28 and 32% of DS events were partially annotated with respect to the reference transcriptome, with splicing events including novel donor or acceptor splice sites, novel exon skip and novel combination events (Supplementary Fig. 13a, b). We were, however, able to detect these events in larger control cohorts suggesting that they represent biologically relevant splicing (Supplementary Note, Supplementary Fig. 13c, d).

DS genes showed a significant enrichment in oligodendrocytes across comparisons of all disease groups with the control group (i.e. these genes had higher expression in oligodendrocytes than expected by chance), an observation that we replicated using the same external PD case–control bulk-tissue RNA-sequencing dataset used in replication of deconvolution results (Fig. 6a, Supplementary Note, Supplementary Fig. 15a, Supplementary Table 10). In contrast, enrichments in other cell types appeared to be disease specific (Fig. 6a). For example, only genes found DS in comparisons of PD with control or DLB with PD enriched in astrocytes. Notably, as the only pairwise comparison, DS genes from DLB compared with PDD consistently enriched in all excitatory neuron annotations. Pathway enrichments were observed across 4 of 6 pairwise comparisons (no enrichments were observed in comparisons of PD or PDD with control; Supplementary Fig. 12b, Supplementary Table 11). Pathways that were shared across comparisons of DLB with control, PD and PDD, included terms related to endosomes and enzyme activity (in particular, GTPase activity), mirroring terms highlighted both by replication analyses and by pathway analysis of single-nucleus DE genes (Fig. 6b, Supplementary Note, Supplementary Fig. 12b, Supplementary Fig. 15b).

Fig. 6

Cell-type enrichments of differentially spliced genes and pathway sharing across analyses. a Enrichment of the top 100 differentially spliced genes (FDR < 0.05, |∆PSI| ≥ 0.1, with rank determined by |∆PSI|) in cell types derived from each disease group. Enrichments were determined using expression-weighted cell-type enrichment (EWCE). The x-axis denotes the disease status of the cell type in question, while the y-axis denotes the groups compared in the differential splicing analysis. Pairwise comparisons have been grouped by whether diseased individuals are compared with control individuals (Ref: control) or other diseased individuals (Ref: disease). Tiles were coloured by standard deviations (s.d.) from the mean, which indicate the distance (in s.d.) of the target list from the mean of the bootstrapped samples. Multiple test correction was performed across EWCE results using FDR. Non-significant results (FDR > 0.05) were coloured white. ***FDR < 0.001; **FDR < 0.01; *FDR < 0.05. All results available in Supplementary Table 10. b Clustering of shared pathway enrichments using genes identified across the three main analyses (represented by grey bar entitled, “Analysis”). These included: bulk-tissue differential splicing (“Bulk DS”, Supplementary Fig. 12); gene contributions to bulk-tissue gene expression PC1 (“Bulk PC”, Supplementary Fig. 6); and single-nucleus differential expression (“snRNA DEG”, Fig. 3). Pathways (in rows) from all three analyses were filtered to include only those that appear across more than one type of analysis. Pathways are ordered from highest to lowest by the number of gene sets in which they are enriched (as displayed in the bar plot on the right-hand side). Gene sets (in columns) are clustered using hierarchical clustering on the Pearson correlation between gene sets (pathways were encoded with a binary 1 for “Present” or 0 for “Absent”, represented on the plot by black and white, respectively). Gene sets derived from differential splicing (Bulk DS) were collapsed across our own dataset and the replication dataset, resulting in one gene set (column) per pairwise comparison. Likewise, gene sets derived from up- and down-regulated single-nucleus DE gene sets were collapsed across cell types (represented by the coloured bar entitled, “Cell type”), such that each cell type was represented by a single column. Pathway overlaps using pairwise comparisons between disease groups are displayed in Supplementary Fig. 16. ∆PSI delta percent spliced in, GO gene ontology, OPC oligodendrocyte precursor cell

Visualisation of pathway sharing across gene sets derived from the three analyses (bulk-tissue differential splicing, gene contributions to bulk-tissue gene expression PC1 and single-nucleus differential expression) demonstrated limited sharing between the two bulk-tissue analyses (the exceptions being “presynapse”, “transport vesicle”, “coated vesicle”, and “endosome membrane”; Fig. 6b; Supplementary Fig. 16). Notably, pathway analysis of DS genes from DLB compared with PDD implicated a much wider breadth of pathways compared with pathway analysis of single-nucleus DE genes from the same comparison, and indeed, no pathways overlapped between the two analyses in this pairwise comparison (Supplementary Fig. 16). This observation suggests that differences between PDD and DLB are not sufficiently captured by consideration of gene expression alone.

Patterns of pathway sharing between each of the bulk-tissue analyses and single-nucleus differential expression highlighted highly shared terms related to synaptic function, unfolded protein binding, and vesicle transport. Of note, RNA splicing was (i) jointly implicated by differential splicing and single-nucleus differential expression derived from excitatory neurons, oligodendrocytes, astrocytes and microglia in comparisons of DLB with control and (ii) separately implicated by single-nucleus differential expression derived from excitatory neurons and oligodendrocytes in comparisons of PDD with control (Fig. 6b). Together with the abundant differential splicing observed, these results indicated that dysregulation of splicing factors may play a role in the pathogenesis of LBDs.

To further investigate this observation, we used a catalogue of known RNA-binding protein (RBP) binding motifs from the ATtRACT database [45], and defined introns by their proximal intronic regions (the 50 nt of an exon and 500 nt of an intron flanking the 5′ and 3′ splice sites), which are an important region for splicing regulation [80]. Proximal intronic regions from DS introns were compared with non-DS introns across each pairwise comparison, identifying a total of 4 RBP binding motifs with a significant enrichment in DS proximal intronic regions from at least one pairwise comparison (Supplementary Table 12). Among these was the consensus sequence GGGGGGG in DS proximal intronic regions from PDD comparisons with control (Bonferroni-adjusted p value = 0.000601; Supplementary Table 12). This sequence is targeted by 17 RBPs from the ATtRACT database (including several members of the hnRNP family, such as HNRNPC and FUS), as well as RBPs not included in the database, such as RBM25 [25, 36]. Notably, RBM25 was found DS across comparisons of PDD with control in our own dataset and the replication dataset (in-house, clu_26788, FDR-adjusted p value = 0.00653; SRP058181, clu_12260, FDR-adjusted p value = 0.0499; Supplementary Table 9). Furthermore, the consensus sequence GAAGGAA, targeted by HNRNPM, was enriched in DS proximal intronic regions from comparisons of DLB with control and PD (Bonferroni-adjusted p values, control vs DLB = 0.0141, PD vs DLB = 0.00133). Finally, two consensus sequences, CUGGAUU and CUAACCCUAA targeted by SRSF9 and PCBP2, respectively, were enriched in DS proximal intronic regions from comparisons of DLB with PDD (Bonferroni-adjusted p values, CUGGAUU = 0.000958, CUAACCCUAA = 0.0174). Of note, SRp30c (encoded by SRSF9) has been shown to interact with hTRA2-β (encoded by TRA2B) [110, 113], which targets the consensus sequence AAGAAGAAGAA, which we also found to be nominally enriched in DS proximal intronic regions from comparisons of DLB with PDD (Bonferroni-adjusted p value = 0.0865).

Overall, these results highlighted (i) the abundant levels of alternative splicing, particularly in PDD and DLB, with evidence to suggest that certain splicing factors may play a role in orchestrating these disease-related splicing changes and (ii) that differential splicing, particularly in comparisons of DLB with PDD, captures additional features of disease-related perturbations, which were not captured by single-nucleus differential gene expression.

Discussion

Here, we applied paired bulk-tissue and single-nucleus RNA-sequencing to transcriptomically profile PD, PDD and DLB. Using this approach, we (i) found transcriptional differences relative to controls for multiple cell types across the LBDs, with PDD and DLB more severely affected than PD; (ii) observed high levels of alternative splicing, particularly in PDD and DLB; and (iii) identified splicing factors, with links to other dementia-related neurodegenerative diseases, that may coordinate these disease-related splicing changes. Together, these results highlight transcriptomic commonalities and distinctions between the LBDs, which can be used to inform our understanding of the relationship between these three clinical disorders.

Existing transcriptomic studies of the LBDs have relied on bulk-tissue analyses and profiled each disease separately, limiting our understanding of the molecular landscape of these diseases individually and in relation to one another. In addition, few initiatives have addressed genome-wide assessment of splicing in this context, despite studies implicating alternative splicing as a disease mechanism in monogenic and sporadic forms of PD [31, 66], and complex disease, in general [64]. Using multiple sequencing and analytic approaches, our analyses had the potential to identify differences between the LBDs attributable to changes in cell-type proportions, cell-type-specific gene expression and bulk-tissue splicing. While we found that increases in microglial and vascular cell-type proportions were a feature of LBDs, these increases did not distinguish among the LBDs. Importantly, the observed microglial increase was consistent with results from: (i) an RNA-sequencing-based study of PD modelling cellular composition in the frontal cortex, where microglial and oligodendrocyte marker gene profiles were increased in PD compared to control [79] and (ii) a study of cell numbers and DNA content in LBD-affected brain regions, which showed an increased number of large-sized and all nuclei (implying gliosis) in the anterior cingulate cortex of LBD cases compared to controls [82]. In contrast to cell-type proportions, cell-type-specific differential gene expression and bulk-tissue differential splicing distinguished PD from the Lewy body dementias, with PDD and DLB demonstrating a higher degree of commonality. These results suggest that irrespective of when dementia onset occurs in the disease process it gives rise to similar end-stage, post-mortem transcriptomic signatures in the anterior cingulate cortex.

It is notable that bulk-tissue differential splicing (i) was a prominent feature of the LBDs; (ii) discriminated between PD and the Lewy body dementias; and (iii) provided evidence of relationships with other neurodegenerative diseases clinically associated with dementia. Enrichment analyses using DS genes associated with each of the three LBDs revealed shared cell-type associations, such as the differential splicing of genes highly expressed in oligodendrocytes, as well as disease-specific cell type and pathway associations. Indeed, splicing analyses highlighted pathways relating to GTPase activity and regulation across several pairwise comparisons involving DLB, perhaps due to their role in a range of cellular processes that have been implicated in PD, such as clearance of Golgi-derived vesicles through the autophagy–lysosome system, mitochondrial fission and fusion, and p38 MAPK signalling [12, 82]. RNA splicing was additionally associated with the Lewy body dementias, by both differential splicing and single-nucleus differential expression. To further investigate these observations, we assessed RBP binding motif enrichment to identify potential upstream regulators of splicing. All four significantly enriched RBP binding motifs were targeted by RBPs that have been implicated to varying degrees in neurodegenerative diseases, with HNRNPC implicated in AD [90], and FUS, HNRNPC, HNRNPM and PCBP2 associated with frontotemporal dementia (FTD) [11]. Furthermore, not only has PCBP2 (encoding hnRNP E2) been found to colocalise with TDP-43 pathology in specific pathological subtypes of FTD [57], but SRSF9 together with TRA2B are implicated in tau splicing [110]. Given that both Lewy body dementias are characterised by co-pathology [92, 99], including tau and TDP-43 pathology, we speculate whether dysregulation of splicing might be one of the drivers of this co-pathology. Further studies will be required to understand whether this is the case.

Looking at cell-type-specific differential gene expression, the most prominent difference between the LBDs was the widespread down-regulation of genes and pathways in the Lewy body dementias, as compared with PD. In genetic association analyses, these genes did not enrich for genetic determinants of PD age of onset or PD risk, suggesting that this down-regulation is a consequence of the disease process, as opposed to a cause. In contrast, up-regulated genes (identified primarily in comparisons of PD with control) enriched for genetic determinants of PD age of onset and PD risk, highlighting known (OPCs/oligodendrocytes [5, 22]) and new (astrocytes) cell types in PD pathogenesis. In fact, common to all three LBDs was the presence of transcriptional alterations across multiple cell types. While DE genes were found to be largely cell-type-specific (i.e. DE in only one cell type), these genes converged on similar pathways, with GO terms found to be perturbed across multiple cell types in a given case–control comparison. Restricting to genes and pathways genetically associated with PD (which arguably are more likely to be causal), we similarly saw multiple cell-type involvement across all three LBDs, albeit with some suggestion of a hierarchy of increasing perturbation in excitatory neurons and glia (i.e. PD < PDD < DLB). Together, these results suggest the involvement of multiple cell types in LBD pathogenesis, and potentially indicate a common regulatory response across cell types in each disease.

While we observed transcriptional alterations in multiple cell types, some cell types, such as excitatory neurons and oligodendrocytes, were more strongly impacted than others (most notably, excitatory neurons), implying some degree of selective vulnerability. In support of this observation, expression of SNCA (encoding α-synuclein, the major component of Lewy bodies [102]) in excitatory neurons from the Lewy body dementias, as compared with the control group, was marked by a decrease in the proportion of SNCA-expressing nuclei in PDD and a shift in the expression range of the top 10% highest expressing nuclei to lower values. While we recognise that this is an observational study, it is tempting to speculate that (i) variability in physiological levels of SNCA may impact on pathogenesis, an area of research that has received far less attention as compared with increased SNCA dosage [26, 66, 97, 100] and (ii) that the absence of cells expressing high physiological levels of SNCA may contribute to the selective vulnerability of subpopulations of excitatory neurons to Lewy body pathology.

There are several limitations to this work. Some of these, including the use of post-mortem tissue and the subsequent inability to distinguish differences that arise early in the disease course from those that arise later, are natural limitations. Others, however, emphasise key areas for future work; the most important are the study of one brain region in diseases that gradually affect multiple brain regions and the small size of the cohort used. Where possible, we attempted to validate results in larger independent control and case–control studies, but larger studies covering more brain regions will be needed in the continuing assessment of the LBDs.

Among technological limitations, a known issue in single-nucleus RNA-sequencing is the depletion of transcripts that preferentially enrich in the cytoplasmic compartment, such as transcripts that localise to neuronal dendrites [10] and signatures of microglial activation [104]. This limitation has implications both for differential gene expression, but also downstream deconvolution and indeed, the use of single-nucleus RNA-sequencing as a reference was found to decrease the performance of three deconvolution algorithms (including Scaden) on post-mortem human brain data [76]. This limitation stresses the importance of relating cell types defined by single-nucleus RNA-sequencing back to their spatial phenotypes, a process for which the emerging field of spatial transcriptomics will be instrumental in resolving [70]. Our results provide clear hypotheses to test using spatial transcriptomics both for cell-type-specific DE analysis and analysis of differential cell-type proportions.

Among methodological limitations, we recognise that RBP binding motif enrichment oversimplifies the biology of RBPs. A common feature of RBPs is the presence of multiple RNA-binding domains, which are thought to interact with repeating motifs spaced apart on pre-mRNA transcripts [36, 43]; this feature is not captured in the current analysis. Similarly, our analyses do not account for sequence context [36] (e.g. flanking nucleotide composition, repeated motifs, RNA structure) and thus cannot distinguish between RBPs that bind similar motifs. Developing tools that could address this in silico represents an opportunity to identify additional regulators of splicing in the LBDs.

In summary, our comprehensive transcriptomic analysis of all three LBDs highlights the complex, multi-cell-type transcriptional response to Lewy body pathology and LBD co-pathologies. Furthermore, it identifies post-mortem molecular signatures in the anterior cingulate cortex that distinguish PD from the two Lewy body dementias, such as perturbation of RNA splicing, a mechanism linked to several dementia-related neurodegenerative diseases. Together, these findings have important implications for the design of RNA-targeted therapies for these diseases and highlight a potential molecular “window” of therapeutic opportunity between the initial onset of PD and subsequent development of Lewy body dementia.

Materials and methods

Sample selection

Individuals with clinical parkinsonism and/or dementia with Lewy bodies (DLB) and pathologically confirmed PD were obtained from the Parkinson’s UK Tissue Bank. Clinical assessment of individuals was carried out on clinical notes collated retrospectively using records from movement disorder neurologists, neurosurgeons, psychiatrists, geriatricians, PD nurse specialists and general practitioners. Clinical parkinsonism was defined using the current MDS task force criteria [88] and Lewy body dementia by the most recent clinical diagnostic criteria for PDD and DLB [38, 74]. The 1-year rule, alongside positive clinical features for DLB (spontaneous parkinsonism, REM sleep behaviour disorder, fluctuating cognition and complex visual hallucinations) were used to separate individuals with PDD and DLB. Pathologic assessment was performed on representative tissue sections from recommended brain regions in the Braak α-synuclein [20] and Braak tau [19] staging systems as part of the routine diagnostic process for the Parkinson’s UK Tissue Bank. A maximum Braak tau stage of 3 was used to filter out individuals with excessive Alzheimer’s pathology, thus ensuring that dementia in these individuals arose from α-synucleinopathy. PD without cognitive impairment was defined either by (i) a lack of evidence of positive cognitive features, such as memory impairment, executive dysfunction and visuo-spatial dysfunction in retrospective clinical case notes or (ii) where positive cognitive features were reported present, cognitive impairment was ruled out based on objective cognitive testing or positive cognitive features were proven to be adverse effects of medication. In addition, where possible, individuals were selected based on a post-mortem interval of less than 24 h to ensure optimal tissue quality for nuclear extraction. In total, 7 PD, 7 PDD and 7 DLB individuals were selected, matched where possible for demographic and pathologic factors, along with 7 age-matched non-neurological control individuals. Control individuals were defined by a lack of clinical neurological features and no definitive pathological diagnoses. To ensure consistency, a cutoff of Braak tau stage 3 was also used for control individuals. The severity of α-synuclein pathology in the anterior cingulate was graded semi-quantitatively from 0 to 3 based on the validated scoring system from Alafuzoff et al. [6] Furthermore, Lewy pathology (i.e. Lewy bodies and Lewy neurites) was scored using the most recent LP consensus criteria [9]. For each individual, a tissue block of cortical grey matter from the anterior cingulate was sectioned at 80 µm thickness. Adjacent sections were subsequently used for bulk-tissue RNA isolation (2 sections per sample) or isolation of nuclei for single-nuclei RNA-sequencing. Clinical, pathological and sample measures for the cohort are available in Supplementary Table 1.

Isolation of nuclei

Nuclei were isolated using buffers prepared as in Krishnaswami et al. [60], including nuclei isolation medium #1 (NIM1), nuclei isolation medium #2 (NIM2), Homogenisation Buffer (HB), 29% and 50% vol/vol iodixanol dilutions. Briefly, brain tissue sections were suspended in 800 µL HB and homogenised in a pre-cooled 2 mL dounce homogeniser, with five strokes of the loose pestle, followed by 10–15 strokes with the tight pestle. The homogenate was filtered through a BD Falcon tube with a cell strainer cap (35 µm) and centrifuged at 1000g for 8 min. Thereafter, nuclei were subjected to an additional clean-up step (density gradient centrifugation), as detailed in Krishnaswami et al., albeit with centrifugation of the layered nuclei/29% iodixanol solution at 13,000g for 40 min at 4 °C. The supernatant was carefully removed, and the nuclei pellet washed with PBS buffer (PBS + 1% BSA + 0.2 U/ml RNAseIn), filtered through a BD Falcon tube with a cell strainer cap, centrifuged at 500g for 5 min at 4 °C and washed again. Nuclei were counted using an LUNA-FL Dual Fluorescence Cell Counter (Logos Biosystems, L20001) using Acridine orange dye to stain nuclei.

Nuclei encapsulation and single-nucleus RNA-sequencing data generation

All samples were processed as per 10× Genomics Chromium Single Cell Reagent Kits Protocol (chemistry: Single Cell 3′ v2). Following manufacturer’s guidelines, the samples were processed to target 10,000 nuclei per sample. Briefly, we performed 8 cycles of cDNA amplification and 14 cycles of final indexing PCR. cDNA concentrations were measured using Qubit dsDNA HS Assay Kit (ThermoFisher, Q32851), and cDNA and library preparations were assessed using the Bioanalyzer High-Sensitivity DNA Kit (Agilent, 5067-4627). All samples were pooled to equimolar concentration and sequenced together across 28 lanes on an Illumina Hi-Seq 4000.

Single-nucleus RNA-sequencing data processing

Sequenced reads were demultiplexed and processed using Cell Ranger (v 3.0.2) and thereafter mapped to the GRCh38 human reference genome using gene annotations from Ensembl v93 [35, 116]. Across each of the 28 sequenced samples, reads mapped to primary transcripts were summarised as counts. Droplets containing nuclei were distinguished from empty droplets (containing ambient RNA) using the EmptyDrops algorithm, as implemented in the R package DropletUtils (v 1.6.1) [69]. An ambient profile threshold of 300 UMI was used to determine the background RNA content of the empty droplets. Thereafter we removed nuclei with > 5% mitochondrial content and genes expressed in < 5 nuclei. Once low-quality nuclei had been filtered out, the dataset was normalised using the NormalizeData() function in Seurat (v 3.2.0) [103]. The default normalising method used by Seurat (version 3) is a global-scaling normalisation method, “LogNormalize”. The method normalises the gene expression values in each cell (n) by multiplying n by the total expression of the cell (a size factor of 10,000 for each cell is used by default) and log-transforming the result. After this normalisation step, we used Seurat’s pipeline to cluster the nuclei. First, distances were calculated between two nuclei with similar gene expression patterns using Euclidean algorithm and edges were drawn. Second, a Louvain algorithm was used to cluster the nuclei. Finally, clustering was carried out using the FindClusters() function using 30 principal components (PCs) and a resolution parameter of 2. The clustered cells were tested to remove barcodes with more than 1 nuclei encapsulated in the droplet using DoubletFinder (v 2.0.2), with the expected proportion of doublets set at ~ 7% [72].

Cell-type identification

The remaining nuclei were visualised using a non-linear dimensionality reduction algorithm known as Uniform Manifold Approximation and Projection (UMAP, v 0.1.10) [73]. We then used the Wilcoxon rank sum test (FDR < 0.05) implemented in the Seurat function FindAllMarkers() to identify genes differentially expressed in one cluster compared with all other clusters. Cell types were assigned by testing genes differential to a particular cell type for enrichment (Fisher’s exact test) for cell-type markers from two human single-cell datasets [61, 109]. Nuclei classified as endothelial cells and pericytes were merged into one class referred to as vascular cells.

A joint graph of 205,498 nuclei from across all individuals from each of their respective filtered datasets (referred to as the panel of datasets) was generated using the R package, Clustering On Network Of Samples (Conos, v 1.1.2) [13]. This was done to bring panel datasets into a common expression space accounting for technical differences between datasets, which could be used for downstream cell-type-specific differential expression analyses between disease groups. buildGraph() was used to construct a graph with parameters for nearest neighbour parameters set at k = 30, k.self = 5, in space of 30 CPCA (common principal component). The embedGraph() function was used to partition cells into 7 clusters for the 7 broad cell types.

Bulk-tissue RNA-sequencing data generation

RNA isolation was performed by the commercial company, BioXpedia A/S. Samples were lysed with QIAzol and RNA extracted using the RNeasy 96 Kit (Qiagen) with an optional on-membrane DNase treatment, as per manufacturer instructions. Samples were thereafter quantified by absorption on the QIAxpert (Qiagen) and their RNA integrity number (RIN) assessed using the Agilent 4200 Tapestation (Agilent). RIN ranged from 1.6 to 7.8, with a median of 6.5. Only samples derived from tissue-sections with a RIN ≥ 4.2 were included in downstream RNA sequencing. As a result, only 24 samples were sequenced (5 controls, 7 PD, 6 PDD and 6 DLB; Supplementary Table 1). 250 ng of total RNA was used as input for cDNA library construction with the TruSeq Stranded mRNA Sample Preparation Kit (Illumina), as per manufacturer instructions. To minimise read mis-assignment in downstream sample de-multiplexing, xGen UDI-UMI Adapters (Integrated DNA Technologies, Inc.) were used. Libraries were multiplexed on the NovaSeq S2 Flow Cell (the same 24 libraries were run across both lanes) for paired-end 100 bp sequencing on the NovaSeq 6000 Sequencing System (Illumina) to obtain an average read depth of ~ 180 M paired-end reads per sample.

Bulk-tissue RNA-sequencing data processing

Fastp (v 0.20.0), a fast all-in-one FASTQ pre-processor, was used for adapter trimming, read filtering and base correction [28]. Fastp default settings were used for quality filtering and base correction. Processed reads were mapped to the GRCh38 human reference genome via STAR (v 2.7.0a) using gene annotations from Ensembl v97 [35, 116]. Multi-sample 2-pass mapping was used, wherein two rounds of mapping were performed to improve the sensitivity of novel splice junction detection. ENCODE standard options for long RNA-seq were used, with the exception of (i) -outFilterMultimapNmax, which was set to 1, thus retaining only uniquely mapped reads and (ii) -alignSJDBoverhangMin, which was set to the STAR default of a minimum 3 bp overhang required for an annotated spliced alignment. Processed reads were also quantified with Salmon (v 0.14.1) using the mapping-based mode, with sequence-specific, fragment GC-content and positional bias correction options enabled (-seqBias, -gcBias, -posBias) [86]. A decoy-aware transcriptome file based on GRCh38 and Ensembl v97 was generated using MashMap2 (v 2.0) [53] and used as a reference together with the appropriate option for the sequencing library type (-libType ISF). The R package tximport (v 1.14.2) was used to transform Salmon transcript-level abundance estimates to gene-level abundance estimates [101]. Genes found to overlap ENCODE blacklist regions were removed from downstream analyses (“Key resources”) [7]. Pre-alignment quality control metrics were generated using Fastp and FastQC (v 0.11.8) [8], and post-alignment quality control metrics using RSeQC (v 2.6.4) [111]. Pipeline source code can be found in https://github.com/RHReynolds/RNAseqProcessing.

Processing of PD case–control replication dataset

Replication of several downstream bulk-tissue RNA-sequencing analyses were performed using a PD case–control bulk-tissue RNA-sequencing dataset provided by Dumitriu et al. [37] and processed for re-use by recount2 [32]. The dataset was accessed via recount2 (recount accession ID: SRP058181). The original study contained RNA-sequencing of prefrontal cortical samples (Brodmann Area 9) derived from 44 control individuals and 29 individuals with PD. Paired-end 101-bp sequencing was applied to each sample, with a mean depth of 83.3 million read pairs per sample. All samples were of a reasonably high quality, with RIN values ranging from 5.8 to 9.1 and a median of 7.6. Accessed samples were checked for any mismatch between the reported sex of brain donors and the sex as determined by the expression of sex-specific genes (XIST and DDX3Y). As a result, one control sample was removed (recount sample ID: SRR2015746; study sample ID: C0061); the sample was reported to be male, but notable expression of XIST was observed. Furthermore, as sample demographics from the original study included whether PD patients were diagnosed with dementia, the 29 PD cases were split into those with and without dementia (PD, n = 18; PDD, n = 11).

Deconvolution

Cell-type proportions in bulk-tissue RNA-sequencing samples were estimated using Scaden (v 0.9.2), a deep-learning-based deconvolution algorithm [76]. Unlike linear-regression-based deconvolution algorithms, Scaden does not require cell-type-specific gene expression profiles. Instead, Scaden trains on artificial bulk-tissue RNA-sequencing samples simulated from tissue-specific single-cell RNA-sequencing data, after which the model is used to predict cell-type proportions from real bulk-tissue RNA-sequencing samples. In this study, training data was generated separately for each individual with paired single-nucleus RNA- and bulk-tissue RNA-sequencing, allowing Scaden to capture cross-subject heterogeneity. This yielded a total of 24,000 artificial bulk-tissue RNA-sequencing samples (1000 samples per individual). Prior to generation of training data, single-nucleus RNA-sequencing counts per cell were normalised using the total counts over all genes, ensuring that every cell had the same total count after normalisation. Thereafter, artificial bulk-tissue RNA-sequencing samples were simulated using the Scaden bulk_simulation.py script, which sub-samples cells from input single-nucleus RNA-sequencing data and then aggregates expression across sub-sampled cells. Here, 1000 cells were used per simulated sample. Artificial bulk-tissue RNA-sequencing samples were combined and stored in a h5ad file, using the Scaden create_h5ad_file.py script. To ensure generated training data and bulk-tissue RNA-sequencing samples (in the form of counts normalised by library size) for prediction shared the same features (genes) and feature scale, both datasets were pre-processed with scaden process (the two datasets shared a total of 13,191 genes following processing). Following this, each of the three Scaden ensemble models was independently trained (scaden train) for 5000 steps, as recommended by the developers to prevent overfitting, using the default values for batch size and learning rate [76]. Finally, predictions for cell-type proportions were made with scaden predict.

Replication of predicted cell-type proportions was performed using a second independent PD case–control dataset accessed from recount2 (see “Processing of PD case–control replication dataset”). As the Scaden algorithm requires that training data and prediction data have a perfect overlap of features, it was necessary to re-perform pre-processing with scaden process (using library-normalised counts from the replication dataset; the two datasets shared a total of 14,094 genes following processing) and to train a new model (using the same parameters as previously). In both datasets, significant differences in cell-type proportions between disease groups were a two-sided Wilcoxon rank sum test, with FDR-correction for multiple testing.

Bulk-tissue RNA-sequencing covariate selection

Sources of variation in bulk-tissue RNA-sequencing data were identified using principal component analysis (PCA) performed on gene-level expression filtered to include only genes with count > 0 in all samples (28,692 genes) and transformed with DESeq2’s vst(), which applies a variance stabilising transformation. RIN and age of death were significantly correlated with the first and second PC, respectively. Furthermore, cell-type proportions for excitatory and inhibitory neurons, microglia and astrocytes were significantly correlated with the first, third and fourth PC, respectively. Thus, the final model for differential expression and splicing (referred to as the “cell-type- and covariate-corrected” model) consisted of the disease group and the top 4 PCs (which collectively explained 52.6% of the total variance).

To explore the effect of accounting for cell-type proportions, vst-transformed gene expression was batch-corrected using the final “cell-type- and covariate-corrected” model or a minimised “covariate-corrected” model consisting of disease group, age of death, RIN and sex. Samples were thereafter plotted by their first two principal components to determine how well disease groups separated (Supplementary Fig. 6). Batch correction was performed using the removeBatchEffect() function from the R package, limma (v 3.42.2) [91]. Prior to correction, covariates to be used in the model were scaled to ensure that variables that are measured on different scales (e.g. age of death vs RIN) are comparable.

As in the original study [37], the final model for the replication dataset (see “Processing of PD case–control replication dataset”) included disease group and the covariates age of death, RIN and post-mortem interval (PMI). In addition, cell-type proportions for all cell types were included in the final model, as these were significantly correlated with several of the top 8 PCs.

Differential gene expression

Single-nucleus RNA-sequencing

We used Model-based Analysis of Single-cell Transcriptomics (MAST, v 1.12.0), a method specifically designed to carry out differential expression analysis, on our single-nucleus RNA-sequencing data [40]. MAST is a two-part, generalised linear model. The first part of the model uses logistic regression to model whether a gene is expressed i.e. the discrete rate of expression of each gene over the background of other transcripts. The second part of the model models the level of expression (conditional on whether a gene is expressed in a cell) using a Gaussian linear model. Information from both parts of the model are combined to model changes in gene expression levels and with control for multiple sources of variation, such as cell–cell variation. MAST also models the cellular detection rate, which is defined as the fraction of genes that are detectably expressed in each cell. The cellular detection rate acts as a substitution for both technical and biological factors, such as dropout, cell volume and other extrinsic factors that could influence gene expression. Controlling for the cellular detection rate improves the sensitivity (true positive rate) and specificity (true negative rate) of MAST in the presence of confounding between the cellular detection rate and true biological signals.

To perform differential expression, cell-type-specific nuclei from each of the 28 filtered sample count matrices (see “Single-nucleus RNA-sequencing data processing”) were merged to create 7 cell-type count matrices. Genes that were expressed in ≤ 3 nuclei were removed from the analysis. Following this, differential expression analysis was performed separately for each cell type, across all pairwise combinations of the disease groups (n = 6). A likelihood ratio test was used, with age of death, post-mortem interval (PMI), and sex included as covariates. Genes with FDR < 0.05 and absolute fold-change > 1.5 were considered significant.

Bulk-tissue RNA-sequencing

Bulk-tissue differential gene expression was assessed using the DESeq2 R package (v 1.26.0) and gene-level expression filtered to include only genes with count > 0 in all samples (28,692 genes) [68]. With one exception (the maximum number of iterations allowed for convergence, maxit = 1000), default parameters were used, including the default Wald test of significance. Differentially expressed genes were identified in a pairwise manner, controlling for covariates identified using gene-level expression (see “Bulk-tissue RNA-sequencing covariate selection”). Multiple testing was performed by FDR-correction, with a cutoff of FDR < 0.05 applied for significance.

Differential splicing analysis

Differential splicing was assessed using Leafcutter (v 0.2.8), which detects splicing variation using sequencing reads with a gapped alignment to the genome (here, termed junction reads) [65]. Junction reads, which are presumed to represent intron excision events, are used to quantify intron usage across samples without any reliance on existing reference annotation. Importantly, Leafcutter does not estimate isoform abundance or exon inclusion levels, but rather captures changes in local splicing events through construction of intron clusters, wherein overlapping introns are connected by the splice junction(s) they share. As input, splice junctions outputted by STAR (SJ.out.tab) were first filtered to remove any regions that overlapped ENCODE blacklist regions (“Key resources”) [7] and thereafter converted to the .junc files used by Leafcutter for intron clustering. The conversion was performed using custom R code (convert_STAR_SJ_to_junc() in https://github.com/RHReynolds/RNAseqProcessing). Intron clusters were defined using Leafcutter’s leafcutter_cluster.py with thresholds ensuring the removal of: (i) introns supported by < 30 junction reads across all 24 samples or < 0.1% of the total number of junction read counts for the entire cluster and (ii) introns of more than 1 Mb. This yielded a total of 43,544 clusters encompassing 152,298 introns that were used for further analysis. Differentially spliced (DS) clusters were identified in a pairwise manner, controlling for covariates identified using gene-level expression (see “Bulk-tissue RNA-sequencing covariate selection”), and annotated to genes using exon files generated from GRch38 Ensembl v97 (with the Leafcutter helper script gtf_to_exons.R). As per Leafcutter default filters, only introns detected in ≥ 5 samples were tested and an intron cluster was only tested if detected in ≥ 3 individuals in each comparison group with an overall coverage of ≥ 20 junction reads. p values were FDR-corrected for multiple testing and an intron cluster and its overlapping gene were considered differentially spliced if (i) FDR < 0.05 and (ii) the intron cluster contained at least one intron with an absolute delta percent-spliced-in value (|∆PSI|) ≥ 0.1. The latter filter was applied to improve the specificity of Leafcutter [106].

Annotation of differential splicing events

Introns within intron clusters were annotated using annotate_junc_ref() from the R package Detecting Aberrant Splicing Events from RNA-sequencing (dasper, v 1.1.4) [117], which categorises junctions based on (i) whether the junction is present within the entire set of annotated introns or (ii) whether both, one of, or neither the donor and acceptor splice site precisely overlap the boundary of a known exon. For both checks, Ensembl v97 was used. When defining and clustering introns, leafcutter_cluster.py adds 1 bp to the end of a junction read; thus, to ensure optimal mapping to reference annotation, 1 bp was removed from all intron ends prior to use of annotate_junc_ref() using custom code (convert_leafcutter.R from https://github.com/RHReynolds/LBD-seq-bulk-analyses). Junctions (and the introns they represent) were then classified into one of the following categories: annotated, novel exon skip, novel combination, novel acceptor, novel donor, ambiguous gene and unannotated (“none”) (Supplementary Fig. 13). Annotated junctions are those that match the boundaries of an existing intron. Unannotated junctions have neither end overlapping a known exon. Novel acceptors and novel donors are junctions, where one end (acceptor or donor) matches the boundary of a known exon. Novel exon skip and novel combination junctions have both ends overlapping known exon boundaries, which are not part of the set of annotated introns. They are distinguished by whether their start or end overlaps exons derived from the same transcript. That is, for an event to be a novel exon skip, both the start and end must overlap an exon contained in the same transcript, whereas to be a novel combination, the start and end overlap exons are from different transcripts. Junctions that mapped to more than one gene (“ambiguous gene”) were not considered in downstream analyses.

Gene set enrichment

Functional enrichment of cell-type-specific differentially expressed genes

Functional term enrichment analysis for cell-type-specific differentially expressed genes from each pairwise comparison was performed using the overrepresentation analysis module from the R package implementation of WEB-based Gene SeT AnaLysis Toolkit (WebGestaltR, v 0.4.4) [67]. Two separate analyses were performed using (i) only non-redundant Gene Ontology (GO) terms (which are generated by selecting the most general terms in each branch of the GO directed acyclic graph structure from all terms with 20–500 genes) and (ii) 46 biological pathways associated with PD risk in a large-scale pathway-specific polygenic risk analysis [12]. For both analyses, default values for WebGestalt parameters were used, which include a minimum and maximum overlap of 10 and 500, respectively. FDR-correction for multiple testing was performed, and significant pathways were those with FDR < 0.05.

Functional enrichment of differentially spliced genes

Gene set enrichment for GO terms was performed using enrichGO() and clusterCompare() from clusterProfiler (v 3.14.3), which permit GO enrichment analysis (based on a hypergeometric distribution) and comparison across multiple gene lists [115]. Two separate analyses were run using (i) all differentially spliced genes (FDR < 0.05, |∆PSI|> = 0.1) across each pairwise comparison in the discovery dataset and (ii) genes overlapping validated intron clusters with ≥ 1 intron that shared the same direction of effect. In both analyses, default parameters were used; these included FDR-correction for multiple testing and filtering for terms with FDR < 0.05.

Functional enrichment of genes associated with bulk-tissue gene expression principal components

Genes contributing to PC1, following batch correction of cell-type proportions (as described in “Bulk-tissue RNA-sequencing covariate selection”), were extracted using get_pva_var() from the R package, factoextra (v 1.0.7). The top 100 genes contributing to gene-expression-derived PC1 were used for gene set enrichment with enrichGO() from clusterProfiler [115]. Default parameters were used, which included FDR-correction for multiple testing and filtering for terms with FDR < 0.05.

Visualisation of GO term overlaps between analyses

Overlapping GO-derived pathway enrichments from each of the three analyses (i.e. single-nucleus differential expression, bulk-tissue differential splicing, and gene expression contributions to bulk-tissue PC1) were visualised using the ComplexHeatmap R package (v 2.7.7) [47]. Pathways from all three analyses were filtered to include only those that were shared across more than one type of analysis. Pathways were encoded by a binary 1 and 0 for present and absent, respectively, permitting clustering of gene sets by Pearson correlation. Gene sets derived from differential splicing were collapsed across our own dataset and the replication dataset, resulting in one gene set per pairwise comparison. Likewise, gene sets derived from up- and down-regulated single-nucleus DE gene sets were collapsed across cell types, resulting in 7 gene sets per pairwise comparison.

Reduction of GO terms using semantic similarity

To reduce redundancy across GO-derived pathway enrichment analyses derived from various analyses (i.e. single-nucleus differential expression, bulk-tissue differential splicing, genes contributing to bulk-tissue PC1), two steps were taken. First, GO terms were filtered to exclude terms with ≥ 20 genes or ≤ 2000 genes. Second, semantic similarity of all enriched GO terms was calculated using mgoSim() from the GOSemSim R package (v 2.17.1) [114] and a graph-based measure of semantic similarity (measure = “Wang”) [108]. Thereafter, reduceSimMatrix() from the rrvgo R package (v 1.1.4) was used to reduce terms [94]. This function reduces terms by generating a distance matrix from the semantic similarity scores, which is hierarchically clustered using complete linkage (a “bottom-up” clustering approach). Both steps were combined into the function go_reduce(), available at: https://github.com/RHReynolds/rutils. The hierarchical tree was then cut at a threshold of 0.9 (leading to fewer groups), and the term with the highest semantic similarity score was used to represent each group of terms. This reduction was performed separately for each of the three analyses.

Cell-type enrichment of differentially spliced genes

Expression-weighted cell-type enrichment (v 0.99.2) was used to determine whether differentially spliced genes demonstrate higher expression in certain cell types than would be expected by chance [98]. EWCE requires two inputs: a gene list and gene cell-type specificity values derived from single-cell/nucleus data (here, termed a specificity matrix). Two sets of gene lists were run. The first set of gene lists included the top 100 differentially spliced genes (FDR < 0.05, |∆PSI|> = 0.1, ranked by p value) across each pairwise comparison in the discovery dataset. In the case, where a gene had multiple significant intron clusters, the most significant cluster with the highest |∆PSI| was used for ranking. The second set of gene lists included genes overlapping validated intron clusters with ≥ 1 intron that shared the same direction of effect. Both sets of gene lists were run together with gene cell-type specificity values separately derived from each disease group (i.e. control, PD, PDD and DLB); specificity matrices were generated for cell types in each disease group using the generate.cell.data() function of the EWCE package. For each combination of gene list and specificity matrix, 100,000 bootstrap replicates were used. Transcript length and GC-content biases were controlled by selecting bootstrap replicates with comparable properties to the target gene lists. Data are displayed as standard deviations from the mean, which indicate the distance of the mean expression of the target gene list from the mean expression of the bootstrap replicates.

RNA-binding protein binding motif analysis

Generating sequences

Two sets of sequences were generated per pairwise comparison. These sets included all differentially spliced introns (FDR < 0.05, |∆PSI|) and non-differentially spliced introns (FDR > 0.05), as defined by their 5′ and 3′ proximal intronic regions (500 nucleotides of proximal intron and 50 nucleotides of exon flanking the 5′ and 3′ splice sites). A 5′ or 3′ splice site could be associated with more than one intron (e.g. in the case of two introns with the same 5′ splice site, but varying 3′ splice sites), and thus could be associated with more than one |∆PSI| value. In these cases, the highest |∆PSI| was assigned to the proximal intronic region.

Enrichment of RBP binding motifs

The position weight matrices (PWMs) of RBP binding motifs in humans were collected from the ATtRACT database (v 0.99β) [45]. Motifs < 7 nucleotides in length and with a quality score of < 1 were removed to reduce false positives in the motif matches (quality score estimates the binding affinity between RBPs and binding sites). Furthermore, to remove redundancy between multiple motifs for one RBP, the longest available motif was selected. Finally, RBPs that had a median TPM of 0 in GTEx (v 8) anterior cingulate cortex samples were removed (e.g. RBMY1A1) [46]. This resulted in 82 unique PWMs, which were used to identify enrichment of RBP binding motifs. Analysis of Motif Enrichment (AME, v 5.1.1) [75] was used with default parameters (-scoring avg) to compare enrichment of RBP binding motifs between differentially spliced and non-differentially spliced proximal intronic regions. RBP binding motifs with an enrichment-optimised and Bonferroni-adjusted p < 0.05 were considered to be significantly over-represented in differentially spliced proximal intronic regions compared with non-differentially spliced proximal intronic regions.

Integration with GWAS

To test for enrichment of genetic association of a gene set to a trait we employed two orthogonal methods, Hi–C-coupled Multi-marker Analysis of GenoMic Annotation (H-MAGMA) [96] and stratified LD score regression (sLDSC) [41]. Both methods were run with two sets of annotations: (i) the top 10% most cell-type-specific genes, as determined using specificity values derived from EWCE (see “Cell-type enrichment of differentially spliced genes”) and (ii) cell-type-specific differentially expressed genes (FDR < 0.05, |log2(fold change)|> log2(1.5)). These annotations were run with 3 genome-wide association studies (GWASs), including Alzheimer’s disease (AD), Parkinson’s disease (PD) and Parkinson’s disease Age of Onset (PD AOO) (Table 1) [15, 54, 77] In both analyses, p values were FDR-corrected for the number of cell types tested.

Table 1

Summary of GWAS datasets

Disease	First author, year	N cases	N controls	PMID	References
AD	Jansen, 2019	71,880	383,378	30617256	[54]
PD—risk	Nalls, 2019 (excluding 23 and Me contributions)	33,674 (18,618 proxy cases from UK Biobank)	449,056	31701892	[77]
PD—age of onset	Blauwendraat, 2019	17,415		30957308	[15]

AD Alzheimer’s disease, PD Parkinson’s disease

H-MAGMA

Hi–C-coupled MAGMA (H-MAGMA) (v 1.08b of MAGMA [63]) was used to carry out gene-set enrichment analysis using three GWAS summary statistics. Gencode v26 (“Key resources”) was used to assign exonic SNPs and promoter SNPs, which is defined as 2 kb upstream of the transcription start site (TSS), to their target genes based on their genomic location. Chromatin interactions to exons and promoters generated from Hi–C performed on adult dorsolateral prefrontal cortex, were used to assign intergenic and intronic SNPs to their cognate genes [96]. Gene-level association statistics were computed using window coordinates of 10 kb downstream and 35 kb upstream.

sLDSC

Stratified LDSC (v 1.0.1) was used to test whether cell-type-specific DE genes or the top 10% most cell-type-specific genes contributed to the common SNP heritability of AD, PD or PD AOO [24, 42]. To ensure gene lists were sufficiently large, only gene lists with more than 20 genes were run. Gene coordinates (Ensembl v97, GRCh38) were extended by 100 kb upstream and downstream of their transcription start and end site, to capture regulatory elements that might contribute to disease heritability [42]. All annotations were constructed in a binary format (1 if the SNP was present within the annotation and 0 if not), using all SNPs with a minor allele frequency > 5%. Annotations were then added individually to the baseline model of 53 annotations provided by Finucane et al. (v 1.2, GRCh38), comprising genome-wide annotations reflecting genetic architecture. As annotations and the baseline model were mapped to GRCh38, all GWAS summary statistics were converted from GRCh37 to GRCh38 using the R implementation of the LiftOver tool, which is available from the rtracklayer package (v 1.46.0) [62]. HapMap Project Phase 3 (HapMap3) SNPs and 1000 Genomes Project Phase 3 European population SNPs were used for the regression and LD reference panels, respectively [1, 51]. The MHC region (chr6: 25,000,000–34,000,000, GRCh38) was excluded from all analyses owing to the complex and long-range LD patterns in this region. For all stratified LDSC analyses, we report a one-tailed p value (coefficient p value) based on the coefficient z-score outputted by stratified LDSC. A one-tailed test was used as we were only interested in annotation categories with a significantly positive contribution to trait heritability, conditional upon the baseline model.

Key resources

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 3224 KB)

Supplementary file2 (XLSX 52 KB)

Supplementary file3 (XLSX 22 KB)

Supplementary file4 (XLSX 16 KB)

Supplementary file5 (XLSX 20 KB)

Supplementary file6 (XLSX 2405 KB)

Supplementary file7 (XLSX 453 KB)

Supplementary file8 (XLSX 60 KB)

Supplementary file9 (XLSX 74 KB)

Supplementary file10 (XLSX 3116 KB)

Supplementary file11 (XLSX 22 KB)

Supplementary file12 (XLSX 32 KB)

Supplementary file13 (XLSX 9 KB)

Resource	Source/reference	Identifier/URL
Biological Samples
Frozen human anterior cingulate cortex samples	Parkinson’s UK Tissue Bank
Critical Commercial Assays
Chromium Single Cell 3’ Gene Expression Kit, v2	10 × Genomics	PN-120237
Qubit dsDNA HS Assay Kit	ThermoFisher	Q32851
Bioanalyzer High-Sensitivity DNA Kit	Agilent	5067-4627
QIAzol	Qiagen	79306
RNeasy 96 Kit	Qiagen	74181
TruSeq Stranded mRNA Library Prep Kit	Illumina	20020594
xGen UDI-UMI Adapters, 1–96	Integrated DNA Technologies	10005903
Deposited Data
ATtRACT database (v 0.99β)	Giudice et al., 2016 [45]	https://attract.cnic.es/index
Cell-type marker genes	Wang et al., 2018	http://resource.psychencode.org/ (DER-21_Single_cell_markergenes_UMI.xlsx)
ENCODE blacklist regions (v 2)	Amemiya et al., 2019 [7]	https://github.com/Boyle-Lab/Blacklist/blob/master/lists/hg38-blacklist.v2.bed.gz
Ensembl GRCh38 Ensembl v97	Ensembl genome browser	ftp://ftp.ensembl.org/pub/release-97/gtf/homo_sapiens/Homo_sapiens.GRCh38.97.gtf.gz
H-MAGMA: Hi-C gene-SNP pairs for adult dorsolateral prefrontal cortex	Sey et al., 2020 [96]	https://github.com/thewonlab/H-MAGMA/blob/master/Input_Files/Adult_brain.genes.annot
Gencode v26		https://www.gencodegenes.org/human/release_26lift37.html
GTEx portal (v 8)	GTEx Consortium, 2015 [46]	https://www.gtexportal.org/
LDSC baseline annotations (v 1.2)	Finucane et al., 2015 [41]	https://data.broadinstitute.org/alkesgroup/LDSCORE/
PD-associated genes	Blauwendraat et al., 2020 [16]
PD-associated pathways	Bandres-Ciga et al., 2020 [12]	https://pdgenetics.shinyapps.io/pathwaysbrowser/
Recount2	Collado-Torres et al., 2015 [32]	https://jhubiostatistics.shinyapps.io/recount/
Software and Algorithms
Analysis of Motif Enrichment (AME, v 5.1.1)	McLeay et al., 2010 [75]	http://meme-suite.org/doc/ame.html?man_type=web
Bulk-tissue RNA-sequencing pipeline		https://github.com/RHReynolds/RNAseqProcessing
Cell Ranger (v 3.0.2)	10 × Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation
clusterProfiler (v 3.14.3)	Yu et al., 2012 [115]	https://github.com/YuLab-SMU/clusterProfiler
Conos (v 1.1.2)	Barkas et al., 2019 [13]	https://github.com/kharchenkolab/conos
ComplexHeatmap (v 2.7.7)	Gu et al., 2016 [47]	https://github.com/jokergoo/ComplexHeatmap
DESeq2 (v 1.26.0)	Love et al., 2014 [68]	https://github.com/mikelove/DESeq2
Detecting Aberrant Splicing Events from RNA-sequencing (dasper, v 1.1.4)	Zhang et al., 2021 [117]	https://github.com/dzhang32/dasper
DoubletFinder (v 2.0.2)	McGinnis et al., 2019 [72]	https://github.com/chris-mcginnis-ucsf/DoubletFinder
DropletUtils (v 1.6.1)	Lun et al., 2019 [69]	https://github.com/MarioniLab/DropletUtils
EWCE (v 0.99.2)	Skene et al., 2016 [98]	https://github.com/NathanSkene/EWCE
Factoextra (v 1.0.7)		https://github.com/kassambara/factoextra
Fastp (v 0.20.0)	Chen et al., 2018 [28]	https://github.com/OpenGene/fastp
FastQC (v 0.11.8)	Andrews et al., 2010 [8]	http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
GoSemSim (v 2.17.0)	Yu et al., 2010 [114]	https://github.com/YuLab-SMU/GOSemSim
ggplot2 (v 3.3.2)		https://ggplot2.tidyverse.org/
LDSC (v 1.0.1)	Bulik-Sullivan et al., 2015 [23]	https://github.com/bulik/ldsc
Leafcutter (v 0.2.8)	Li et al., 2018 [65]	https://github.com/davidaknowles/leafcutter/
Limma (v 3.42.2)	Ritchie et al., 2015 [91]	https://github.com/cran/limma
MAGMA (v 1.0.8b)	de Leeuw et al., 2015 [63]	https://ctg.cncr.nl/software/magma
MashMap2 (v 2.0)	Jain et al., 2018 [53]	https://github.com/marbl/MashMap
MAST (v 1.12.0)	Finak et al., 2015 [40]	https://github.com/RGLab/MAST/
recount (v 1.11.8)	Collado-Torres et al., 2015 [32]	https://github.com/leekgroup/recount
rrvgo (v 1.1.4)	Sayols et al., 2020 [94]	https://ssayols.github.io/rrvgo/
RSeQC (v 2.6.4)	Wang et al., 2012 [111]	http://rseqc.sourceforge.net/
rtracklayer (v 1.46.0)	Lawrence et al., 2009 [62]	https://github.com/lawremi/rtracklayer
rutils (v 0.99.2)		https://github.com/RHReynolds/rutils
Salmon (v 0.14.1)	Patro et al., 2017 [86]	https://salmon.readthedocs.io/en/latest/index.html
Seurat (v 3.2.0)	Stuart et al. 2019 [103]	https://github.com/satijalab/seurat/
Scaden (v 0.9.2)	Menden et al., 2020 [76]	https://github.com/KevinMenden/scaden
STAR (v 2.7.0a)	Dobin et al., 2013 [35]	https://github.com/alexdobin/STAR
Tximport (v 1.14.2)	Soneson et al., 2015 [101]	https://github.com/mikelove/tximport
UMAP (v 0.1.10)	McInnes et al., 2018 [73]	https://github.com/lmcinnes/umap
WebGestaltR (v 0.4.4)	Liao et al. [67]	https://github.com/bzhanglab/WebGestaltR