Failure of diet-induced transcriptional adaptations in alpha-synuclein transgenic mice

Abstract

Nutritional influences have been discussed as potential modulators of Parkinson’s disease (PD) pathology through various epidemiological and physiological studies. In animal models, a high-fat diet (HFD) with greater intake of lipid-derived calories leads to accelerated disease onset and progression. The underlying molecular mechanisms of HFD-induced aggravated pathology, however, remain largely unclear. In this study, we aimed to further illuminate the effects of a fat-enriched diet in PD by examining the brainstem and hippocampal transcriptome of alpha-synuclein transgenic mice exposed to a life-long HFD. Investigating individual transcript isoforms, differential gene expression and co-expression clusters, we observed that transcriptional differences between wild-type (WT) and transgenic animals intensified in both regions under HFD. Both brainstem and hippocampus displayed strikingly similar transcriptomic perturbation patterns. Interestingly, expression differences resulted mainly from responses in WT animals to HFD, while these genes remained largely unchanged or were even slightly oppositely regulated by diet in transgenic animals. Genes and co-expressed gene groups exhibiting this dysregulation were linked to metabolic and mitochondrial pathways. Our findings propose the failure of metabolic adaptions as the potential explanation for accelerated disease unfolding under exposure to HFD. From the identified clusters of co-expressed genes, several candidates lend themselves to further functional investigations.

Introduction

Parkinson’s disease (PD) is a complex neurodegenerative disease with the demise of dopaminergic neurons in the substantia nigra leading to dopamine depletion, which, in turn, causes the characteristic clinical phenotype with progressive motility impairment with resting tremor and postural instability (1). While rare familial forms of the disease exist, the preponderance of PD cases occurs seemingly sporadic. Familial as well as sporadic forms share misfolded amyloid fibrils of alpha-synuclein (αSYN) that are accumulated in Lewy bodies, the pathological hallmark of the disease (2).

For αSYN, encoded by the SNCA locus, rare point mutations (3–7) and genomic multiplications have been identified in familial cases (8). In addition, several single nucleotide polymorphisms in SNCA have been associated with sporadic cases in genome-wide association as well as candidate gene studies (9–15).

Besides genetic contributions and age as key risk factor, an array of environmental factors and lifestyle choices seem to modulate the risk for PD (16). Diet is one of these factors that may—positively and negatively—impact the onset and progression of the disease (17). In humans, epidemiological studies suggest a higher risk for PD among individuals with greater intake of total animal fat (18–21), whereas other studies find no significant associations (19,22–24). These controversial results are also addressed in a recent meta-analysis that suggests high total energy intake rather than total fat intake being relevant for an increased risk of PD (25).

In animal models, the commonly used high-energy diets are so-called HFDs (26) with lipid-derived calories from saturated fat (lard), more unsaturated fat as well as increased levels of sucrose. It is important to note that an HFD is known to induce insulin resistance and to mimic diabetes mellitus type 2 in animal models, which, in humans, is associated with a more aggressive PD phenotype (27–29). Studies in toxin-induced PD models show HFD to exacerbate PD progression by exhibiting increased dopamine depletion in the substantia nigra, striatum and nigrostriatal pathway and an aggravation of vascular pathology (30–33). Furthermore, in a genetic mouse model expressing the human mutant h[A30P]αSYN (34), long-term HFD accelerates the onset of the locomotor phenotype, accompanied by earlier alpha-synucleinopathy and astrogliosis (35).

While several hypotheses have been proposed as to how HFD and/or diabetic conditions aggravate PD pathology, including increased oxidative stress, mitochondrial dysfunction and neuroinflammation (reviewed in (36)), the molecular and cellular mechanisms underlying these effects remain to be elucidated. To better understand associated pathways, we here investigated the brainstem and hippocampal transcriptome of αSYN transgenic mice (TG) exposed to a life-long HFD.

Results

HFD-dependent differential expression in brainstem and hippocampus

