RNA Sequencing of Parkinson’s Ventral Midbrain to Estimate Gene Expression Reveals Marked Heterogeneity and GADD45A as a Potential Therapeutic Target

James P. Bennett, Jr.1* and Paula M. Keeney1

1Neurodegeneration Therapeutics, Inc
3050A Berkmar Drive
Charlottesville, VA 22901-3450
PH: 434-529-6457
FAX: 434-529-6458
aging.mitochondria@gmail.com
www.NDTherapeutics.org

Abstract

Parkinson’s disease (PD) is the second most prevalent neurodegenerative disease of humans, is a progressive disorder that can begin outside of the central nervous system and involves multiple brain areas. We used high-density Illumina® RNA sequencing (RNA-seq) to examine the transcriptomes of ventral midbrain samples from 14 PD and 9 CTL age-, postmortem interval- and RNA quality-matched postmortem brain samples. We aligned de-multiplexed sequences against the latest UCSC-hg38 human genome using tophat2/bowtie2, then processed the resulting .bam files with either Cufflinks (FPKM normalization) or Qlucore Omics Explorer® (QOE, FPKM or TMM normalization). Gene expression data were processed using FDR (false discovery rates, Cufflinks and QOE) to compare multiple samples, PCA (principal component analysis-QOE) and GSEA (gene set expression analysis-QOE).

PD samples showed >65% loss of dopaminergic markers. At FDR<5%, 153/177 genes (86.4%) were under-expressed in PD by Cufflinks and 199/351 (56.7%) by QOE (FPKM normalization). 26 genes were co-identified by both Cufflinks and QOE (FPKM normalization) at FDR<5%, and their PD/CTL ratios (25/26 < 1.0) were linearly related (r2=0.992). GSEA analysis of PD samples showed greatest enrichment in cholesterol homeostasis (31/73 genes), mTORC1 signaling (87/198 genes) and oxidative phosphorylation (123/199 genes). PD samples were markedly heterogenous as reflected in plots of gene expression variance in PD vs CTL samples, and in PD samples against gene expression ratios (PD/CTL). The stress-response gene GADD45A/α showed small variance and was under-expressed in PD. Our systems biology approaches suggest that increased expression of GADD45A/α could benefit the majority of PD subjects.

Introduction

Parkinson’s disease (PD) is a common neurodegenerative disease of adults and ranks second behind Alzheimer’s disease as an age-related affliction of societies. 7-10 million persons worldwide have PD, which is commonly (>90%) sporadic and afflicts men about 1.5 times as often as women. Older age is the major known risk factor for PD. Rare mutations of chromosomal genes are believed to be causal for autosomal dominantly or recessively inherited PD that can have variable postmortem neuropathology. [1-5]
Motor symptoms of PD are currently the basis of clinical diagnosis and are characterized by progressive bradykinesia of speech and limbs, rigidity, postural instability that can lead to increased falls and, at least once during the illness, a 4-6 Hz resting tremor of face/tongue and/or limbs. These motor deficits of tremor, rigidity and bradykinesia are believed to arise from progressive loss of dopaminergic (DA) neurons in the substantia nigra pars compacta. PD is commonly complicated by “secondary” symptoms, which include deficits of neuropsychiatric function, autonomic disturbances, and sleep disorders. [6-9]
Frequently, affected nigral neurons possess eosinophilic inclusions containing aggregated and post-translationally modified alpha-synuclein (A-SYN), among other proteins [5, 10-15]. Nigral DA neurons in PD show ultrastructural changes of both apoptosis and autophagy [16].

Neuropathological examinations of many hundreds of autopsies carried out by the Braak group have contributed significantly to our collective understanding of the correlations of PD progression and neuropathology. These studies revealed that aggregated A-SYN appears initially in gastrointestinal enteric neurons, followed temporally by rostral spread into dorsal medulla, midbrain, limbic system and frontal cortex. Braak, et al have divided PD into 6 stages, with the first 3 stages representing pre-clinical (in terms of motor deficits) progression, and the latter 3 stages representing progression of motor deficits and appearance of important secondary symptoms such as mood disorders and other neuropsychiatric symptoms and cognitive impairment. These latter 3 stages of PD correlate with spread of A-SYN aggregation pathology into rostral limbic and frontal lobe regions. [7-12, 17-40]

PD motor deficits are commonly treated by dopaminergic agents, either L-DOPA, the immediate precursor of lost dopamine (DA), or DA-ergic agonists, sometimes accompanied with agents that interfere with L-DOPA methylation or DA catabolism. These approaches are typically effective for 3-5 years, but commonly are followed by erratic responses to DA-ergic agents and troublesome side effects. Clinical responses can be improved, sometimes with reduced and/or more frequent DA-ergic drug dosing, or by continuous deep brain stimulation of basal ganglia nuclei, such as globus pallidus interna (GPi) or subthalamic nucleus (STN) [2, 3, 5, 41-44].

There are presently no known approaches that retard progression of underlying PD pathology. Should such approaches become available, then management of early clinical motor stages of PD, which generally are straightforward and successful, could be temporally extended, and later disabilities frequently derived from secondary symptoms [45] might be avoided.

Analyses of PD brain transcriptomes have yielded complex results. [46] Transcriptomes in postmortem brain represent attempts of cells to respond to stresses present at the time of death. These stresses may or may not reflect processes causal to neurodegeneration. In addition, individuals may vary in their molecular responses to similar stresses, further complicating the picture.

