Brain region-specific and systemic transcriptomic alterations in a human alpha-synuclein overexpressing rat model

Overexpression of human SNCA transcripts differs between brain regions in rat

To better understand effects of SNCA overexpression in brain regions, the striatal, frontocortical, and cerebellar transcriptomes of transgenic rats (TG) overexpressing full-length human SNCA, and wild type (WT) controls were interrogated using RNA-sequencing. By including 5- and 12-month-old rats of both genotypes a particular emphasize was put on revealing longitudinal transcriptomic changes (Figure 1A). After preprocessing the raw data using well-established, community approved pipelines [24, 25], high quality of the data was ensured by stringent quality controls. Principal component analysis showed that gene expression differed strongly between brain regions (Supplementary Figure 1A), age and genotype (Supplementary Figures 1B, 3A, 5A, 7A). In line, cell type-specific single-cell reference data [26] indicate cell type composition differences between brain regions, most prominent with respect to microglia, pyramidal neurons, and oligodendrocytes. Within each brain region, experimental groups reflect great homogeneity suggesting that transcriptomic differences are unlikely to result from compositional shifts (Supplementary Figure 2).

Figure 1. Brain region-specific overexpression of endogenous and human SNCA in rats. (A) Schematic illustration showing the experimental design of 5- and 12-months-old wildtype (WT) and transgenic (TG) BAC SNCA rats. RNA was isolated from striatum, cortex and cerebellum for 5 animals per group and used for RNA-seq to identify differential expression. (B) Expression level of endogenous rat Snca (left) and human SNCA (right) as individual nRPKM data points per rat across experimental groups with mean and standard error of the mean. (C) Composition and expression level of endogenous and human SNCA transcript isoforms across experimental groups.

Since aSyn protein load differs between brain regions in this model [19], expression of endogenous and human SNCA were analysed. Endogenous Snca showed the highest expression in cortex and lowest in cerebellum and remained largely stable with respect to age and genotype. In striatum, however, endogenous Snca showed a significant age-dependent expression increase in both WT and TG animals (Figure 1B). On top of endogenous Snca levels, a region-specific expression of human SNCA was detected with highest expression in cortical samples of TG rats. Like the endogenous Snca, expression levels of human SNCA were not affected by age in cortex and cerebellum, whereas a significant age-dependent increase was observed in striatum (Figure 1B). Furthermore, expression of SNCA transcripts was analysed, revealing an intriguing addition of human transcripts on top of endogenous Snca in TG rats. Although no compositional differences of SNCA isoforms were observed, relative isoform expression differed significantly across brain regions, dominated by ENST00000508895.5 and ENST00000394991.8 transcripts (Figure 1C).

These results indicate that human SNCA is expressed in an age- and region-specific manner that follows the expression of endogenous Snca on both gene and transcript level, lending the model to learn more about aSyn biology under normal and disease-like conditions.

Striatal transcriptomic changes under SNCA overexpression largely occur already in young animals

Since the striatum is notably affected by aSyn pathology, particularly in PD [22, 27], the striatal transcriptome was analysed first. Globally, samples partitioned according to experimental conditions (Supplementary Figure 3A), and differential expression was determined with respect to age, genotype, and their interaction. Surprisingly, more differential expression was detected in younger than older TG rats, with 578 differentially expressed genes (DEGs) (369 up- and 209 downregulated) in 5-month-old and 225 DEGs (127 up- and 98 downregulated) in 12-month-old rats, of which 125 DEGs overlapped (Figure 2A and Supplementary Figure 3B). The counter-intuitive numbers of differential genes for an age-related disease also extended to transcript level (Supplementary Figure 3E) and cannot be attributed to higher inter-individual variability of the data in older animals (Supplementary Figure 3F). Also, differential splice event counts were higher in young TG rats, dominated by skipped exons, and occurring in genes associated with cytoplasm and synapse. Despite splicing events were largely distinct, there were few reoccurring genes with respect to genotype and age, respectively (Supplementary Figure 18A–18D).