Analogous to the experimental paradigm used in our previous study (35), h[A30P]αSYN (TG) as well as wild-type (WT) mice were fed either a standard chow diet (SD, 3.8% total fat, 3.1 kcal/g) or an HFD (22.8% total fat, 4.6 kcal/g) from weaning till 12 months of age (Fig. 1A). As in our previous study, at 12 months of age, mice on HFD weighed about 50 g (48 – 61 g) in contrast to 35 g (33 – 41 g) on a standard diet (SD). For both diets, there was no difference in weight between TG and WT mice. Afterwards, RNA from brainstem and hippocampus was collected from six animals per group and subjected to RNA-sequencing (RNA-seq). After alignment and stringent quality control, the principal component analysis showed no outliers in any of the experimental groups (Supplementary Material, Fig. S1). We further computationally estimated the cell-type composition across samples in order to assess such confounding effects underlying differential expression. Cell type-specific expression profiles based on single-cell reference data (37,38), however, suggested great homogeneity and showed no significant compositional shifts between samples or groups (Supplementary Material, Figs S2 and S3).

Figure 1

Diet-dependent gene expression changes in brainstem and hippocampus of WT and TG mice. (A) Schematic diagram showing the experimental design and groups of WT and TG mice fed either an SD or a HFD from weaning till 12 months of age. From six animals per group, RNA was isolated from hippocampus and brainstem and subjected to RNA-sequencing. (B) Expression levels of human SNCA in nRPKM shown for brainstem (left) and hippocampus (right) as individual data points for each animal with the mean ± standard error of mean (SEM) per experimental group. (C) Expression levels of murine Snca for brainstem (left) and hippocampus (right) per animal as individual nRPKM data points with the mean ± SEM. (D) Number of DEGs for the main contrasts in brainstem between experimental groups in a 2 × 2 factorial design. Gene expression was modeled as a function of genotype, diet and their interaction. Union of DEGs indicated at the center. (E) Analogous to (D) with DEG counts for hippocampus.

After establishing this ground truth in the data, we examined in the first step the expression of endogenous and human SNCA. In TG mice, human SNCA showed a prominent and similar expression of about 700 normalized Reads Per Kilobase per Million (nRPKM) in both brain regions (Fig. 1C). Endogenous Snca expression remained largely unaffected by the transgene, but showed about 10-fold higher levels in hippocampus compared to brainstem (Fig. 1C). Hence, the relative overexpression was much stronger in brainstem, potentially contributing to the pronounced pathology observed in the region (34).

In a next step, expression changes along all primary contrasts were determined resulting in 264 differentially expressed genes (DEGs) for brainstem and 609 for hippocampus (Fig. 1D and E). Besides SNCA, few genes showed diet-independent differential expression in TG mice, among them Auh and Gabra2 in brainstem as well as Cntn3 and Zfp932 in hippocampus (Supplementary Material, Fig. S4). In contrast, the much larger fraction of DFGs was identified in both brain regions when mice were fed the HFD (Fig. 1D and E), suggesting interactions between genotype and diet. As partially captured directly through the interaction term of the statistical model for genes like Ppp1r14a and Grip2 in brainstem and hippocampus, respectively, these interactions seemingly encompassed lacking or even opposite responses in TG_HFD animals compared to WT_HFD (Supplementary Material, Fig. S4E and F).

Differential expression results mainly from genotype–diet interactions

In order to put expression changes into perspective, we visually examined the profiles of differential genes from each contrast across the groups (Supplementary Material, Fig. S5). As also apparent from the union of these DEGs, there were two dominating expression patterns that partitioned into clusters C1 and C2 (Fig. 2A and D). Interestingly, this applied to both brain regions in striking similarity. In cluster C1, genes showed lower expression in the WT_HFD group but higher in TG animals, largely irrespective of the diet (Fig. 2B and E). Nearly mirror-imaged, genes in cluster C2 showed higher expression in WT_HFD mice but lower expression in TG mice, again largely irrespective of diet, a pattern already alluded through genes identified earlier (Supplementary Material, Fig. S4E and F).

Figure 2

Genotype–diet interactions led to highly similar gene expression patterns in brainstem and hippocampus. (A) Heatmap of hierarchically clustered expression profiles (log₂ expression change relative to WT_SD) of 262 DEGs (union of brainstem DEGs indicated in Fig. 1E) across all experimental groups. Number of DEGs in clusters C1 and C2 indicated on the right. (B) Average expression profile and standard deviation of C1 and C2 across experimental groups in brainstem. (C) Enriched pathways (Gene Ontology, KEGG, Reactome) among C1 (left) and C2 (right) genes. Five most significant terms, their adjusted P-values and DEG count are shown. KG: KEGG; BP: biological process, CC: cellular component, MF: molecular function (D–F) analogues to (A–C) for hippocampus.