Because one goal of transcriptomic studies of PD postmortem brain is to attempt defining therapies that may alter deficits or overactivities of pathways that are defective in the majority of PD subjects, not a selective minority (e.g., those with a rare genetic cause), we undertook this study to try and find such a therapeutic approach. We used RNA sequencing of postmortem ventral midbrain, which consists mainly of glia and residual dopaminergic (and other phenotypic) neurons. We found significant heterogeneity of gene expression but were able to define underexpression in the majority of PD samples of a stress-activated gene (GADD45A/α) that alters many intracellular events. We propose that augmentation of expression of that gene and similar genes offers the potential to restore protective pathways in PD neurodegeneration.

Methods

Tissue Samples

Blocks of snap frozen human substantia nigra/ventral mid brain were obtained from the University of Virginia Brain Resource Facility, the Virginia Commonwealth University Brain Tissue Resource Facility or the National Disease Research Interchange (NDRI), Philadelphia, PA. Demographics of the samples are presented in Supplemental Table 1.

Isolation of total RNA

Briefly, the area of interest was dissected from fresh post mortem midbrain, immediately frozen in liquid nitrogen-chilled 2-methly butane and stored at -80 oC until sectioned. At the time of sectioning, tissue blocks were trimmed to about 1 by 1.5 cm, embedded in Neg-50 embedding compound (Richard-Allan Scientific, Kalamazoo, MI) and sectioned at -16oC. 25 sections (20 um) were collected from each sample and placed in Qiazol (Qiagen, Germantown, MD). Total RNA was isolated according to miRNeasy instructions with sonication of the tissue in Qiazol, on-column DNase digestion and an additional RPE wash. RNA quantity and quality was assessed using an Agilent Tape Station. Total RNA was sent to CoFactor Genomics for preparation and sequencing, see below.

RNA-seq and analyses

Total RNA was rRNA depleted, multiplex sequencing libraries were constructed and quantitated, and Illumina RNA paired-end 150 base sequencing carried out at CoFactor Genomics as described in [47, 48]. All bioinformatics analyses were performed “blind”, with the samples identified only as numbers and not by disease. Sequencing fastq files were examined for quality in FastQC, then Illumina sequencing adaptors removed using Trimmomatic®. Resulting sequencing files with minimum Phred scores of 33 were aligned against the UCSC-hg38 version of the human genome using Tophat2/Bowtie2®. The resulting *.bam files were processed either by Cufflinks® to yield gene expression normalized and expressed by the FPKM (Fragments Per Kilobase of exon sequenced per Million reads) method, or processed using Qlucore Genomics Explorer® with the NGS plug-in (Qlucore Genomics, Inc., www.qlucore.com) and normalized by either the TMM (Trimmed Mean of M-values) or FPKM methods. The same genome.fa and genes.gtf files were used by all programs.

Graphing and Statistics

All graphing and statistics were performed in GraphPad Prism version 7 (for Mac OSX) (www.graphpad.com). QOE (www.qlucore.com) was version 3.3(80) and included the NGS (Next Generation Sequencing) plug-in. All bioinformatic analyses were carried out on a Macintosh Pro desktop computer operating with OSX v 11 and possessing a 12-core processor and 64 GB of RAM. Sequencing data and output files were stored on a 32TB server.

Results

PD ventral midbrain samples show loss of expression of multiple dopaminergic markers
By the time an individual dies with advanced PD, s/he has lost the majority of dopaminergic neurons in the substantia nigra [2, 13, 14]. Dopaminergic neuron loss in the more dorsal ventral tegmental area is variable in PD [13]. Our PD ventral midbrain samples (Supplemental Table 1) showed loss of expression of tyrosine hydroxylase (TH; PD/CTL (ratios of mean PD FPKM/mean CTL FPKM) = 0.332, q (false discovery rate metric) = 0.058); dopamine transporter (DAT/SLC6A3; PD/CTL = 0.309, q=0.062); VMAT2/SLC18A2; (PD/CTL = 0.360; q= 0.063). Expression of dopa-decarboxylase, which can decarboxylate multiple L-aromatic amino acids in addition to L-DOPA, was also reduced (DDC; PD/CTL= 0.414; q=0.087).

Expression of multiple autosomal genes mutated in rare familial PD were either reduced or not statistically present in our ventral midbrain samples.

Expression of alpha-synuclein (ASYN/SNCA), the initial gene identified as mutated in familial PD [49, 50], was reduced (PD/CTL = 0.557; q=0.044), as was that of UCHL1 (PD/CTL = 0.62; q=0.056). Additional autosomal genes, mutations of which have also been associated with familial PD, were not found in our list of genes with FDR’s < 10% (q<0.1). These included LRRK2, PINK1, parkin/PARK2, DJ-1/PARK7, UPS35 and EIF4G1.

PCA plots demonstrate separable PD and CTL RNA-seq populations

PCA (principal component analysis) plots were constructed in QOE after processing the Tophat2/BowTie2-generated *.bam RNA-seq files using the NGS (Next Generation Sequencing) plug-in. Figure 1a shows the PCA plot using TMM normalization and Figure 1b shows the PCA plot obtained using FPKM normalization for the same *.bam files. Both plots showed that at a FDR (q) of <5% it was possible to separate the two populations (PD, CTL). Networks could also be constructed in QOE using TMM (Figure 1a) or FPKM (Figure 1b) normalization. Note that the same PD sample (160) overlapped in network analysis with the CTL group in both TMM and FPKM PCA plots.

 