Figure 2. Striatal gene expression changes under SNCA overload largely occur in young TG rats. (A) Number of differentially expressed genes (DEGs) between experimental groups in the striatum, along the genotype (WT and TG) and age axes (5 and 12 months) and their interaction with the indicated significance cut-offs. (B) Heatmap of hierarchically clustered striatal expression changes relative to WT_5m of 2776 DEGs (union of striatal DEGs shown in Figure 2A) across experimental groups. (C) Average gene expression changes and standard deviation of all striatal DEGs plotted as centroids clustered in eight groups. Numbers of DEGs shown in brackets. (D) Five most significant enriched biological processes for the DEGs in cluster 6 with indicated adjusted p-value, enrichment ratio and DEG count in brackets. (E) Expression level for Gria2 as individual nRPKM data points across experimental groups with mean and standard error of the mean.

The differential genes in 5-month-old TG rats were most significantly enriched for monovalent inorganic cation transport and sodium ion transport (Supplementary Figure 3C), whereas the DEGs in 12-month-old TG rats were neither significantly enriched for Gene Ontology terms, KEGG pathways, nor related annotations (Supplementary Figure 3D).

Trying to relate the differential entities with other models, even the striatal transcriptome of 6-months old mice overexpressing the same BAC SNCA construct showed limited overlap in differential genes (Supplementary Figure 4). Similar for RNA-seq data from putamen of PD patients [28], while indeed 83 orthologues were shared between TG rats and patients, similarity in directionality and magnitude was restricted to few genes only, including LILRA5 and TNNI3 as the most significantly up- and downregulated genes (Supplementary Figure 3G, 3H).

With respect to age, 1988 DEGs (971 up- and 1017 downregulated) were identified in WT and 1399 DEGs (482 up- and 917 downregulated) in TG rats (Figure 2A), suggesting that aging under aSyn overload differs between genotypes, which is statistically modeled by the interaction term that captured 98 DEGs (35 up- 63 downregulated) showing an interaction between age and genotype.

Borrowing reference data for striatal aging in mice [23], the expression patterns of these signature genes in the data formed four dominant clusters (Supplementary Figure 5). Cluster A2 and A3 contained genes that were either up- or downregulated specifically in TG rats, indicating an age-independent effect of the overexpressed transgene. In addition, the largest cluster (A4) contained genes showing an age-dependent upregulation largely irrespective of genotype. Furthermore, genes showing the strongest age-dependent upregulation, clustered separately in A1 because young TG rats reflected premature upregulation of those genes.

Broadening the view on these distinct expression patterns, DEGs from all striatal contrasts were merged, visualized in a heatmap and hierarchically clustered, revealing eight distinct groups (Figure 2B), summarized as centroids (Figure 2C). These clusters contained genes of analogous expression patterns described before (C4, C5 and C8). Furthermore, the largest clusters contained genes showing a premature up- (C6) and downregulation (C1 and C3) meaning they are differential in 5-month-old rats only (Supplementary Figure 6). Thus, premature alterations largely explain the high DEG count in young animals. Functionally, the 782 DEGs identified in cluster C6 were enriched for biological processes linked to synaptic functions including >modulation of chemical synaptic transmission and regulation of trans-synaptic signalling (Figure 2D), including Gria2 which encodes for a glutamate ionotropic receptor with an important role in excitatory synaptic transmission [29] (Figure 2E). In addition, premature upregulated genes were enriched for targets of the transcription factor Pax6, which showed increased activity and is known to regulate neuronal development and neuronal plasticity [30] (Supplementary Figure 7).

Taken together, these results indicate greater perturbation of the striatal transcriptome in young TG animals. At 12 months, expression differences of these genes have largely converged again such that DEGs are fewer.

Modest impact of SNCA overexpression on cerebellar transcriptome