In brainstem, C1 and C2 genes were most significantly enriched for calcium signaling pathway and mitochondrion, respectively (Fig. 2C), whereas C1 and C2 genes in hippocampus were overrepresented for transmitter-gated ion channel activity and nervous system development, respectively (Fig. 2F).

Transcriptome perturbations through genotype–diet interactions encompass shifts in transcript isoform composition

In order to explore whether dysregulation in TG animals also extended to other aspects of the transcriptome, we next investigated alternative splicing. Detecting alternatively spliced parts in transcripts still is inherently difficult from RNA-seq data because reads are typically much shorter than the transcripts themselves. While a greater sequencing depth, larger sample size and suitable statistical methods allow mitigating this obstacle, available tools still vary considerably in accuracy and robustness of results (39). Since here, the hippocampal samples were sequenced deeper, we restricted the splicing analyses to this brain region and applied two complementary tools to identify robust effects.

In line with the gene count for differential expression, the largest number of differential splicing events identified by both tools was observed comparing TG_HFD with WT_HFD animals (Supplementary Material, Fig. S6A). Among them was Nab2, a gene for which HFD led to increased usage of exon 6 (of seven) in WT but not in TG animals (Fig. 3A, Supplementary Material, Fig. S6B). An increased usage of exon 6 agreed with a shift toward higher expression of the longer transcript isoform ENSMUST00000026469 over the shorter ENSMUST00000099157 in WT_HFD animals (Fig. 3B). Intriguingly, the shift in transcript composition without entailing a significant gene-level expression change was also observed in brainstem (Supplementary Material, Fig. S6C).

Figure 3

Genotype–diet dependent perturbations extend to transcript isoform usage. (A) Sashimi plot illustrating the transcript isoform structure and expression of Nab2 toward the 3′ end. Depicted are per-group read coverage and read count spanning junctions around exon 6. Coordinates of affected transcript isoforms at the bottom. MISO-estimated per cent spliced-in (PSI/Ψ) values with 95% confidence intervals to the right. (B) Transcript isoform-specific expression levels of Nab2 in hippocampal samples. Plotted are the mean TPM per group obtained with Salmon.

These results indicated that in both brainstem and hippocampus, interactions between genotype and diet gave rise to highly similar dysregulation patterns on a gene as well as transcript level.

Co-expression analyses point at a shared regulatory network in brainstem and hippocampus

The similarity of perturbances in both brain regions under HFD in WT and TG animals—from individual transcript isoforms to general expression patterns—led us to assume that underlying regulatory principles might be shared between brainstem and hippocampus. Co-expression analyses lend themselves to such investigations as they identify expression similarities between genes beyond comparing DEGs with little information regarding expression relations between them. In co-expression analyses, groups of genes with similar expression patterns across conditions are derived that likely function in the same pathways or regulate the same biological processes, the so-called guilt-by-association principle (40).

Here, we used weighted gene correlation network analysis (WGCNA) (41) to explore co-expression towards identifying gene modules and hub genes underlying the dysregulated expression in TG_HFD mice. As we examined different brain regions and dietary conditions, we opted for a consensus analysis that allows deriving co-expression modules across data sets. Following this approach, the consensus network partitioned the expression space into 32 modules (Supplementary Material, Fig. S7). Of those, 12 correlated significantly with the genotype in at least one brain region for either diet (Fig. 4A). Consistent with results from the differential analysis, the correlation was typically stronger under HFD and very similar overall between the brain regions (Fig. 4A). By considering the DEG ratio in the modules, M1 and M21 stood out from the rest as they captured the strongest expression changes (Fig. 4B). Although 1093 genes in module M1 were enriched for metabolic processes, 126 genes in M21 pointed toward mitochondrial pathways (Fig. 4C).

Figure 4

Co-expression analysis identifies gene modules with significant diet-dependent genotype correlations that are highly similar across brain regions. (A) Diet-dependent genotype correlations for brainstem and hippocampus with corresponding P-values for all significant modules. Cells are color-coded according to the correlation between a module’s expression and the TG genotype. Numbers in cells indicate the nominal P-value that is marked with an asterisk when still significant after FDR 0.05 correction. Module colors to the left as per WGCNA partitioning (Supplementary Material, Fig. S7) and their gene count on the right. (B) Significant co-expression modules plotted according to their relative DEG content with respect to TG_HFD/WT_HFD in hippocampus (x-axis) and brainstem (y-axis). Median module expression change color-coded per brain region with the diameter reflecting the total gene count. (C) Enriched pathways (Gene Ontology, KEGG, Reactome) for genes in modules M1 (left) and M21 (right). Five most significant terms, their adjusted P-values and DEG counts are shown. BP: biological process; RC: Reactome.