Figure 1. Principal component analysis (PCA) plots generated by QOE with q<5% and fold-change > 1.5. 1a. shows PCA plot generated with TMM normalization and 1b shows PCA plot generated with FPKM normalization. Note that in both 1a and 1b PD sample 160 generates an edge with one CTL sample.

Cufflinks and QOE processing showed small overlap among DEG’s (differentially expressed genes) but with nearly identical expression ratios

Processing of the *.bam sequencing files against the same hg38 genes.gtf file yielded mean FPKM values (Cufflinks) or “fold-change” values (QOE). Controlling for multiple comparisons (q<5%) yielded 177 differentially expressed genes (DEG’s) using Cufflinks and 351 DEG’s (by QOE) (Supplemental Table 2). With Cufflinks 153/177 genes (86.4%) had PD/CTL ratios < 1.0, and with QOE using FPKM normalization 199/351 genes (56.7%) had PD/CTL ratios < 1.0. There were 26 genes in common in these two datasets (Supplemental Table 3), and Figure 2 shows that the PD/CTL fold-changes in expression values derived from each analytical approach were highly correlated with each other (r2=0.992). 25/26 genes in this group were under-expressed in the PD datasets.

Figure 2. Linear relationship among 26 genes co-identified with q < 5% using TMM normalization (QOE, y-axis ) compared to q < 5% in Cufflinks (x-axis). These genes are listed in Supplemental Table 3.

Heatmaps also demonstrated separability of the PD and CTL populations

Figure 3a (top, left) shows a QOE-generated heatmap of genes hierarchically clustered with q <5% and FPKM normalization, and Figure 3b shows a similar plot with q<5% and TMM normalization. In both of these heatmaps, hierarchical clustering “mis-labeled” one PD sample as belonging to the CTL cluster group. This sample (160) was the same PD sample that interacted with CTL samples in network analysis of the PCA plots (Figures 1a, 1b) Lowering of the q metric resulted in proper heirarchical cluster separation of PD and CTL samples (Figures 3c, 3d).

Figure 3. Heatmaps generated in QOE using q < 5%, hierarchical clustering and either FPKM (3a) or TMM (3b) normalization. Note that PD160 (black arrow) was mis-classified in both clusterings. When q is reduced to < 4.8% (3c, FPKM normalization) or < 2.9% (3d, TMM normalization), PD160 is properly clustered in the PD group.

GSEA analyses show gene enrichment in cholesterol homeostasis, mTORC1 signaling and oxidative phosphorylation

Figure 4 shows heatmaps derived from gene set enrichment analysis (GSEA) analysis (using QOE) of cholesterol homeostasis, mTORC1 signaling and oxidative phosphorylation, the three gene sets having the highest relative enrichment in the PD samples. We used the most current “hallmark” GSEA file for these calculations. GSEA delineated that genes were represented in cholesterol homeostasis (31/73 genes), mTORC1 signaling (87/198 genes) and oxidative phosphorylation (123/199 genes).