To investigate whether transcriptomic changes occurring largely in young TG rats is specific to the striatum, analyses were extended to cortical and cerebellar samples. While the PCA of the cerebellar samples also indicates a genotype effect similar to striatum, samples separate predominantly by age (Supplementary Figure 8A). In line, differential analysis revealed 2475 DEGs (1719 up- and 756 downregulated) in WT and 1873 DEGs (1236 up- and 637 downregulated) in TG animals with respect to age, of which 1405 overlapped (Figure 3A and Supplementary Figure 8B). With respect to genotype, a similar DEG count of 367 (203 up- and 164 downregulated) and 311 (166 up- and 145 downregulated) were detected in 5- and 12-month-old TG rats, respectively (Figure 3A), with Sncg, Gstt4, and Col27a1 being the most significantly upregulated and P2rx4, Rtn4ip1, and Evc the most significantly downregulated DEGs in TG animals. Of note, the DEGs of TG rats could neither be characterized by significant enrichments for Gene Ontology terms, KEGG pathways, and related annotation data. Differential splice events were similar to striatal observations, in which skipped exons dominated all contrasted conditions, occurred in genes annotated for cytoplasm, and were otherwise largely disjunct with the exception of a few reoccurring genes with respect to age and genotype (Supplementary Figure 18E–18H).

Figure 3. Modest impact of SNCA overexpression on cerebellar transcriptome. (A) Number of DEGs between experimental groups in the cerebellum, along the genotype (WT and TG) and age axes (5 and 12 months) and their interaction with the indicated significance cut-offs. (B) Heatmap of hierarchically clustered cerebellar expression changes of 3402 DEGs relative to WT_5m (union of cerebellar DEGs shown in Figure 3A) across experimental groups. (C) Average gene expression changes and standard deviation of all cerebellar DEGs plotted as centroids clustered in eight groups. Numbers of DEGs are shown in brackets.

Due to the lack of cerebellar PD patient data, relating differential genes to the human condition remains for future studies. With respect to other models, genotype-driven DEGs partially overlapped with genes identified in cerebellum of PLP-alpha-synuclein transgenic mice, with Sncg being most significantly upregulated and Chst5 being the most significantly downregulated in both models at 5 and 2 months of age and Sorc2 and Casq1 being most significantly up- respectively downregulated in both models and 12 months of age (Supplementary Table 1).

The similar count of differential genes in the cerebellum of BAC SNCA rats at 5- and 12-months of age also extends to transcript level (Supplementary Figure 8C). Together, the similarity of DEGs with respect to age in both genotypes and the stable impact of the transgene agree with only 21 genes (7 up- and 14 downregulated) that show significant interaction between age and genotype, suggesting only modest deviation on how WT and TG rats aged (Figure 3A). Yet again, signature genes of cerebellar aging showed similar patterns that were observed for striatal reference genes (Supplementary Figure 9), most importantly a cluster of genes showing premature upregulation in TG rats (A1 in Supplementary Figure 9B).

Investigating this further, the union of all cerebellar DEGs were visualized across experimental groups (Figure 3B), revealing expression patterns that clustered into eight distinct groups (Figure 3B, 3C and Supplementary Figure 10). Similar to striatum, two clusters (C5 and C3) comprise DEGs that show age-dependent up- and downregulation, respectively, largely independent of genotype. In addition, there are two cluster (C4 and C8) which show stable genotype driven perturbations (C4 and C8). Lastly, the largest number of genes were partitioned into clusters C6 and C7, showing similar activation with age in both genotypes, but premature upregulation in young TG animals. However, despite of clear and homogenous expression profiles, these DEGs were not significant enriched for any biological processes and pathways.

It should be noted that despite stronger overexpression of the transgene in cerebellum and the large DEG count in cerebellar contrasts their effect sizes tended to be weaker compared to striatal DEGs (Supplementary Figure 8D), suggesting that SNCA overexpression has a weaker impact on the aging cerebellar transcriptome.

Myelination-linked genes were differential in cortex of BAC SNCA rats and PD patients

In a next step, cortical samples of both time points and genotypes were examined (Supplementary Figure 11A). Whereas in 5-month-old TG rats 610 DEGs (362 up- and 248 downregulated) were identified, 255 genes (169 up- and 86 downregulated) were differential in 12-month-old TG animals. Having more perturbations in young TG animals—similar to results in striatum (Figure 2)—also extended to transcript level and splice event counts (Supplementary Figures 11C, 18I) and was, again, not due to increased data variability in older animals (Supplementary Figure 11D). Also similar to striatal observations, differential splice events occurred mostly disjunct across contrasts, and were enriched for cytoplasm-linked genes (Supplementary Figure 18J, 18K).