Following the reasoning of WGCNA that higher ranked genes and their relations within the modules are biologically more meaningful (42), we focused on the top 20% of genes in M1 and M21 that were differentially expressed. From the resulting 150 genes, more than one-third were linked to the SNCA interactome (Fig. 5A), indicating that the filtering steps led to relevant αSYN-related biology. Intriguingly, for all these genes, WT mice responded to HFD, whereas TG animals showed no or even opposite regulation (Fig. 5B). This failure of adaptation to HFD apparent in different brain regions might explain the previously described accelerated pathology in TG mice under long-term HFD (35).

Figure 5

Interaction network of differentially expressed hub genes around Snca. (A) Network of Snca and 57 differentially expressed as well as highly ranked (top 20%) co-expressed genes from modules M1 and M21. Genes color-coded according to their type/function. Line width indicating number of curated interactions between each gene pair. Information derived from IPA. (B) Expression changes (as log₂ fold changes per group relative to WT_SD) for 57 alphabetically sorted network genes in (A) for hippocampus and brainstem.

Discussion

As shown in our previous study, mice expressing human mutant αSYN (TG) as a model for PD have increased neuroinflammation as well as a significantly accelerated onset of neurodegeneration and terminal phenotype under a lifelong HFD (35). Towards a better molecular understanding of the interactions between genetic predisposition and diet, we here profiled the transcriptome in brainstem and hippocampus of this mouse model. Intriguingly, we observed a transcriptional adaptation of WT animals in response to an HFD that was largely missing in TG animals. Thereby, both brain regions exhibited great similarity in transcriptomic perturbations, ranging from individual transcript isoforms, to common differential genes and overall patterns of expression. Differential as well as co-expression analyses highlighted dysregulated genes that are involved in metabolic processes and mitochondria-related pathways.

The enrichment of perturbed genes for metabolic pathways agrees with epidemiological and physiological studies in PD research (43–45). Besides dietary effects that possibly modulate disease risk and influence development, disturbances in metabolic processes are also observed in the course of the disease (46).

With respect to mitochondria, our results agree with numerous findings in animal models as well as patients that suggest a pivotal role of mitochondrial dysfunction in PD pathophysiology (47). The central role of mitochondria-related processes was recognized early on as PD patients show a decreased activity of the mitochondria electron transport chain, primarily complex I. This dysfunction is thought to lead not only to increased generation of mitochondrial-derived reactive oxygen species and subsequent oxidative damage but also to the energy failure because of the inability of neurons to compensate the lack of ATP generation. Such bioenergetic perturbations are thought to be central in neurodegeneration as they are linked directly to important homeostatic processes in dopaminergic cells such as neurotransmitter release, axonal vesical transport, protein quality control and cell metabolism (46).

Our study points at a decrease of mitochondrial gene expression in the context of PD that is well reflected on protein level. Focus has, for example, been put on defective mitochondrial protein import showing reduced protein levels of the mitochondrial translocases TIM23 and TOM20 as well as nuclear-encoded mitochondrial proteins such as NDUFS3 and COX IV in post-mortem substantia nigra samples of PD patients (48). In this context. αSYN seems to play a crucial role as the aberrant αSYN–TOM20 interaction in nigrostriatal dopaminergic neurons is associated with such loss of imported mitochondrial proteins (49), thereby confirming decreased protein levels in the human disease.

Besides the relation between mitochondrial (dys-)function and PD, there are also links between mitochondria and HFD. Typically, HFDs are known for a negative impact on mitochondrial biogenesis, dynamics and function that, ultimately, lead to glucose intolerance and insulin resistance (50). However, feeding on HFD has also been shown to increase protein levels of mitochondrial genes in muscle tissue despite inducing glucose intolerance and insulin resistance (51). A number of earlier studies provide further evidence that HFD has beneficial effects on mitochondria (52–55) and promotes their biogenesis with increased protein levels of several mitochondrial genes as well as their oxidative capacity of fatty acids in skeletal muscle (56). Hence, some of the observed expression changes in WT_HFD might reflect adaptations of the metabolic regulatory network in the cell to cope with altered fatty acid levels.

