Gene Expression Profiling of Post Mortem midbrain of Parkinson’s Disease patients and Healthy Controls

1. Introduction

Parkinson's disease (PD) stands out as the most prevalent movement disorder and neurodegenerative disease after Alzheimer’s dementia, affecting approximately seven million people globally [1,2]. Clinically, PD is a heterogeneous condition primarily characterized by resting tremor, bradykinesia, and rigidity. Additionally, numerous non-motor symptoms, equally debilitating and often preceding motor manifestations, contribute to the complexity of the disease [3,4,5,6]. Despite considerable research efforts, encompassing both preclinical and clinical studies, PD remains incurable. The progressive loss of dopaminergic (DA) neurons in the substantia nigra pars compacta (SNpc), a crucial part of the midbrain regulating movement tone and velocity, persists without effective intervention, and the complete pathomechanisms behind this degeneration remain incompletely elucidated [1,7].

Pathophysiologically, PD has a multifactorial origin involving complex interactions between various genetic and environmental factors. Some factors affect the elderly in general, while others seem more specifically tied to this disorder. It's noteworthy that post-mortem studies consistently report that over 60% of DA neurons in the SNc are already degenerated when overt clinical signs manifest. This indicates that PD is neuropathologically evident long before its clinical onset, suggesting a period during which the human brain can compensate for dopaminergic loss until reaching a "clinical threshold" for PD [8,9,10].

Histological studies reveal that PD is characterized by the abnormal deposition of the insoluble alpha-synuclein protein, forming aggregates known as Lewy bodies [11]. These protein aggregates progressively accumulate throughout the brainstem and various neocortical and limbic regions, reflecting the progressive degeneration of the entire central nervous system (CNS) [12]. More recently, a neuroinflammatory state has been observed in the brains of PD patients, particularly evident in the SNc [13,14]. Both neuroinflammation and dysfunctional activation of the immune system within the CNS significantly contribute to PD pathology and pathophysiology [15]. The inflammasome, a crucial complex of immune-modulating receptors and sensors, plays a role in recruiting proteins associated with apoptotic mechanisms through caspase-1 activation [16,17]. Caspase-1, in turn, activates the pro-inflammatory cytokines interleukin (IL)-1β and IL-18, perpetuating the neuroinflammatory state in PD brains [17,19,20,21]. A recent proposed model suggests that alpha-synuclein activates the inflammasome in the SN, leading to IL-activated pro-inflammatory profiles, neuronal death, and clinical symptoms [22,23,24].

Genetically, several genes or gene variants, including leucine-rich repeat kinase 2 (LRRK2), synuclein alpha (SNCA), glucosylceramidase beta-1 (GBA1), Parkin RBR E3 ubiquitin protein ligase (PARKIN), and PTEN-induced kinase 1 (PINK1), have been implicated in causing PD [25,26,27,28]. Molecular profiling studies of post-mortem SNpc samples, aimed at identifying differential molecular expression changes specific to PD compared to controls, have been conducted [29]. For instance, Simunovic et al. [30] used RNA microarrays to analyze SNpc gene expression in PD samples, identifying dysregulation of known molecular regulatory pathways in PD, including dysfunction in mitochondrial and oxidative stress-induced cellular responses [31,32]. In a recent comparative gene expression analysis on laser-dissected neurons from SNpc, Zaccaria et al. [33] revealed 52 dysregulated genes in PD samples compared to controls.

In our study, we conducted mRNA analysis and subsequent enrichment using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) to assess mRNAs extracted postmortem from the SN of subjects with PD and healthy controls.

2. Results

In the examination of mRNA deregulation in the SNc of PD patients, we conducted gene expression profiling on 8 post-mortem SNpc samples from PD patients and 6 from healthy control subjects (CTRL) using next-generation sequencing (RNA-Seq). The aim was to identify specific and differential changes in molecular expression. Following the removal of low-quality reads and adapter sequences, the high-quality reads were aligned against the human genome reference (hg38). In detail, we identified a total of 16,148 transcripts (Supplementary Table 1, sheet A), with 92 mRNAs differentially expressed (DEGs) between the two groups (PD vs. CTRL). Among these, 33 mRNAs were significantly upregulated (Table 1; Figure 1), while 59 mRNAs were significantly downregulated in PD compared to CTRL (Table 2; Figure 1). The normalized count of mRNAs is available at ArrayExpress (E-MTAB-13295).

The heatmap (Figure 1A) illustrates statistically significant differences in mRNA expression profiles between PD and CTRL. The volcano plot (Figure 1B) depicts the distribution of differentially expressed transcripts by their fold change and p-values. The most upregulated genes are toward the right, the most downregulated toward the left, and the most statistically significant genes at the top.

We employed the pathfindR tool to analyze significant DEGs in PD based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. Enrichment analysis was performed to explore functional variations between the two groups and investigate pathways potentially associated with PD. Statistically significant signaling pathways included (hsa04260) cardiac muscle contraction, (hsa04727) GABAergic synapse, (hsa04140) Autophagy, and (hsa04666) Fc gamma R-mediated phagocytosis (Figure 2A and Supplementary Table 1).