With respect to age, 1385 DEGs (1016 up- and 369 downregulated) were detected in WT animals, whereas only 526 genes (436 up- and 90 downregulated) were differential in TG rats, with 315 DEGs overlapping (Figure 3A and Supplementary Figure 11B). As suggested by these numbers cortical aging under SNCA overexpression differed significantly from WT animals, most prominently captured in 197 DEGs (79 up- and 118 downregulated) showing a significant interaction between age and genotype, the highest number of all three examined brain regions (Figure 3A).

In addition, overall higher effect sizes were detected compared to striatum and cerebellum (Figure 4B and Supplementary Figure 8D). Clustering of this expression changes revealed six distinct groups with similar perturbation patterns as seen in the other regions (Figure 4B and Supplementary Figure 12). Also similar to striatum and cerebellum the largest number of DEGs (C4 and C6 in Figure 4B) show the strongest age-dependent upregulation and premature upregulation in young TG animals. This phenomenon was also reflected in the cortical age signature (Supplementary Figure 13). Functionally, the DEGs identified in cluster C4 and C6 were most significantly enriched for biological processes linked to synaptic functions including regulation of postsynaptic membrane potential, regulation of membrane potential, and modulation of chemical synaptic transmission (Supplementary Figure 14A), suggesting that synapse-related pathways are affected early on.

Figure 4. Cortical transcriptomic perturbations under SNCA overload include myelination-linked genes in rat and PD patients. (A) Number of DEGs between experimental groups in the cortex, along the genotype (WT and TG) and age axes (5 and 12 months) and their interaction with the indicated significance cut-offs. (B) Heatmap of hierarchically clustered cortical expression changes relative to WT_5m of 2061 DEGs (union of cortical DEGs shown in Figure 4A) across experimental groups. (C) Five most significant enriched biological processes for the 255 cortical DEGs of 12-month-old TG rats with indicated adjusted p-value, enrichment ratio, and DEG count in brackets. (D) Venn diagram comparing 255 DEGs identified in the cortex of 12-month-old TG rats and 979 DEGs identified in the cortex of PD patients [31]. (E) Scatter plot of 28 overlapping DEGs identified in the cortex of 12-month-old TG rats and in the cortex of PD patients. Oligodendrocyte associated DEGs are labelled.

Further exploration of cortical DEGs under SNCA overload revealed significant enrichment for central nervous system myelination and axon ensheathment in central nervous system in TG rats, primarily at 12 months of age (Figure 4C and Supplementary Figure 14B). These genes were also significantly enriched for genes attributed to oligodendrocytes (Supplementary Figure 14C), suggesting that myelination-associated and oligodendrocyte-specific genes are significantly perturbed in old TG rats.

To relate these finding in the animal model to human, previously published post-mortem cortical transcriptomes of PD patients were analysed analogous to rat samples [31]. Among 979 genes differentially expressed in the cortex of PD patients, 28 orthologues overlapped with DEGs identified in 12-month-old TG rats (Figure 4D). Intriguingly, 28 shared differential genes are specifically expressed in oligodendrocytes including MOG, UGT8 and MYRF, that showed increased and nearly identical expression in old TG rats and patients (Figure 4E). This is in line with previous transcriptomic analyses of cortical samples against other PD patient data [20] but did not extend to protein level (Supplementary Table 2).

Taken together, these results indicate that—compared to striatum and cerebellum—the cortex exhibits the largest effect sizes of age-dependent expression changes. Perturbations under the transgene not only occurred early on but also led to the highest number of genes with significantly different age trajectories under SNCA overload. Lastly, upregulation of myelination-associated genes in cortex of TG rats was shared between animal model and PD patients.

Core set of genes exhibits common perturbation patterns across brain regions in SNCA overexpressing rats