In particular, we observed increased expression of several mitochondrial genes in WT_HFD animals such as Ndufa2, Cox6a1, Uqcrq. These genes were slightly decreased in TG_SD mice already and, in contrast to WT animals, showed no increase in TG animals under HFD.

Besides genes directly involved in mitochondrial function, others previously linked to SNCA also showed failed adaptation in TG_HFD mice including Stub1, Uchl1, Atp13a2/Park9 and Mif. Interestingly, increase of these proteins is typically considered protective in the context of PD. Hence, lower expression of these genes as observed in TG_HFD animals might contribute to the detrimental phenotype of these mice (35).

Stub1 encodes an E3 ubiquitin-protein ligase, also known as CHIP in human, which targets misfolded chaperone substrates toward proteasomal degradation. A beneficial role of CHIP in the pathogenesis of several neurological diseases has been shown (57). CHIP co-localizes with αSYN in Lewy bodies and inclusions, and its overexpression reduces αSYN aggregation and increases its degradation via proteasomal as well as lysosomal degradation pathways (58). Moreover, CHIP preferentially binds and degrades toxic oligomeric forms of αSYN (59) and ubiquitinates αSYN itself (60).

Uchl1, also known as PARK5 in human, encodes a ubiquitin-protein hydrolase that plays an important role in maintaining a stable pool of monoubiquitin for ubiquitination reactions. Mutations of the gene are associated with PD and its protein is found in Lewy bodies (61,62). It has been proposed that modulation of UCHL1 activity may serve as a therapeutic tool to enhance the autophagy pathway and induce clearance of aggregated αSYN (63). Here, failure of TG mice to increase Uchl1 expression under HFD agrees with the worse phenotype (35).

A similar hypothesis can be proposed for Atp13a2 or PARK9 in human, which regulates the autophagy–lysosome pathway and plays a role in lipid homeostasis. ATP13A2 promotes the secretion of exosomes as well as secretion of αSYN via exosomes (64,65) and its overexpression suppresses toxicity of αSYN.

Lastly, Mif, encoding the macrophage migration inhibitory factor, is a pluripotent pro-inflammatory cytokine. MIF reduces apoptosis and induces autophagy in an in vitro model of PD and, thereby, might mediate neuroprotective effects in PD (66).

For these example genes, the expression increase in WT_HFD animals potentially helps to cope with challenges imposed on the cell through higher energy levels. It is known that a HFD results in neurochemical adaptations including altered neurotransmission and bioenergetics in the hippocampus of rodents (67). The failure to adaptation in TG_HFD mice might be a key principle that explains the aggravated phenotype.

Some of the highlighted genes such as Stub1 and Uchl1 are linked to ubiquitination and autophagy mechanisms and point at a crucial relationship between these mechanisms as well as metabolic- or mitochondria-related pathways. The quality of mitochondria regulating numerous metabolic pathways is known to be strictly monitored to maintain cell homeostasis (68). Impaired mitochondrial quality is readily identified to eliminate damaged mitochondria, a process relying on the orchestrated crosstalk between ubiquitin signaling and autophagy. The loss of mitochondrial quality control systems is known to be associated with many types of neurodegenerative diseases including PD (68). In PD, progressive accumulation of dysfunctional mitochondria ultimately impairs cellular metabolism and causes neuronal death. As there is increasing evidence that energy metabolism plays a fundamental role in the pathomechanism of neurodegenerative diseases (69–73), it is recently becoming a potential target for preventing and treating PD (72). However, the notion of PD as a disease initiated by dysfunctions of energy metabolism just barely started (72). Our approach highlights the failure of adaptations to metabolic challenges in the context of PD and points at candidate genes that might be addressed therapeutically.

However, although these candidate genes seem plausible and the transcriptomic analyses pointed at affected pathways in the αSYN interactome, they do not allow identifying an entry point for the perturbations or deriving further causalities between genes. Such functional aspects need to be addressed by future studies. We are, however, convinced our results provide a basis to select targets for such investigations.

Taken together, our work shows that a long-term HFD leads to gene expression adaptations of several genes linked to mitochondrial and metabolic biology in the brain of WT mice. In αSYN TG animals, in contrast, there is a widespread failure of these diet-induced transcriptional adaptations, indicating that the previously observed aggravation of phenotype might be based on a failed response of related genes.

Materials and Methods

Animals and diets