For Gene Ontology (GO) enrichment, we identified both down-regulated and up-regulated genes involved in various Molecular Functions or Biological Processes (Figure 2B and Supplementary Table 2). Notable enrichments encompassed chemical synaptic transmission, cellular response to hypoxia, autophagosome assembly, GABA receptor binding, and cell communication by electrical coupling in both heart and brain activity.

To explore Diseases and Biological Functions significantly enriched in DEGs and assess potential associations with PD susceptibility in the SNc, we employed Ingenuity Pathway Analysis (IPA). The analysis revealed significant enrichment in Neurological Disease, with multiple annotations related to disorders of basal ganglia, movement disorders, neuromuscular disease, dyskinesia, progressive motor neuropathy, familial neurological disorder, Parkinson's disease, progressive neurological disorder, abnormal morphology of the nervous system, and tauopathy (Figure 3 and Supplementary Table 3). The DEGs IPA Network Analysis identified 7 networks with nodes and interactions associated with the top Diseases or Functions Annotations, such as Neurological Disease, Organismal Injury and Abnormalities, and Psychological Disorders (Figure 4 and Table 3).

Among the genes of interest, GABARAPL1 was found to be overexpressed, as indicated by both KEGG and GO analyses, involving a distinct set of pathways (5 and 3, respectively): GABAergic synapse, autophagy – animal, NOD-like receptor signaling pathway, autophagy – other, and mitophagy – animal (according to KEGG), as well as autophagosome, GABA receptor binding, and Tat protein binding (according to GO). Autophagy, a highly conserved cellular degradation process regulated by specific autophagy-related (Atg) factors, entails the formation of double-membrane autophagosomes that engulf cytoplasmic components for degradation. In mammals, this process is complex due to the presence of six Atg8 homologues, categorized into the GABA type A receptor-associated protein (GABARAP) and microtubule-associated protein 1 light chain 3 (MAP1LC3) subfamilies [39]. GABARAPL1/GEC1, a member of the GABARAP subfamily, exhibits the highest mRNA expression among Atg8 homologues in the CNS. Notably, GABARAPL1 brain expression is observable as early as embryonic day 11, increasing progressively to peak in adulthood. Significantly, GABARAPL1 expression in the adult brain is particularly intense in neurons involved in motor and neuroendocrine functions, notably in the SNc [40]. Dysregulation of Atg8 homologues has been observed in other synucleinopathies, such as Lewy-body dementia and multiple system atrophy [41]. In PD, alterations in autophagic mechanisms are evident, as demonstrated by transcript levels of several autophagy genes in blood cells. A recent study found overexpression of autophagy-related genes, including MAP1LC3B, GABARAP, GABARAPL1, GABARAPL2, and sequestosome 1 (P62/SQSTM1) in PD patients, with potential implications for predicting markers and therapeutic responses [42].

Similarly, the SLC38A2 gene was overexpressed and associated with two relevant pathways according to KEGG: GABAergic synapse and protein digestion and absorption. This finding strengthens the connection between GABA overexpression, PD pathology, and neural degeneration. In a recent study on a rotenone-induced PD rat model [43], nardosinone, a biochemical compound enhancing NGF-mediated neurite outgrowth and synaptogenesis, demonstrated anti-PD efficacy. Transcriptome and proteome analysis suggested that the anti-PD target of nardosinone is the SLC38A2 gene, potentially involving the GABAergic synaptic pathway. This underscores the SLC38A2 gene as a potential target for PD treatment and the modulatory effects of this gene as anti-PD agents through the GABA system [43].

Diseases and Function analysis in the IPA software highlighted that DEGs correlated with Neurological Disease (Figure 3, Figure 4A, Figure 4B). Network IPA Analysis (Figure 4B) revealed downregulation of many mitochondrial genes in the SN of PD subjects, including MT-CO1, MT-CO2, MT-CO3, MT-CYB, and Cytochrome-oxidase. This further supports mitochondrial dysfunction in PD, where mitochondria, involved in crucial functions, primarily energy generation, are essential for nearly all cellular activities. Alterations in mitochondrial functioning lead to insufficient energy production, particularly affecting the CNS [45]. Evidence indicates that mitochondrial respiratory chain dysfunction plays a primary role in various neurodegenerative diseases, including PD [45,46]. Impaired elements of the respiratory chain, such as defects in complex I, have been associated with PD and frontal cortex dysfunction [47,48,49,50]. Damage to the electron transport chain increases oxidative stress and neuronal dysfunction, potentially contributing to the onset and progression of PD. Progressive mitochondrial damage results in the accumulation of non-functional mitochondria, further contributing to neuronal degeneration [53].