To investigate whether the common patterns observed in all examined brain regions are based on the same underlying genes, the overlap of DEGs across brain regions was examined (Figure 5A). Indeed, among the most significant cross-region DEGs was Sncg, encoding for gamma-synuclein, another member of the synuclein family, Ddc encoding for dopa decarboxylase as well as genes like Fam111a and Gjc3 (Supplementary Figure 15). In addition, there were several cross-region DEGs agreeing with known aSyn biology such as Orai2 encoding for a calcium channel and Gpr157, both showing age-independent upregulation, or P2rx4 a ligand-gated channel with high calcium permeability significantly downregulated across brain regions in TG animals (Supplementary Figure 15).

Figure 5. Common transcriptomic changes in SNCA overexpressing rats across brain regions and tissues. (A) Venn diagram comparing DEGs identified along the genotype axis in 5- and 12-month-old WT and TG rats between the striatum, the cortex and the cerebellum. (B) Heatmap of striatal, cortical and cerebellar expression profiles of 122 overlapping genes (shown in A) as log₂ expression change relative to striatal WT_5m per experimental group. (C) Circos plot of the spatial distribution of 122 common DEGs found in TG animals across all examined brain regions. Genes are color coded based on their up- or downregulation in at least two brain regions. (D) Venn diagram comparing 122 DEGs identified as core set of DEGs in three examined brain regions in TG animals and 1697 DEGs identified in the gut of the same rat model [32]. (E) Scatter plot of 122 cross-regional DEGs in TG animals in striatum and gut with 37 overlapping DEGs between brain regions and gut highlighted in purple. DEGs with the same regulation in brain and gut are labelled.

Together, the common genes extend to a total of 122 (Figure 5A and Supplementary Table 3). Characterizing them further with respect to pathway and gene set enrichments, and known interactions did not yield significant or conclusive results, emphasizing that these DEGs can be understood only in context of the other region-specific DEGs (Supplementary Figures 3D, 14B). Transcription factor analysis of the 122 genes revealed enrichment of targets of the transcription factor SP1, where SP1 itself was downregulated in an age-dependent manner (Supplementary Figure 16). However, gene regulatory network analysis did not point to specific signalling pathway, even under relaxed criteria (Supplementary Figure 17). Yet, the expression pattern of the core genes is very clear across regions, with largely similar dysregulation in TG animals, and in which most genes show significantly altered expression in both young and old TG animals (Figure 5B). The cross-regional DEGs were also evenly distributed across the genome, suggesting the differential genes not to result from effects in the genomic vicinity of transgene integration sites (Figure 5C).

Further support for non-random effects underlying the core gene set is given by similarities in differential splice events. While the preponderance of them occurs region-specifically, some affected genes are shared across brain regions with nearly identical impact on transcript isoform usage. These overlaps exist specifically with respect to genotype and are exemplified by candidates like Ptprn, Zfyve28, and Ifi44 (Supplementary Figure 19).

To investigate whether the core set of differential genes across brain regions are also reflected in another tissue, published transcriptome data from gut epithelial cells of the same rat model [32] were analysed analogous to brain samples. 37 DEGs in the gut overlapped with the core set of differential genes identified in brain (Figure 5D). Again, these common DEGs displayed largely similar expression across both tissues in TG animals (Figure 5E and Supplementary Figure 20), suggesting a systemic effect of SNCA overexpression.

These results indicate that although most transcriptomic changes under SNCA overexpression are brain region-specific, there is a core set of genes in the brain which share very similar perturbations under aSyn overload. Part of these cross-regional genes are similarly dysregulated in the gut, a tissue known to play an important role in synucleinopathies, particularly in PD [33].

Experimental animals and tissue preparation

Male homozygous transgenic rats overexpressing full-length human SNCA including its regulatory elements [19] as well as WT rats with the same genetic background (Sprague-Dawley) were housed in a standard environment until the age of 5 and 12 months. Experimental animals were obtained by crossing heterozygous male with heterozygous female rats. Homozygous or WT status was confirmed by genotyping with quantitative PCR using DNA from ear biopsies with primer sequences specifically for human SNCA (F: 5′ccgctcgagcggtaggaccgcttgttttagac-3′; R: 5′-ctctttccacg ccactatc-3′) and normalized to β-actin as reference (F: 5′agccatgtacgtagccatcca-3′; R: 5′-tctccggagtccatcacaatg-3′). Animals were anaesthetized, decapitated, and the brain immediately dissected on ice and snap frozen in liquid nitrogen. The tissue was stabilized with RNAlater (Qiagen) at 4° C.