TG of C57BL/6 background expressing human mutant h[A30P]αSYN under the control of the CNS neuron-specific Thy1 promotor (34,74) were maintained as a homozygous colony. WT controls were derived from the same transgenic outcross with C57BL/6 mice and maintained as a parallel colony. Only male mice were used in the study. For grouping the mice into either a standard or HFD, we used randomization within blocks representing litters. Thereby, we sought to mitigate any litter bias due to genetic effects that might confound differential expression analyses later on. From the age of 5 weeks onward till 12 months of age, homozygous (Thy1)-h[A30P]AS and WT mice were either kept on standard chow diet (SD, n = 56) (3.8% total fat, 3.1 kcal/g, ssniff R/M H Extrudat; ssniff Spezialitäten GmbH, Soest, Germany) or HFD (22.8% total fat, 4.6 kcal/g, TD.06415 Adjusted Calories Diet 45/Fat; Teklad Custom Research Diets, Harlan Laboratories, Boxmeer, The Netherlands). Groups of three to four male mice were housed in standard cages (365 × 207 × 140 mm, Typ II long) with normal light/dark cycle (12 h light/12 h dark) and free access to food and water. To avoid gene expression changes during the preparation process, WT and TG mice were sacrificed with cervical dislocation followed by head decapitation within 2 min from disturbing the home cage. Brain regions were immediately dissected on ice and snap frozen in liquid nitrogen. All animal procedures were approved by local government authorities for animal research (file references N13/16) according to the guidelines for laboratory animal care.

RNA isolation and sequencing