Figure 4. Heatmaps from GSEA analyses (using the “hallmark” group h.all.v6.1.symbols.gmt, http://software.broadinstitute.org/gsea/downloads.jsp) of differentially expressed genes (DEG) using q < 5% and QOE with TMM normalization. Shown are heatmaps from the three groups with the highest “normalized enrichment”: oxidative phosphorylation (left), cholesterol homeostasis (center) and mTORC1 signaling (right). The vertical black lines define the genes which met enrichment criteria. Note that in the PD samples most genes that met enrichment criteria are under-expressed.

Gene expression in PD brain samples exhibited heterogeneity.

Figure 5 shows plots of variance of gene expression (variance = standard deviation squared) in PD samples (from Cufflinks, y-axis) compared on a gene-by-gene basis to variance of expression in CTL samples (from Cufflinks, x-axis). The lower plot is an amplification of axes used in the upper plot. As can be seen, there was significant scatter in the variance values, with no clear relationships between PD and CTL.

Figure 5. Plots of gene-by-gene variance (std dev squared) in the PD population (y-axis) compared to variance in the CTL population (x-axis). Data are from Cufflinks (FPKM) analysis. The bottom plot is from the same dataset as the top plot except the axes are enlarged.

Figure 6 shows a gene-by-gene plot of gene expression variance in the PD samples (from Cufflinks, y-axis) versus mean PD/mean CTL gene expression ratio (from Cufflinks, x-axis). The lower plot is an amplification of the axes shown in the upper plot. There is significant scatter in the values, with a suggestion of a Gaussian-type distribution. The vertical blue lines represent PD/CTL expression ratios of 1.0 (equal expression in PD and CTL).

The orange rectangle in the lower plot encompasses the genes with lowest PD variance (< 50) and expression ratios (PD/CTL) <0.5 (under-expressed in PD > 2.0-fold). These 33 genes (Supplemental Table 4), if over-expressed, should theoretically represent therapeutic targets for the majority of PD subjects.

Figure 6. Plots of gene-by-gene variance in the PD population (y-axis) vs expression ratio (mean PD FPKM/mean CTL FPKM) (x-axis). Data are from Cufflinks (FPKM) analysis. The bottom plot is the same as the top plot except the axes are enlarged. The orange rectangle represents genes whose variance is <50 and whose expression ratio is < 0.5 (>2.0-fold under-expressed in PD). These 33 genes are listed in Supplemental Table 4.

The stress-response gene GADD45α stands out in this group for several reasons. 1.It has among the lowest variances in the PD group (5.60) and is under-expressed in PD (PD/CTL = 0.39) ; 2. It possesses pleotropic actions on multiple cell processes. 3.It may be involved in PD neurodegeneration (see Discussion for more details)

Discussion

In this study we used RNA sequencing (RNA-seq) to characterize gene expression in PD and CTL ventral midbrain samples that were matched for age, postmortem intervals and RNA purity. RNA-seq, although still expensive, is rapidly becoming the state-of-the-art for gene expression studies, having replaced microarray platforms and quantitative PCR (qPCR), both 20th century technologies. There are almost 200,000 PubMed citations for this technique, and publicly available software for bioinformatic analysis of RNA-seq data has been available for ~ a decade. Commercial software is also becoming available, and we compared publicly available bioinformatic software (Tophat2/Bowtie2/Cufflinks) to that available commercially (QOE). Note that the commercial software we used still required generation of binary sequencing files (.bam files), which we did using Tophat2/Bowtie2 and the most current version (hg38) of the human genome.

Although RNA-seq is powerful and generates large datasets (which can be challenging to analyze), when applied to postmortem human brain samples from individuals that shared a common clinical phenotype, RNA-seq data must be interpreted with caution. There are multiple reasons for this caution, three of which are discussed below:

First, most individuals whose brains were donated for storage and study likely died at end-stage of disease. In the case of neurodegenerative diseases, RNA-seq may reveal molecular strategies of “survivors”, not the vulnerable neurons that have already died. Observed gene expression likely represents events mostly occurring in astroglia and represents little contribution from vulnerable neuronal populations. “Genesis events” for neurodegeneration may also not be revealed. This limitation may be addressed by use of cell lines derived from living donors, but even this approach must be used cautiously [51].

Second, although humans share much genetic material in common, each person is genetically/epigenetically unique, and there is no a priori reason to assume that common molecular pathogenic abnormalities will emerge from RNA-seq studies. In fact, our findings indicate great degrees of heterogeneity of gene expression across tissue samples that included quantitation of expression of thousands of genes.

Third, alterations in RNA expression may represent causal changes or may represent adaptive responses to causal events. One cannot assume that changes in RNA expression belong to either group.
These arguments suggest that our collective attempts in the past to define common sets of altered genes, even after corrections for multiple samples, is problematic in terms of defining both biological origins/changes and therapeutic relevance. In fact, there is little agreement across multiple gene expression studies as to alterations of specific genes (see [46] for analysis of recent Parkinson’s transcriptomic studies).

Our findings demonstrate the limited overlap among analyses that used the same *.bam sequencing files to estimate expression of genes. Even with use of common FPKM normalization, one must be cautious, as we observed very limited overlap of individual genes identified by what should be very similar approaches.
We do not believe there is a “right” or “wrong” approach. Rather, one must be aware of potential differences in estimation of gene expression. For comparative RNA-seq studies, one should use raw sequencing files and process them in the same manner.

GADD45A/α

The stress-response gene GADD45A/α (Growth Arrest and DNA Damage-inducible 45 A/α) is a rapidly-induced member of the GADD45 family that itself is non-enzymatic but serves multiple critical cell survival functions through nucleic acid and protein partner binding. Among these is the response to genotoxic DNA damage mediated by p53 induction of GADD45a expression.

GADD45α also serves an important epigenetic role by physically interacting with ten-eleven-translocation (TET) proteins such as TET1 [52]. These dioxygenase proteins convert methyl cytosine to hydroxymethyl cytosine, which generally results in increased gene expression. GADD45α can also demethylate DNA in the course of p53 stimulation of DNA repair [53, 54] providing another theoretical mechanism to increase gene expression. GADD45 members are also involved in apoptosis, inflammation, resistance to oxidative stress and cellular senescence. [55-57]

GADD45α can be induced by vitamin D3 (1,25(OH)2D3) [58], by exposure to the DA-ergic neurotoxin MPP+ [57], or by fucoxanthin, a marine carotenoid found in brown algae [59], which is available as a dietary supplement. We speculate that increased expression of GADD45A/α will counter neurodegenerative mechanisms in the majority of PD subjects and slow disease progression.

mTORC1

Mechanistic (formerly mammalian) target of rapamycin (mTOR) is a protein kinase member of the PI3K family that combines with other protein subunits to form mTOR complexes 1 (mTORC1) and 2 (mTORC2). mTORC1 is sensitive to inhibition by rapamycin and analogues (“rapalogues”), whereas mTORC2 is relatively insensitive to inhibition by rapamycin and rapalogues. Both complexes are believed to act by sensing external nutrients and to control cell growth and differentiation by modifying downstream processes such as protein translation, cell division and growth, and autophagy. [60-64]

Our finding of decreased mTORC1 signaling in the majority of PD samples following GSEA analyses suggests that this crucial regulatory system is impaired in PD (at least in end stage disease). Reduced mTORC1 signaling would be predicted to increase autophagy, which would likely take place since we observed that all autophagy genes (ATGx) and beclin-1 (BECN1) were normally expressed in PD samples (PD/CTL ratios ~0.9-1.1). Activating autophagy has been proposed as a treatment for PD [65], although consideration must be given to the alternative viewpoint, that autophagy is contributing to PD pathology [16, 66]

The combination of impaired mitobiogenesis signaling [67, 68] and increased oxidative stress (which would further damage mitochondria) [46, 69, 70], sets the stage for reduced mitochondrial mass (impaired mitobiogenesis signaling, increased oxidative damage to mitochondrial components). Thus, our findings of decreased oxidative phosphorylation in GSEA analyses of the PD samples is not surprising, although one can postulate multiple origins of this energy deficit.

The above speculations do not identify a single “genesis event” for PD neurodegeneration, nor can we account for all abnormalities observed. Rather multiple interacting deficits could appear slowly over time, interacting to amplify each other’s impacts on cell survival. That concept is a major thrust of systems biology, to identify a series of interacting processes rather than individual molecules that, when combined, interact pathogenically.

In addition, it must be recalled that our samples are devoid of the majority of PD vulnerable neurons (nigral dopaminergic neurons). These neurons likely have increased oxygen consumption compared to other neurons due to mitochondrial respiration, and activities of MAO and tyrosine hydroxylase, which in turn have higher oxygen consumption compared to non-neuronal cells that likely comprise the majority of our samples. It appears likely that increased oxidative stress could be an early event in the death spiral of dopaminergic neurons in PD, but that by itself increased oxidative stress cannot likely account completely for the pathophysiology of neurodegeneration.

Summary

We propose that gene variance, or some other measure of gene variation across samples, be utilized in transcriptomic studies of human diseases, so as to define genes with both minimal variability in expression across individuals and significantly reduced expression in the disease population. Manipulation (specifically, increased expression) of these genes would be predicted to have the greatest therapeutic benefit across multiple persons, and thus they are appealing therapeutic targets. Alteration of expression of these “low variance” genes could also serve as biomarkers to predict therapeutic utility of gene expression approaches. Such approaches could seek different pathways for altering expression, including direct expression using gene vectors, administration of drugs that increase expression, or alteration of non-coding RNA’s (ie microRNA’s, etc) that indirectly alter gene expression.

Using this approach, we identified the important stress-response gene GADD45A/α as a viable candidate for therapeutic use in PD. We speculate that increased brain expression of GADD45A will slow disease progression in most PD patients.

In addition, we further speculate that sporadically occurring PD neurodegeneration is both highly heterogeneous in etiology and unlikely to be successfully managed by a traditional single protein targeting approach. Rather, we propose that low variance/low expression genes be identified in PD populations and tested for their efficacy in altering disease progression.

Figure Legends

Figure 1. Principal component analysis (PCA) plots generated by QOE with q<5% and fold-change > 1.5. 1a. shows PCA plot generated with TMM normalization and 1b shows PCA plot generated with FPKM normalization. Note that in both 1a and 1b PD sample 160 generates an edge with one CTL sample.

Figure 2. Linear relationship among 26 genes co-identified with q < 5% using TMM normalization (QOE, y-axis ) compared to q < 5% in Cufflinks (x-axis). These genes are listed in Supplemental Table 3.

Figure 3. Heatmaps generated in QOE using q < 5%, hierarchical clustering and either FPKM (3a) or TMM (3b) normalization. Note that PD160 (black arrow) was mis-classified in both clusterings. When q is reduced to < 4.8% (3c, FPKM normalization) or < 2.9% (3d, TMM normalization), PD160 is properly clustered in the PD group.

Figure 4. Heatmaps from GSEA analyses (using the “hallmark” group h.all.v6.1.symbols.gmt, http://software.broadinstitute.org/gsea/downloads.jsp) of differentially expressed genes (DEG) using q < 5% and QOE with TMM normalization. Shown are heatmaps from the three groups with the highest “normalized enrichment”: oxidative phosphorylation (left), cholesterol homeostasis (center) and mTORC1 signaling (right). The vertical black lines define the genes which met enrichment criteria. Note that in the PD samples most genes that met enrichment criteria are under-expressed.

Figure 5. Plots of gene-by-gene variance (std dev squared) in the PD population (y-axis) compared to variance in the CTL population (x-axis). Data are from Cufflinks (FPKM) analysis. The bottom plot is from the same dataset as the top plot except the axes are enlarged.

Figure 6. Plots of gene-by-gene variance in the PD population (y-axis) vs expression ratio (mean PD FPKM/mean CTL FPKM) (x-axis). Data are from Cufflinks (FPKM) analysis. The bottom plot is the same as the top plot except the axes are enlarged. The orange rectangle represents genes whose variance is <50 and whose expression ratio is < 0.5 (>2.0-fold under-expressed in PD). These 33 genes are listed in Supplemental Table 4.

Supplemental Tables:

TABLE 1. Demographics of PD and CTL ventral midbrain samples used in this study.

TABLE 2. Genes identified at q<5% by Cufflinks or QOE (using FPKM normalization)

TABLE 3. Common genes identified by both Cufflinks and QOE (FPKM normalization) using q<5%

TABLE 4. 33 genes identified by Cufflinks with variance (Std.Dev.2) <50 and PD/CTL expression ratios < 0.5

Acknowledgements. This research was supported by internal funds from Neurodegeneration Therapeutics, Inc., an IRS-registered 501(c)3 non-profit medical research company. We are grateful to the many patients and families for their generous donations of postmortem brains. We particularly acknowledge the help provided by Dr. S. Churn of the VCU Brain Tissue Resource Facility that is supported by institutional funds. We also acknowledge Dr. Sara Strandberg for providing assistance in use of Qlucore software.

Author Contributions. JPB designed the study, performed all bioinformatics analyses and wrote the paper. PMK supervised the acquisition of all brain samples, sectioned all samples, extracted RNA and supervised their analyses. Both authors have reviewed and approved the final manuscript, and neither author declares any conflicts of interest.

Study data. Compressed (.gz) raw sequencing files in fastq format are the property of Neurodegeneration Therapeutics, Inc and are available from the corresponding author following execution of a Material Transfer Agreement and provision of either a FTP site URL or memory storage device sufficient for ~200GB of data.

REFERENCES CITED

[1] Delamarre A, Meissner WG (2017) Epidemiology, environmental risk factors and genetics of Parkinson’s disease. Presse Med 46, 175-181.
[2] Erkkinen MG, Kim MO, Geschwind MD (2017) Clinical Neurology and Epidemiology of the Major Neurodegenerative Diseases. Cold Spring Harb Perspect Biol.
[3] Srivanitchapoom P, Pitakpatapee Y, Suengtaworn A (2018) Parkinsonian syndromes: A review. Neurol India 66, S15-S25.
[4] Tysnes OB, Storstein A (2017) Epidemiology of Parkinson’s disease. J Neural Transm (Vienna) 124, 901-905.
[5] Radhakrishnan DM, Goyal V (2018) Parkinson’s disease: A review. Neurol India 66, S26-S35.
[6] Albers JA, Chand P, Anch AM (2017) Multifactorial sleep disturbance in Parkinson’s disease. Sleep Med 35, 41-48.
[7] Braak H, Rub U, Del Tredici K (2006) Cognitive decline correlates with neuropathological stage in Parkinson’s disease. J Neurol Sci 248, 255-258.
[8] Del Tredici K, Braak H (2013) Dysfunction of the locus coeruleus-norepinephrine system and related circuitry in Parkinson’s disease-related dementia. J Neurol Neurosurg Psychiatry 84, 774-783.
[9] Wolters E, Braak H (2006) Parkinson’s disease: premotor clinico-pathological correlations. J Neural Transm Suppl, 309-319.
[10] Braak E, Sandmann-Keil D, Rub U, Gai WP, de Vos RA, Steur EN, Arai K, Braak H (2001) alpha-synuclein immunopositive Parkinson’s disease-related inclusion bodies in lower brain stem nuclei. Acta Neuropathol 101, 195-201.
[11] Braak H, Braak E, Yilmazer D, Schultz C, de Vos RA, Jansen EN (1995) Nigral and extranigral pathology in Parkinson’s disease. J Neural Transm Suppl 46, 15-31.
[12] Del Tredici K, Braak H (2016) Review: Sporadic Parkinson’s disease: development and distribution of alpha-synuclein pathology. Neuropathol Appl Neurobiol 42, 33-50.
[13] Dickson DW (2018) Neuropathology of Parkinson disease. Parkinsonism Relat Disord 46 Suppl 1, S30-S33.
[14] Dickson DW, Braak H, Duda JE, Duyckaerts C, Gasser T, Halliday GM, Hardy J, Leverenz JB, Del Tredici K, Wszolek ZK, Litvan I (2009) Neuropathological assessment of Parkinson’s disease: refining the diagnostic criteria. Lancet Neurol 8, 1150-1157.
[15] Savica R, Bradley BF, Mielke MM (2018) When Do alpha-Synucleinopathies Start? An Epidemiological Timeline: A Review. JAMA Neurol.
[16] Anglade P, Vyas S, Javoy-Agid F, Herrero MT, Michel PP, Marquez J, Mouatt-Prigent A, Ruberg M, Hirsch EC, Agid Y (1997) Apoptosis and autophagy in nigral neurons of patients with Parkinson’s disease. Histol Histopathol 12, 25-31.
[17] Braak H, Bohl JR, Muller CM, Rub U, de Vos RA, Del Tredici K (2006) Stanley Fahn Lecture 2005: The staging procedure for the inclusion body pathology associated with sporadic Parkinson’s disease reconsidered. Mov Disord 21, 2042-2051.
[18] Braak H, Braak E (2000) Pathoanatomy of Parkinson’s disease. J Neurol 247 Suppl 2, II3-10.
[19] Braak H, Braak E, Yilmazer D, de Vos RA, Jansen EN, Bohl J, Jellinger K (1994) Amygdala pathology in Parkinson’s disease. Acta Neuropathol 88, 493-500.
[20] Braak H, de Vos RA, Bohl J, Del Tredici K (2006) Gastric alpha-synuclein immunoreactive inclusions in Meissner’s and Auerbach’s plexuses in cases staged for Parkinson’s disease-related brain pathology. Neurosci Lett 396, 67-72.
[21] Braak H, Del Tredici K (2017) Neuropathological Staging of Brain Pathology in Sporadic Parkinson’s disease: Separating the Wheat from the Chaff. J Parkinsons Dis 7, S73-S87.
[22] Braak H, Del Tredici K (2009) Neuroanatomy and pathology of sporadic Parkinson’s disease. Adv Anat Embryol Cell Biol 201, 1-119.
[23] Braak H, Del Tredici K (2008) Cortico-basal ganglia-cortical circuitry in Parkinson’s disease reconsidered. Exp Neurol 212, 226-229.
[24] Braak H, Del Tredici K, Bratzke H, Hamm-Clement J, Sandmann-Keil D, Rub U (2002) Staging of the intracerebral inclusion body pathology associated with idiopathic Parkinson’s disease (preclinical and clinical stages). J Neurol 249 Suppl 3, III/1-5.
[25] Braak H, Del Tredici K, Rub U, de Vos RA, Jansen Steur EN, Braak E (2003) Staging of brain pathology related to sporadic Parkinson’s disease. Neurobiol Aging 24, 197-211.
[26] Braak H, Ghebremedhin E, Rub U, Bratzke H, Del Tredici K (2004) Stages in the development of Parkinson’s disease-related pathology. Cell Tissue Res 318, 121-134.
[27] Braak H, Muller CM, Rub U, Ackermann H, Bratzke H, de Vos RA, Del Tredici K (2006) Pathology associated with sporadic Parkinson’s disease–where does it end? J Neural Transm Suppl, 89-97.
[28] Braak H, Rub U, Sandmann-Keil D, Gai WP, de Vos RA, Jansen Steur EN, Arai K, Braak E (2000) Parkinson’s disease: affection of brain stem nuclei controlling premotor and motor neurons of the somatomotor system. Acta Neuropathol 99, 489-495.
[29] Braak H, Sandmann-Keil D, Gai W, Braak E (1999) Extensive axonal Lewy neurites in Parkinson’s disease: a novel pathological feature revealed by alpha-synuclein immunocytochemistry. Neurosci Lett 265, 67-69.
[30] Braak H, Sastre M, Bohl JR, de Vos RA, Del Tredici K (2007) Parkinson’s disease: lesions in dorsal horn layer I, involvement of parasympathetic and sympathetic pre- and postganglionic neurons. Acta Neuropathol 113, 421-429.
[31] Braak H, Sastre M, Del Tredici K (2007) Development of alpha-synuclein immunoreactive astrocytes in the forebrain parallels stages of intraneuronal pathology in sporadic Parkinson’s disease. Acta Neuropathol 114, 231-241.
[32] Del Tredici K, Braak H (2008) A not entirely benign procedure: progression of Parkinson’s disease. Acta Neuropathol 115, 379-384.
[33] Del Tredici K, Braak H (2012) Spinal cord lesions in sporadic Parkinson’s disease. Acta Neuropathol 124, 643-664.
[34] Del Tredici K, Braak H (2012) Lewy pathology and neurodegeneration in premotor Parkinson’s disease. Mov Disord 27, 597-607.
[35] Del Tredici K, Hawkes CH, Ghebremedhin E, Braak H (2010) Lewy pathology in the submandibular gland of individuals with incidental Lewy body disease and sporadic Parkinson’s disease. Acta Neuropathol 119, 703-713.
[36] Ghebremedhin E, Del Tredici K, Langston JW, Braak H (2009) Diminished tyrosine hydroxylase immunoreactivity in the cardiac conduction system and myocardium in Parkinson’s disease: an anatomical study. Acta Neuropathol 118, 777-784.
[37] Halliday GM, Del Tredici K, Braak H (2006) Critical appraisal of brain pathology staging related to presymptomatic and symptomatic cases of sporadic Parkinson’s disease. J Neural Transm Suppl, 99-103.
[38] Hawkes CH, Del Tredici K, Braak H (2010) A timeline for Parkinson’s disease. Parkinsonism Relat Disord 16, 79-84.
[39] Lionnet A, Leclair-Visonneau L, Neunlist M, Murayama S, Takao M, Adler CH, Derkinderen P, Beach TG (2018) Does Parkinson’s disease start in the gut? Acta Neuropathol 135, 1-12.
[40] Rub U, Del Tredici K, Schultz C, Ghebremedhin E, de Vos RA, Jansen Steur E, Braak H (2002) Parkinson’s disease: the thalamic components of the limbic loop are severely impaired by alpha-synuclein immunopositive inclusion body pathology. Neurobiol Aging 23, 245-254.
[41] Georgiev D, Hamberg K, Hariz M, Forsgren L, Hariz GM (2017) Gender differences in Parkinson’s disease: A clinical perspective. Acta Neurol Scand 136, 570-584.
[42] Fox SH, Katzenschlager R, Lim SY, Barton B, de Bie RMA, Seppi K, Coelho M, Sampaio C, Movement Disorder Society Evidence-Based Medicine C (2018) International Parkinson and movement disorder society evidence-based medicine review: Update on treatments for the motor symptoms of Parkinson’s disease. Mov Disord.
[43] Rughani A, Schwalb JM, Sidiropoulos C, Pilitsis J, Ramirez-Zamora A, Sweet JA, Mittal S, Espay AJ, Martinez JG, Abosch A, Eskandar E, Gross R, Alterman R, Hamani C (2018) Congress of Neurological Surgeons Systematic Review and Evidence-Based Guideline on Subthalamic Nucleus and Globus Pallidus Internus Deep Brain Stimulation for the Treatment of Patients With Parkinson’s Disease: Executive Summary. Neurosurgery.
[44] Santaniello S, Gale JT, Sarma SV (2018) Systems approaches to optimizing deep brain stimulation therapies in Parkinson’s disease. Wiley Interdiscip Rev Syst Biol Med, e1421.
[45] Kempster PA, O’Sullivan SS, Holton JL, Revesz T, Lees AJ (2010) Relationships between age and late progression of Parkinson’s disease: a clinico-pathological study. Brain 133, 1755-1762.
[46] Borrageiro G, Haylett W, Seedat S, Kuivaniemi H, Bardien S (2018) A review of genome-wide transcriptomics studies in Parkinson’s disease. Eur J Neurosci 47, 1-16.
[47] Brohawn DG, O’Brien LC, Bennett JP, Jr. (2016) RNAseq Analyses Identify Tumor Necrosis Factor-Mediated Inflammation as a Major Abnormality in ALS Spinal Cord. PLoS One 11, e0160520.
[48] Bennett J, Keeney P (2017) Micro RNA’s (mirna’s) may help explain expression of multiple genes in Alzheimer’s Frontal Cortex. Journal of Systems and Integrative Neuroscience 3, 1-9.
[49] Polymeropoulos MH (2000) Genetics of Parkinson’s disease. Ann N Y Acad Sci 920, 28-32.
[50] Polymeropoulos MH, Higgins JJ, Golbe LI, Johnson WG, Ide SE, Di Iorio G, Sanges G, Stenroos ES, Pho LT, Schaffer AA, Lazzarini AM, Nussbaum RL, Duvoisin RC (1996) Mapping of a gene for Parkinson’s disease to chromosome 4q21-q23. Science 274, 1197-1199.
[51] Brohawn DG LA, O’Brien LC, Bennett JP Jr <Neuralized human embryonic or induced pluripotential stem cell-derived motor neurons are genetically different from those isolated from human adult cervical spinal cord.pdf>. Journal of Stem Cell Research and Medicine 2(2), 1-5.
[52] Kienhofer S, Musheev MU, Stapf U, Helm M, Schomacher L, Niehrs C, Schafer A (2015) GADD45a physically and functionally interacts with TET1. Differentiation 90, 59-68.
[53] Niehrs C, Schafer A (2012) Active DNA demethylation by Gadd45 and DNA repair. Trends Cell Biol 22, 220-227.
[54] Schuermann D, Weber AR, Schar P (2016) Active DNA demethylation by DNA repair: Facts and uncertainties. DNA Repair (Amst) 44, 92-102.
[55] Li XJ, Li ZF, Wang JJ, Han Z, Liu Z, Liu BG (2017) Effects of microRNA-374 on proliferation, migration, invasion, and apoptosis of human SCC cells by targeting Gadd45a through P53 signaling pathway. Biosci Rep 37.
[56] Moskalev AA, Smit-McBride Z, Shaposhnikov MV, Plyusnina EN, Zhavoronkov A, Budovsky A, Tacutu R, Fraifeld VE (2012) Gadd45 proteins: relevance to aging, longevity and age-related pathologies. Ageing Res Rev 11, 51-66.
[57] Wang XF, Zeng QG, Zeng Y, Man RY, Lu BX, Luo YF (2014) Induction of GADD45alpha protects M17 neuroblastoma cells against MPP*. IUBMB Life 66, 786-792.
[58] Ishizawa M, Akagi D, Yamamoto J, Makishima M (2017) 1alpha,25-Dihydroxyvitamin D3 enhances TRPV6 transcription through p38 MAPK activation and GADD45 expression. J Steroid Biochem Mol Biol 172, 55-61.
[59] Satomi Y (2012) Fucoxanthin induces GADD45A expression and G1 arrest with SAPK/JNK activation in LNCap human prostate cancer cells. Anticancer Res 32, 807-813.
[60] Dunlop EA, Tee AR (2013) The kinase triad, AMPK, mTORC1 and ULK1, maintains energy and nutrient homoeostasis. Biochem Soc Trans 41, 939-943.
[61] Rabanal-Ruiz Y, Korolchuk VI (2018) mTORC1 and Nutrient Homeostasis: The Central Role of the Lysosome. Int J Mol Sci 19.
[62] Rabanal-Ruiz Y, Otten EG, Korolchuk VI (2017) mTORC1 as the main gateway to autophagy. Essays Biochem 61, 565-584.
[63] Tatebe H, Shiozaki K (2017) Evolutionary Conservation of the Components in the TOR Signaling Pathways. Biomolecules 7.
[64] Weichhart T (2018) mTOR as Regulator of Lifespan, Aging, and Cellular Senescence: A Mini-Review. Gerontology 64, 127-134.
[65] Fowler AJ, Moussa CE (2018) Activating Autophagy as a Therapeutic Strategy for Parkinson’s Disease. CNS Drugs 32, 1-11.
[66] Bento CF, Ashkenazi A, Jimenez-Sanchez M, Rubinsztein DC (2016) The Parkinson’s disease-associated genes ATP13A2 and SYT11 regulate autophagy via a common pathway. Nat Commun 7, 11803.
[67] Shin JH, Ko HS, Kang H, Lee Y, Lee YI, Pletinkova O, Troconso JC, Dawson VL, Dawson TM (2011) PARIS (ZNF746) repression of PGC-1alpha contributes to neurodegeneration in Parkinson’s disease. Cell 144, 689-702.
[68] Zheng B, Liao Z, Locascio JJ, Lesniak KA, Roderick SS, Watt ML, Eklund AC, Zhang-James Y, Kim PD, Hauser MA, Grunblatt E, Moran LB, Mandel SA, Riederer P, Miller RM, Federoff HJ, Wullner U, Papapetropoulos S, Youdim MB, Cantuti-Castelvetri I, Young AB, Vance JM, Davis RL, Hedreen JC, Adler CH, Beach TG, Graeber MB, Middleton FA, Rochet JC, Scherzer CR, Global PDGEC (2010) PGC-1alpha, a potential therapeutic target for early intervention in Parkinson’s disease. Sci Transl Med 2, 52ra73.
[69] Bhattacharjee N, Borah A (2016) Oxidative stress and mitochondrial dysfunction are the underlying events of dopaminergic neurodegeneration in homocysteine rat model of Parkinson’s disease. Neurochem Int 101, 48-55.
[70] Jiang P, Dickson DW (2018) Parkinson’s disease: experimental models and reality. Acta Neuropathol 135, 13-32.