RNA isolation and sequencing

The polyadenylated fraction of RNA isolated from striatum, frontal cortex and cerebellum (n = 5 animals in each of the four experimental groups per brain region) was used for RNA-seq. Total RNA, miRNA, and DNA were simultaneously extracted using the DNA/RNA/microRNA Universal Kit (Qiagen) using the manufacturer protocol. Quality was assessed with an Agilent 2100 Bioanalyzer. Samples with high RNA integrity number (RIN > 7) were selected for library construction. Using the TruSeq RNA Sample Prep Kit (Illumina) and 100 ng of total RNA for each sequencing library, poly(A) selected paired-end sequencing libraries (125 bp read length) were generated according to manufacturer’s instructions. All libraries were sequenced on an Illumina HiSeq 2500 platform at a depth of 20–35 million reads each.

Quality control, alignment and expression analysis

The Nf-core RNA-seq Pipeline (v3.13.0) [24] served as the primary tool for processing the sequencing data. To assess the read quality of the RNA-seq data FastQC (v0.12.0) was used [60]. Reads were aligned against a custom-built reference genome of the Ensembl Rattus norvegicus genome (mRatBN7.2) including the human SNCA transgene using STAR [61]. Normalized read counts were obtained with Rsubread (v2.12.3) [62]. Quality control metrics from all samples were aggregated and reviewed using MultiQC (v1.14) [63], including e.g. total number of mapped reads, read coverage, duplication rates, and read distribution across genomic features, to ensure consistency across samples. To confirm the absence of outlier samples, principal component analysis (PCA) was performed on normalized gene expression data. All samples clustered with their respective experimental groups, and no outliers were identified.

DESeq2 (v1.38.3) [64] was used for the differential gene expression analysis. Transcripts with less than 20 reads in the median of all samples were excluded from subsequent differential analysis resulting in 14,948 genes for striatum, 14,697 for cortex and 14,521 for cerebellum. Gene expression was modeled in a 2 x 2 factorial design as a function of genotype, age and their interaction Supplementary Tables 4–6. Significance thresholds for differentially expressed genes were set to a BH-adjusted p-value ≤ 0.05 and a |log2FC| ≥ 0.5. Surrogate variable analysis (v3.46.0) [65] was applied to remove unwanted variation in the data.

Transcript abundance was quantified with Salmon (v1.10.1, parameters: numGibbsSamples 20, seqBias, gcBias, validateMappings) [25] using a gentrome based on the Ensembl Rattus norvegicus genome and transcriptome including human SNCA. Transcripts per million (TPM) values were imported and scaled (scaleInfReps) with Tximeta (v1.16.1) [66]. Only transcripts with at least 50 counts in at least 5 samples were considered for the differential transcript usage analysis. For differential splicing analysis rMATs (v4.3.0) [67] was used with default parameters. To identify significant splicing events filtering of detected splice events was based on ∆PSI ≥ 0.1 (Percent Spliced In) and a BH-adjusted p-value ≤ 0.01. Sashimi plots to visualize differential splicing events were generated using ggsashimi [68].

Heatmaps and centroids were plotted as log₂ expression changes relative to the mean expression of WT_5m samples. nRPKM values (normalized Reads Per Kilobase per Million total reads) were calculated using read counts from DESeq2 to measure the relative gene expression changes [69]. For gene ontology analysis WebGestalt and gProfiler [70] were used to identify overrepresented biological processes for differentially expressed genes [71]. Gene set enrichment analysis was conducted using FGSEA [72]. Transcription factor prediction was carried out using decoupleR [73]. Gene regulatory network analysis was performed using Neko [74] in conjunction with the Omnipath interaction database [75]. Bimodal interactions were excluded from the analysis and intermediated nodes were introduced to capture potential regulatory signalling pathways by setting the indirect network pathway length to one.

Data availability

RNA-seq data files have been uploaded to GEO database and are available under the accession number GSE281984.