To focus on the same regions as before when describing the pathology in the animal model (35), the entire brainstem more precisely the medulla (oblongata and spinalis) posterior to approximately bregma −8 based on a coronal view was prepared (https://mouse.brain-map.org/experiment/thumbnails/100048576?image_type=atlas). For the hippocampus, the entire structure including all subfields was dissected. Total RNA and DNA from brainstem and hippocampus (n = 6 animals for each of the four experimental groups per brain region) were simultaneously extracted using the AllPrep DNA/RNA Mini Kit (Qiagen) using the manufacturer’s protocol. Quality was assessed with an Agilent 2100 Bioanalyzer. Samples with high RNA integrity number (>8) were selected for library construction; one WT_SD sample in brainstem failed this criterion and was not included. Using the TruSeq RNA Sample Prep Kit (Illumina), poly(A)-selected single-end sequencing libraries (75 bp read length) were generated according to the manufacturer’s instructions. To normalize the volume and amount of each brain regions per animal, 500 ng of the total RNA per sample was used. All libraries were sequenced on an Illumina HiSeq 2500 platform at a depth of 10–15 million reads each. Library preparation and sequencing were performed by the same individual using a design to minimize the batch effects.

To meet blinding strategies throughout the experimental procedure, the experimenters for brain dissection, RNA isolation, library preparation and sequencing were unaware of the animal’s group during experimentation.

Quality control, alignment and expression analysis

Read quality of RNA-seq data was assessed using FastQC (v0.11.9) (75) to identify sequencing cycles with low average quality, adaptor contamination or repetitive sequences from PCR amplification. Reads were aligned using STAR (v2.7.9a) (76) allowing gapped alignments to account for splicing against a custom-built genome composed of the Ensembl Mus musculus genome v104 and the human SNCA transgene. Normalized read counts for all genes were obtained using DESeq2 (v1.32.0) (77). Transcripts covered with <50 reads were excluded from the analysis leaving 13 309 genes in brainstem and 13 251 in hippocampus for determining differential expression in each of the primary contrasts between experimental groups.

The 2 × 2 factorial design of the experiment was captured in a generalized linear model in DESeq2 modeling expression as a function of genotype, diet and their interaction. Surrogate variable analysis (sva, v3.40.0) was used to minimize unwanted variation between samples (78). Given that differences in transcript abundances in brain tissue are often small in magnitude and in vivo RNA-seq data are deemed to be more variable (79), we set |log₂ fold-change | ≥ 0.3 and adjusted P-value ≤ 0.05 to determine DFGs.

Gene-level abundances were derived from DESeq2 as normalized read counts and used for calculating the log₂-transformed expression changes underlying expression heatmaps and k-means clustering with ratios computed relative to the mean expression in WT_SD. Thy1 and human SNCA were excluded here because of their unique and strongly biased profiles.

The sizeFactor-normalized counts provided by DESeq2 also calculated nRPKMs total reads as a measure of relative gene expression as motivated before (80). Transcript-level expression was determined using Salmon (v1.5.2, parameters: numGibbsSamples 20, seqBias, gcBias, validateMappings) (81). Transcripts per million (TPM) values obtained with Salmon were scaled (scaleInfReps) using the tximeta (v1.10.0) R package (82). In addition, JunctionSeq (v1.21.0, default parameters) (39) was used to identify alternative splicing events underlying transcript-level changes.

To identify alternative splicing events (alternative 5′ and 3′ splice site, retained intron, skipped exon, mutually exclusive exon), MISO (v0.5.4) was used for merged hippocampal samples and filtered for events with --num-total 200 --num-inc 10 --num-exc 10 --num-sum-inc-exc 20 --delta-psi 0.2 --bayes-factor 10 (83,84).

Gprofiler2 (v.0.2.1) (85) with P_FDR ≤ 0.05 was used to determine functional enrichments among gene sets against Gene Ontology, KEGG and Reactome. Thy1 and human SNCA were excluded when determining enrichments.

Interactions among genes were derived from curated data in Ingenuity Pathway Analysis (IPA, v01–20-04, Qiagen) and visualized in Cytoscape (86).

Co-expression analysis

Four sets were formed, each containing the expression data from hippocampal or brainstem samples exposed to SD or HFD, respectively. WGCNA (41) was used to establish a consensus network in order to identify gene co-expression across all four data sets. WGCNA was based on pairwise correlation between all gene pairs in each data set. As the correlation method, biweight midcorrelation (87) was used with maxPOutliers = 0.1, thereby minimizing the influence of potential outliers. Correlations were transformed in a signed hybrid similarity matrix where negative and zero correlations equal to zero, whereas positive correlations remain unchanged. To generate the network adjacency, the similarity matrix was raised to the power β = 8, the minimum value that approximated a scale-free topology in all data sets, thus suppressing low correlations. For a measure of interconnectedness, adjacency was transformed into a topological overlap measure (TOM) that is informed by the adjacency of every gene pair plus the connection strength they share with the neighboring genes. Before calculating the consensus TOM (cTOM), the individual TOMs were calibrated using a 0.95-single quantile scaling. TOMs were raised to a power such that the 95th percentile of all other data sets equaled the same quantile of the reference set (hippocampus under SD). Thus, potential bias deriving from different statistical properties were mitigated best. cTOM was created by selecting the component-wise 0th quantile of the individual TOMs, meaning that for each gene pair, the minimum TOM value across all sets was given as input. A hierarchical clustering of a TOM-based dissimilarity measure (1-cTOM) was used to define modules by applying the Dynamic Tree Cut algorithm (88). Each of these modules was summarized by its eigengene, providing a single value for a module’s expression profile. Based on the eigengene correlation matrix, final modules were derived by iteratively clustering eigengenes based on dissimilarity (given by one minus the respective correlation) and cutting the resulting dendrogram at height 0.1, causing all modules with eigengene correlation ≥ 0.9 to be merged. Finally, module eigengenes were correlated with a dichotomous genotype trait to identify modules affected by transgenic effects. A joint Bayesian-frequentistic algorithm combining the Bayes factor (BF) (89) and significance of a correlation was used to identify modules associated with the disease status. Modules with an eigengene-trait correlation with P_FDR ≤ 0.05 | BF ≥ 3 were considered significantly associated with genotype status.

Acknowledgements

We thank Marina Karakhanyan for her animal work assistance, Zinah Wassouf for her help in isolating RNA/DNA as well as Olaf Riess for advice and discussions. In addition, we would like to thank the core facility cATG for preparing the libraries and sequencing the samples. We acknowledge the support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of University of Tübingen.

Conflict of Interest statement. None declared.

Funding

This study was funded by the decipherPD transnational consortium on Epigenomics of Complex Diseases (Bundesministerium für Bildung und Forschung, grant number: 01KU1503). A.K. received a scholarship from the IZKF-Promotionskolleg, University Hospital Tübingen.

Data availability

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

Authors’ contributions

J.M.S.H. and P.J.K. initiated the study and designed the experiments with C.R. C.R. took care of animals and their feeding as well as isolation of RNA/DNA. A.K. and T.H. were responsible for computational analyses of the data. A.K., T.H. and J.M.S.H. wrote the paper. All authors read and approved the manuscript.

Author notes

Alexander Kilzheimer and Thomas Hentrich authors share equal first authorship