Looking at the downregulated genes in the SN of PD patients, a group of genes involved in maintaining the structure and function of glial cells was noted. These genes include MBP (myelin basic protein), CPN (myelin protein cyclic nucleotide phosphodiesterase), CTNNA3 (alpha-T-catenin), AQP1 (aquaporin 1), and GFAP (glial fibrillary acidic protein). Among the upregulated genes, SLC8A1 (solute carrier family 8 member A1) is involved in linking trans-plasmalemmal gradients of sodium and calcium ions to the membrane potential of astrocytes. This outcome suggests the emerging important role of glial cells in neurodegeneration and PD pathogenesis [59].

A limitation of this study is the small number of PD and control subjects recruited; however, the ability to sample human brains and extract mRNA sequences directly from the SN of both patients and controls is a significant strength of this report.

3. Materials and Methods

3.1. Human post-mortem midbrain samples

Human midbrain samples were generously provided by the Parkinson's UK Brain Bank (Imperial College London, London, UK). A total of 8 PD and 6 CTRL samples were acquired, and the specimens were histologically sectioned at the midbrain, encompassing the human SN in all slides. Each section had a thickness of 4 µm. Supplementary Table 4 provides a summary of the clinical features of the included subjects. The study adhered to the principles of the Declaration of Helsinki of 1964 and its subsequent amendments. The Ethics Committee of the Oasi Research Institute—IRCCS of Troina (Italy) granted approval for the protocol on April 5th, 2022 (approval code: 2022/04/05/CE-IRCCS-OASI/52).

3.2. RNA isolation from human midbrain samples

RNA was extracted from 4 µm FFPE slide-mounted sections using the RecoverAll Total Nucleic Acid Isolation Protocol (ThermoFisher Scientific), following the manufacturer's instructions. Subsequently, the RNA was stored at -80°C until further processing.

3.3. RNA Sequencing and Functional Analysis

RNA sequencing and subsequent data analysis were conducted by Genomix4Life Srl (Baronissi, Italy). The quality and quantity of RNA were assessed using a Qubit fluorometer (Thermo Fisher Scientific) and a TapeStation 4200 (Agilent Technologies, 5301 Stevens Creek Blvd, Santa Clara, CA, USA), respectively.

Indexed libraries were prepared from 50 ng of purified RNA each, employing the Illumina Stranded Total RNA with Ribo-Zero Plus Kit (Illumina), following the manufacturer's guidelines. Library quantification was performed using the TapeStation 4200 (Agilent Technologies) and Qubit (Thermo Fisher). Subsequently, the indexed libraries were pooled in equimolar amounts, resulting in a final concentration of 2 nM. The Illumina NovaSeq 6000 System was utilized to sequence the pooled samples in a 2 × 75 paired-end format.

Quality control analysis was executed on the raw fastq files generated earlier using the FastQC tool [http://www.bioinformatics.babraham.ac.uk/projects/fastqc/]. Trim-ming of short reads (<25 bp) and removal of adapter sequences were carried out with cutadapt (v.2.8) [60]. The trimmed fastq files were mapped to the reference genome (GenCode (HG38-Release 37 (GRCh38.p13)) [https://www.gencodegenes.org/human/]) using STAR (version 2.7.3a) [61] with standard parameters. Gene quantification per sample was accomplished using featureCount (version 2.0) [https://bioweb.pasteur.fr/docs/modules/subread/1.4.6-p3/SubreadUsersGuide], and a custom R script was employed for data normalization and differential expression analysis using the Bioconductor DESeq2 [62] package. The threshold for considering genes as differentially expressed was set at Fold-Change ≥ 1.50 or ≤ -1.50 (|FC| ≥ 1.50) with adjusted p-values ≤ 0.05 (padj). Volcano plots and heatmaps were generated using the EnhancedVolcano (10.18129/B9.bioc.EnhancedVolcano) and ComplexHeatmap [63] packages in R.

For functional analysis, KEGG pathway and GO database analyses were conducted using the R package pathfinder [https://doi.org/10.3389/fgene.2019.00858]. Additionally, to gain a deeper understanding of the complex transcriptomics data, Ingenuity Pathway Analysis (IPA) [64] was performed, specifically for investigating diseases and function analysis.

The raw data (.fastq files) and the normalized count of identified mRNAs are available on ArrayExpress (E-MTAB-13295).

4. Conclusions

Despite the limited sample size, the comprehensive analysis conducted on the transcriptome of these samples, along with KEGG and GO enrichment analyses, revealed statistically significant findings associated with the following signaling pathways: cardiac muscle contraction, GABAergic synapse, autophagy, Fc gamma R-mediated phagocytosis signaling pathway, cellular response to chemical synaptic transmission to hypoxia, autophagosome assembly, GABA receptor binding, and cellular communication via electrical coupling involved in cardiac conduction.

Upon examining the genes that were either upregulated or downregulated, these results contribute to a translational perspective and unveil novel pathophysiological insights into PD. This includes implications in cell communication via electrical coupling, specific mitochondrial activity, and neuronal-glia interaction. These findings may pave the way for enhanced diagnostic precision and the development of novel targeted therapies.