Hypothesis and Theory ARTICLE
Evidence for a Pan-Neurodegenerative Disease Response in Huntington's and Parkinson's Disease Expression Profiles
- 1Bioinformatics Program, Boston University, Boston, MA, United States
- 2Department of Neurology, Boston University, Boston, MA, United States
- 3Biostatistics, Boston University School of Public Health, Boston, MA, United States
Huntington's and Parkinson's Diseases (HD and PD) are neurodegenerative disorders that share some pathological features but are disparate in others. For example, while both diseases are marked by aberrant protein aggregation in the brain, the specific proteins that aggregate and types of neurons affected differ. A better understanding of the molecular similarities and differences between these two diseases may lead to a more complete mechanistic picture of both the individual diseases and the neurodegenerative process in general. We sought to characterize the common transcriptional signature of HD and PD as well as genes uniquely implicated in each of these diseases using mRNA-Seq data from post mortem human brains in comparison to neuropathologically normal controls. The enriched biological pathways implicated by HD differentially expressed genes show remarkable consistency with those for PD differentially expressed genes and implicate the common biological processes of neuroinflammation, apoptosis, transcriptional dysregulation, and neuron-associated functions. Comparison of the differentially expressed (DE) genes highlights a set of consistently altered genes that span both diseases. In particular, processes involving nuclear factor kappa-light-chain-enhancer of activated B cells (NFkB) and transcription factor cAMP response element-binding protein (CREB) are the most prominent among the genes common to HD and PD. When the combined HD and PD data are compared to controls, relatively few additional biological processes emerge as significantly enriched, suggesting that most pathways are independently seen within each disorder. Despite showing comparable numbers of DE genes, DE genes unique to HD are enriched in far more coherent biological processes than the DE genes unique to PD, suggesting that PD may represent a more heterogeneous disorder. The complexity of the biological processes implicated by this analysis provides impetus for the development of better experimental models to validate the results.
Transcriptional dysregulation has been observed in both Huntington's disease (HD) and Parkinson's disease (PD) (Cha, 2007; Elstner et al., 2011). Transcription, neuroinflammation, and developmental processes have been shown to be dysregulated in the brains of HD individuals (Labadorf et al., 2015), while inflammation and mitochondrial dysfunction were observed to be altered in the brains of PD individuals (Dumitriu et al., 2012). However, a systematic comparison of the transcriptional signatures of HD and PD has not been performed to date, and those genes and biological processes common to both diseases, if any, remain to be determined. To address this question, we sought to identify genes that are consistently differentially expressed (DE) in the post-mortem brains of HD and PD human subjects compared to neuropathologically normal control brains using mRNA-Seq. We hypothesize that common altered genes and pathways in HD and PD will elucidate the mechanistic underpinnings of the neurodegenerative process.
This study presents the results of a comparison of DE genes for each of HD and PD vs. controls (C) analyzed separately. In addition, in order to identify consistent effects with lower effect size across diseases, an analysis was performed where the HD and PD datasets are concatenated as a single category, neurodegenerative disease (ND), and compared with C. DE genes are determined using a tailored form of logistic regression as previously described (Choi et al., 2017), which improves control of type I errors when compared to negative binomial based DE detection techniques.
Materials and Methods
This study has been designated exempt (Protocol # H-28974) by the Boston University School of Medicine Institutional Review Board, as no human subjects were studied and all data are derived from post-mortem human brain specimens.
Sample Collection, Processing, and Sequencing
The HD, PD, and C samples used in this study are those previously described in our past work (Dumitriu et al., 2012; Labadorf et al., 2015). Nine additional HD brain samples were included in this study beyond those in Labadorf et al. (2015), including two HD gene positive asymptomatic individuals, obtained from the Harvard Brain Tissue Resource Center. All samples underwent the same tissue dissection and RNA extraction sample preparation protocol performed by the same individuals. Briefly, RNA was extracted from the prefrontal cortex of post-mortem brains of HD and PD subjects, as well as neuropathologically normal controls. RNA was poly-A selected and subjected to mRNA sequencing on the Illumina HiSeq 2000 platform. Sample statistics are contained in Table 1. See Labadorf et al. (2015) and Dumitriu et al. (2012) for more detailed information about sample preparation.
mRNA-Seq Data Processing
Each FASTQ file containing mRNA sequences was first trimmed for sequence quality using the sickle software package (Joshi and Fass, 2011) with default arguments. The resulting short reads were aligned against the hg38 build of the human reference genome using the STAR aligner v2.4.0h1 (Dobin et al., 2013) with permissive multimapping parameters (200 maximum alignments—outFilterMultimapNmax 200) and otherwise parameter values suggested in the STAR manual. Multimapped reads were assigned unique alignment locations using the ORMAN software package (Dao et al., 2014). Aligned reads were counted against GENCODE v21 annotation (Harrow et al., 2012) using the HTSeq package v0.6.1p1 (Anders et al., 2014). Read counts for all samples were normalized using the method described in the DESeq2 package v1.10.1 (Love et al., 2014) and outlier counts were trimmed using the strategy as previously described (Labadorf et al., 2015). Since the original mRNAs were poly-A selected, only genes with biotypes known to be polyadenylated (i.e., protein_coding, lincRNA, processed_transcript, sense_intronic, sense_overlaping, IG_V_gene, IG_D_gene, IG_J_gene, IG_C_gene, TR_V_gene, TR_D_gene, TR_J_gene, and TR_C_gene) as annotated by Ensembl BioMart (Kinsella et al., 2011) downloaded on May 27th, 2015. To avoid spurious results due to low abundance, genes were further filtered if more than half of the ND or C samples had zero counts.
Differential Expression and Assessment of Batch Effects
DE genes were identified using Firth's logistic (FL) regression (Firth, 1993; Heinze and Schemper, 2002) applied to mRNA-Seq data as previously described (Choi et al., 2017). Logistic regression is a well-established statistical model that can identify consistent differences in one or more variables between two groups of samples. In contrast to negative binomial regression models like edgeR (Robinson et al., 2010) and DESeq2 (Love et al., 2014), this method models a binomial status variable (e.g., case vs. control) as a function of gene counts and any other potentially confounding variables (RIN value, PMI, etc.). Classical logistic regression has historically not been used to determine DE genes because of the so-called “complete separation” problem, where model parameter estimation fails when there is perfect or nearly perfect separation of a predictor with respect to a binomial variable (e.g., one condition has extremely high read counts and the other has very low read counts). FL regression addresses this issue by using a modified likelihood function to enable reliable parameter estimation, and has other statistical advantages with respect to type I error rates and power. Note the DE statistic from FL regression is log odds ratio (LOR) of case vs. controls, that is, positive LOR indicates greater mRNA abundance in case and negative LOR indicates greater abundance in control. All reported p-values are Benjamini-Hochberg (BH) adjusted (Benjamini and Hochberg, 1995) unless noted otherwise. See Choi et al. (2017) for further information on this method applied to mRNA-Seq data.
The data sets in this study were sequenced in five separate batches. To evaluate whether there was evidence of systemic batch effects confounding the identification of DE genes, we ran separate Firth DE models with and without a categorical variable representing batches and compared the beta estimates using Spearman correlation. Beta estimate ranks were highly correlated for HD vs. C (ρ = 0.84, p ⋘ 0.001), PD vs. C (ρ = 0.99, p ⋘ 0.001), and ND vs. C (ρ = 0.97, p ⋘ 0.001). We therefore did not include a batch variable in the DE models.
Identification of ND DE Genes and Enriched Genesets
DE genes were identified as those with BH adjusted p-values < 0.01 from the Firth's logistic regression models of HD vs. C, PD vs. C, and ND vs. C models, yielding three independent DE gene lists. Read counts for each gene were scaled to have a mean of zero and standard deviation of one to obtain standardized regression coefficients, which makes coefficients comparable across genes. All controls were used in each model. Gene set enrichment analysis was performed on each gene list ranked by read counts beta coefficient using the GSEA (Subramanian et al., 2005) implementation in the DOSE software package (Yu et al., 2015) against the MsigDB Canonical Pathway (C2) geneset database (Subramanian et al., 2005). GSEA enrichment was computed using the complete list of genes irrespective of significance ranked by standardized beta coefficient of the count variable. The robust rank aggregation (RRA) algorithm (Kolde et al., 2012) was used to identify individual genes that were consistently altered across these gene lists. Briefly, RRA is a probabilistic, non-parametric, rank-based method for detecting genes ranked consistently better than expected under the null hypothesis of uncorrelated inputs in the presence of noise and outliers. The genes identified as most significant by RRA are the most likely to be implicated in the general ND phenotype.
In addition to producing independent HD and PD DE gene lists, we sought to functionally characterize the genes that are uniquely significant to each disease as well as those in common. To accomplish this, the DE genes from HD and PD were intersected, partitioning the genes into HD-specific, PD-specific, and DE genes common to the two gene lists. Each of these partitioned gene lists were then subjected to gene set enrichment on the MsigDB Canonical Pathway (C2), Transcription Factor Targets (C3), and Gene Ontology (C5) gene set databases (Subramanian et al., 2005) using a hypergeometric test.
Firth's logistic (FL) regression identified 2,427, 1,949, and 4,843 significantly DE genes for HD, PD, and ND, respectively, at q-value < 0.01. Gene set enrichment analysis of MsigDB C2 gene sets identified 226, 199, and 250 gene sets significantly enriched at q-value < 0.05 for HD, PD, and ND, respectively. Due to the large number of DE genes in each dataset, we focus exclusively on the GSEA enrichment results here. Complete DE gene list and gene set enrichment statistics for HD, PD, and ND are included in Supplemental File 1.
There was a high degree of overlap between the significantly enriched gene sets of HD and PD. 145 gene sets were significantly enriched in both DE gene lists, while 81 and 54 gene sets were uniquely significant in HD and PD, respectively. When a pathway was enriched in more than one list, the pathway was always, without exception, enriched in the same direction, either positively (genes are more abundant in disease) or negatively (genes are less abundant in disease). There were 24 gene sets uniquely significant in ND. Figure 1 depicts the enriched gene sets grouped by high-level biological category for HD, PD, and ND.
Figure 1. (A) Significantly enriched MsigDB C2 Canonical Pathway gene sets for HD, PD, and ND identified by GSEA. Each colored ring segment corresponds to a single enriched gene set. Red (outer), green (middle), and blue (inner) segmented rings indicate whether the HD, PD, or ND DE gene lists, respectively, were significantly enriched for the gene set. Dark and light colored segments indicate up and down regulation (positive, negative GSEA normalized enrichment score), respectively. Black ring around exterior groups gene sets into high level categories as indicated by the two letter code. Gene set name is listed in interior of rings. A tabular form of the data underlying the figure is in Supplemental File 1. (B) Venn diagram illustrating overlap of significantly enriched gene sets for HD, PD, and ND. All but 24 of the ND significant gene sets were significantly enriched in either HD, PD, or both.
We make several observations of Figure 1A. First, the plurality of enriched gene sets across all three data sets are related to immune processes (IM) and are with few exceptions positively enriched. Pathways related to neuronal processes (NE) are largely negatively enriched and there is a subset of these gene sets that are exclusively enriched in HD. With the exception of DNA damage, all remaining biological categories are represented for both HD and PD. DNA damage related pathways (DN) are unique to the PD dataset and are negatively enriched. Multiple apoptosis (AP), developmental (DE), transcription/translation (TR), and transcription factor target (TF) gene sets are also enriched in all three gene lists.
There were 83 gene sets that did not fit cleanly into a single category (OT), which notably include pathways related to endocytosis, signaling, cellular adhesion and extracellular matrix, glycans, and metabolism. A tabular form of the data underlying Figure 1A is in Supplemental File 1.
RRA identified 1,353 genes with a score <0.01. The top ten genes identified by RRA as most highly ranked across all three gene lists are reported in Table 2. The rank of each gene in the individual gene lists are also reported in the table, showing that most genes are relatively highly ranked across all three studies as expected. The list of all significant genes identified by RRA analysis is included in Supplemental File 1.
The most consistently ranked gene is RP1-93H18.7 (ENSG00000272403.1), a lncRNA, which was removed from Ensembl starting at version GRCh38.p2, but shows consistent transcription in these data. The annotation of this gene was likely removed because it was found to be the 3′ UTR end of one transcript (ENST00000452085.7 as of human genome release GRCh38.p10) of the gene dermatan sulfate epimerase (DSE). DSE is also DE in both HD and PD, is involved in embryonic development (Habicher et al., 2015; Stachtea et al., 2015), and has been related to the immune response in cancer patients (Mizukoshi et al., 2012). The significance of ENSG00000272403 is likely attributable to the differential expression of DSE itself. See Supplemental File 2 for a more thorough investigation of this gene.
Deficiencies in the second ranked gene, SPR (sepiapterin reductase), have been linked to DOPA-responsive dystonia (Wijemanne and Jankovic, 2015) and previously implicated in PD (Tobin et al., 2007). The third gene, DDIT4 (DNA-Damage-Inducible Transcript 4), is a multi-functional gene which, via its inhibition of the mammalian target of rapamycin complex 1 (mTORC1), regulates in cell growth, proliferation, and survival (Dennis et al., 2013), controls p53/TP53-mediated apoptosis in response to DNA damage (Vadysirisack et al., 2011; Cam et al., 2014), and plays a role in neurodegeneration, neuronal death and differentiation, and neuron migration during embryonic brain development (Malagelada et al., 2011; Canal et al., 2014; Ota et al., 2014; Romaní-Aumedes et al., 2014). TRIP10 (thyroid hormone receptor interactor 10), another multi-functional gene, is involved in insulin signaling (Chang et al., 2013), endocytosis (Feng et al., 2010), and structures specific to monocyte-derived cells (Linder et al., 2000). TNFRSF10D (Tumor Necrosis Factor Receptor Superfamily, Member 10d, Decoy With Truncated Death Domain) inhibits certain types of apoptosis and may play a role in NfkB pathway (Degli-Esposti et al., 1997).
Figure 2 illustrates the differences in normalized counts for the top genes identified by RRA. With the exception of (12) PITX1, which is driven entirely by HD, all top genes demonstrate substantial differences between both disease conditions and control.
Figure 2. Box plots of normalized counts for top RRA genes split by condition, RRA rank is in parenthesis. Whiskers extend to 25th and 75th percentile counts, white bars are median counts. With the exception of (12) PITX1, which is driven entirely by HD, all top genes demonstrate substantial differences between both disease conditions and control.
Finally, we examined the significant DE genes from HD and PD for intersection. Figure 3 illustrates the overlap of DE genes between diseases and describes gene set enrichment results for the intersection. Figure 4 contains the enrichment results for the HD unique genes.
Figure 3. (A) Venn diagram of HD and PD DE gene list intersection for DE genes adjusted p < 0.01. (B) Bar chart indicating number of MsigDB C2 Canonical Pathway (CP), C3 miRNA Targets (miR), C3 Transcription Factor Targets (TF), and C5 Gene Ontology (GO) gene sets enriched for the HD unique (HD\PD), intersection (HD n PD), and PD unique (PD\HD) genes. For HD\PD enrichment, 17 redundant or uninformative GO gene sets and 7 TF gene sets for motifs with unknown transcription factors were omitted from the figure results but are included in Supplemental File 1. (C) Gene sets enriched for the intersection genes (HD n PD). Adjusted p-values are listed beside gene set name and the originating gene set (CP, miR, TF, or GO) are indicated by color. Gene sets that are grouped into boxes share more than 20% of their DE genes and are therefore listed together.
Figure 4. Enriched gene sets for the HD unique (HD\PD) genes from Figure 3A and reported similarly as in Figure 3C. Note 17 redundant or uninformative GO gene sets and 7 TF gene sets for motifs with unknown transcription factors were omitted from the figure results but are included in Supplemental File 1.
HD specific enrichments are shown in Figure 4, where a broad spectrum of biological processes are implicated in the HD-unique DE genes. The most striking enriched gene set is KEGG_RIBOSOME, with many other related gene sets involved in translation and molecular metabolism similarly enriched. Multiple gene sets that share >20% of their DE genes are associated with Jun Proto-Oncogene (AP1), BTB And CNC Homology 1, Basic Leucine Zipper Transcription Factor 1 and 2 (BACH1, BACH2), and NRF2/TCF11 are also implicated. Other strongly implicated biological processes are ion channel activity, plasma membrane and signaling, apoptosis, immune system and inflammatory processes, developmental genes, neuron-related signaling pathways, many transcription factors, and two families of miRNAs.
RNA-sequencing was performed in 29 HD, 29 PD, and 49 control prefrontal cortex (BA9) samples to assess the common and unique transcriptional profiles of these two protein aggregation related neurodegenerative diseases. Firth's logistic (FL) regression identified 2,427, 1,949, and 4,843 significantly DE genes for HD, PD, and ND, respectively, at q-value < 0.01. Gene set enrichment analysis of MsigDB C2 gene sets identified 226, 199, and 250 gene sets significantly enriched at q-value < 0.05 for HD, PD, and ND, respectively. Gene sets related to immune processes and inflammatory pathways were highly enriched for both diseases. Notably, Figure 1 shows that the overwhelming majority of enriched biological pathways are common to both diseases and that they are invariably perturbed in the same direction in both diseases. To the authors knowledge, this study presents the first comprehensive comparative analysis of DE gene expression from HD, PD, and ND in post-mortem human brains assessed with mRNA-Seq. In one previous study, Capurro et al analyzed HD and PD microarray data using a cell-type deconvolution method to identify cell-type specific differences in gene expression between cases and controls for both diseases (Capurro et al., 2014). Only one gene from the study, doublecortin-like kinase 1 (DCLK1), was found to be differentially expressed in both HD and PD cortex, and this gene also appears as DE in common between diseases in the analysis presented here.
The comparison of HD and PD in particular is motivated by the fact that these diseases can be viewed as mirror-images of each other. GABAergic medium spiny interneurons, which compose most of the neurons in the striatum and selectively die in HD but are spared in PD, project directly into the substantia nigra and coordinate motor activity throughout the brain via dopamine-induced signaling (Yager et al., 2015). Dopaminergic neurons in the substantia nigra, on the other hand, which also are important in coordinating motor activity as well as arousal, reinforcement, and reward (Schultz, 2007), selectively degenerate in PD but are spared in HD. It was observed in a study of 523 HD subjects that the incidence of PD in this cohort was lower than that of the general population, though both HD and PD individuals develop Alzheimer's disease at the same rate (Hadzi et al., 2012), suggesting the selective death of medium spiny neurons might be neuroprotective of dopaminergic neuron death. Given the intimate neurological link between the affected neurons in HD and PD, and the mutual exclusivity of their degeneration, this comparison poses a very interesting contrast to identify common responses to neurodegeneration that are not confounded by neuron type. Unfortunately, a direct comparison of neurons in these regions of post-mortem human brains is not possible, precisely due to this mirror-image pathology. The choice of the BA9 brain region is motivated by the fact that, due to degeneration, the affected neurons are largely missing from the striatum and substantia nigra in HD and PD, respectively, whereas BA9 is generally free of involvement in both diseases (Braak et al., 2003; Halliday et al., 2006; Hadzi et al., 2012). Because the primarily affected neurons in HD and PD do not exist in BA9, the biological processes implicated by this analysis may likely represent the response to, rather than direct cause of, the respective diseases. Nonetheless, the remarkable consistency between HD and PD observed in this analysis points to important mechanisms that further our understanding of neurodegenerative disease as a general process.
The biological processes implicated by DE gene lists identified from each condition separately are compellingly similar. From Figure 1, we see that the majority of enriched biological pathways are common and that they are invariably perturbed in the same direction in both diseases. Furthermore, combining HD and PD data into an ND condition does not yield significantly more novel biological insights. This remarkable consistency between the pathway enrichment results suggest that the underlying molecular responses to neurodegeneration in HD and PD may be more similar than they are different, despite their different pathological underpinnings. Of particular significance is the strong positive enrichment of immune and inflammatory pathways, which have been convincingly implicated in both diseases separately (Dobbs et al., 1999; Jenner, 2003; Kwan et al., 2012; Dexter and Jenner, 2013; Ellrichmann et al., 2013; Crotti et al., 2014; Allen Reish and Standaert, 2015; Labadorf et al., 2015), but the compelling similarity of these signatures between HD and PD revealed by this analysis has not been illustrated to date.
The negative regulation of neuron-related pathways is also noteworthy, since the BA9 brain region, from which these samples are derived, is not known to be heavily involved in either of these diseases. Despite the lack of clear and consistent neurodegeneration in this brain region, the widespread biological pathways shown to be affected in this analysis strongly suggest neurons in BA9 do indeed experience a common set of effects in the neuropathology for HD and PD.
Many of the individual genes identified by RRA as most consistently different in HD, PD, and ND have previously been the focus of studies in neurodegeneration. The second highest ranked gene SPR has been the focus of significant study in PD and is related to the PARK3 gene locus (Tobin et al., 2007; Sharma et al., 2011), but has not been previously implicated in HD. Inhibition of DNA-damage inducible transcript 4 (DDIT4/RTP801/REDD1) has been associated with neuroprotection in PD models and patients (Malagelada et al., 2011) and is involved with mutant Huntingtin-induced cell death (Martín-Flores et al., 2015). Thyroid hormone receptor interactor 10 (TRIP10) has been shown to interact directly with mutant huntingtin (Holbert et al., 2003), and while it is not known to play a role in PD pathology, its elevated mRNA abundance in these PD samples suggest it may indeed be implicated. Other top genes have also been implicated in neurodegeneration: tumor necrosis factor receptor superfamily 10D (TNFRSF10D) (López-Gómez et al., 2011; Frenkel, 2015), protein arginine methyltransferase 6 (PRMT6) (Scaramuzzino et al., 2015), and toll-like receptor 5 (TLR5) (Arroyo et al., 2011). Further investigation of this list of genes is likely to yield novel insights into the mechanisms of neurodegeneration.
The intersection of DE genes between HD and PD also affords important insight into genes relevant to fundamental neurodegenerative processes. Most notably, pathways related to NFkB and transcriptional targets of CREB are prominent in the enrichment results. The NFkB pathway is prominent in both HD (Marcora and Kennedy, 2010; Träger et al., 2014) and PD pathology (Ghosh et al., 2007; Flood et al., 2011) through its central role in inflammatory signaling. CREB is directly targeted by brain derived-neurotrophic factor (BDNF) (Pizzorusso et al., 2000), an important trophic factor in the brain. Both BDNF (Zuccato et al., 2008), and CREB (Obrietan and Hoyt, 2004; Choi et al., 2017) have been directly implicated in HD pathology, while CREB is also believed to play a critical role in dopamine receptor-mediated nuclear signaling (Andersson et al., 2001), and disruption of the protein's function leads to neurodegeneration (Mantamadiotis et al., 2002; Devi and Ohno, 2014). The specific inflammation-related gene sets (HSF1 transcription factor targets, CXCR4, IL12) suggests there is some specificity in the aspects of the pan-neurodegenerative neuroimmune response. Recent studies in both HD and PD have focused on the role of insulin sensitivity and metabolism in patients (Block et al., 2010; Aviles-Olmos et al., 2013; Russo et al., 2013), supporting the role of insulin synthesis as an enriched biological pathway in the common gene list. While the enrichment of apoptosis-related pathways was not surprising, pathways related to angiopoietin, ephrin, and axon guidance suggest that biological processes related to neuronal plasticity are active in both of these diseases and may even indicate that neuroprotective or neuroregenerative processes are a component of the neurodegenerative response.
These data also point to compelling differences between HD and PD. Interestingly, two groups of genes, DNA damage and repair and tRNA related processes, seem to be uniquely perturbed and negatively enriched in PD. The DNA damage and repair gene set enrichment may be a reflection of mitochondrial DNA damage. In PD, dopaminergic neurons of the substantia nigra (though not cortical neurons) were found to be particularly vulnerable to mitochondrial DNA damage (Sanders et al., 2014), and Lewy body pathology, the histological hallmark of PD, is associated with mitochondrial DNA damage (Müller et al., 2013). More generally, mitochondrial DNA damage and oxidative stress are associated with several neurodegenerative diseases including PD, Alzheimer's disease (Moreira et al., 2010), and ALS (Coppedè, 2011). There is evidence supporting the involvement of aminoacyl tRNA synthetases in neurological disease, including ALS, leukoencephalopathy, and PD (Park et al., 2008).
In HD, a number of uniquely perturbed gene sets related to glycan biosynthesis and metabolism are negatively regulated, and these pathways have not been previously implicated in HD. The 1687 HD-unique DE genes are enriched for many gene sets across a broad spectrum of biological processes, including mRNA and protein metabolism, ion channel activity, signaling and kinase activity, apoptosis, immune response, and development. Other, more specific gene sets related to a large number of transcription factors further support the observation of transcriptional dysregulation in HD (Cha, 2000; Labadorf et al., 2015). The specificity of these enriched TF gene sets is quite striking, as the targeted DE genes appear to be largely disjoint between them, suggesting potential, specific causes of the dysregulated transcriptional effects in HD. The enrichment of two miRNA families are also particularly relevant in light of recent reports of miRNA dysregulation in HD (Hoss et al., 2014, 2015).
It is interesting to note the disparity in enrichment between the HD and PD unique DE genes. Though the numbers of unique DE genes are comparable, the large number of enriched gene sets in HD stands in sharp contrast to the almost total absence of enrichment in PD. This result implies that the DE genes in HD are more consistently related to one another than in PD. One possible, and potentially important explanation for this is that HD is a much more homogeneous disease than PD. It is well established that PD has a significant sporadic component (Lesage and Brice, 2012), caused by a combination of genetic and environmental factors. The relative heterogeneity of PD may make finding consistently effective treatments difficult, and the absence of biological enrichment in specific pathways, other than those common to both diseases, from this analysis may be a reflection of an underlying molecular basis for this effect. It may be that, given sufficient sample size, coherent subgroups of patients may be identified by examining patterns in their gene expression using datasets such as those analyzed here. On the other hand, despite extensive molecular characterization of HD, effective, widely available therapies for HD have remained elusive despite the relative homogeneity of the disease process among HD patients.
Despite the strong literature support of the genes and pathways implicated by this analysis, the findings are purely computational and, therefore, remain hypothetical. Experimentally validating these results poses significant challenges due to the complexity of the pathways implicated in whole tissue from human brain in many respects. First, a fundamental finding of this study is that thousands of genes are implicated between the two diseases, and no one individual gene is immediately apparent as a driver justifying specific follow-up study. Even if an expanded set of genes were chosen while still making wet lab interrogation tractable, the fundamental questions raised by the study would still not be meaningfully addressed. Second, the most compelling result—that these diseases share nearly identical inflammatory signatures—was discovered on the pathway level, not the gene level, suggesting that any functional study would have to modulate entire pathways to make a comparison meaningful. Third, the model systems currently available for HD and PD are very poor with respect to the human disease presentation, in particular to the inflammatory pathways implicated here. For instance, HD mouse models require substantially higher CAG expansions than are seen in humans to show any symptoms, and no mouse studies to date have implicated the strong inflammatory response observed in human tissues by this and other studies. Fourth, when using cell culture systems, since the immune response involves interaction between different cell types, monoculture in vitro systems cannot recapitulate this phenotype. While co-culture of neurons and glia is possible, these systems are not advanced to the point where it is possible to invoke anything like the complex response implicated by the human data. And last, related to this last point, the contribution of age to neuroinflammation is not known, and since these individuals are mostly 50 years old and greater, we would have to find an appropriate model to assess age vs. acute inflammatory vs. chronic inflammatory effects to properly interpret a meaningful experiment. Thus, the available tools as of this study are generally inadequate to meaningfully interrogate these findings in a wet lab setting without substantial improvement and significant effort. The results here are thus proposed as important and potentially fruitful avenues of future study, both in the development of model systems to aid in understanding the mechanism of the neuroinflammatory process in ND and in better understanding the neurodegenerated human brain itself.
The data quality of these datasets has been validated in their original studies. A set of genes across the dynamic range of the sequencing results in 20 of the 29 HD datasets when compared to all 49 C were successfully validated by qPCR (Labadorf et al., 2015). In the PD study, an earlier microarray dataset using many of the same samples was conducted to assess gene expression differences between 19 PD and 24 C and was directly compared to the sequencing data. There was a high level of agreement between the DE genes found in both analyses, which are orthogonal platforms and were generated using different RNA extractions of the same brains (Dumitriu et al., 2012). These previously published results show the validity of the data, which is further enhanced by the agreement between these previous studies and this one, despite using a different analytical methodology.
Nonetheless, these findings have important implications on our understanding of the neurodegenerative disease process. The significant involvement of the inflammatory pathways in both diseases in an area not thought to be directly involved in disease pathogenesis suggests the response to neurodegeneration is widespread throughout the brain. NFkB in particular appears to be a major player, which is well supported in the HD and PD literature. It is unclear whether the neuroinflammatory response is protective, deleterious, or both from these data, but investigation into the role these processes play, and the potential therapeutic value of modulating them, should be made a high priority.
Raw FASTQ reads for these samples are available under GEO accession numbers GSE64810 and GSE68719. The code used to perform the statistical analysis and generate figures is available at https://bitbucket.org/adamlabadorf/hdvpd_mrnaseq.
AL: designed the implementation of the experiments, wrote all analysis code, and prepared the manuscript; SC: provided guidance and critical feedback on the use of Firth logistic regression; RM: conceived of the overall study, made available and processed the samples, and provided critical study and manuscript feedback.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supported by grants from US National Institutes of Health (R01-S076843), Characterization of the Role of Cyclin G-associated Kinase in Parkinson Disease (R01-NS073947), Epigenetic Markers in Huntington's Disease Brain, (R01-NS088538) An IPSc based platform for functionally assessing genetic and environmental Risk in PD, (U24-NS072026) National Brain and Tissue Resource for Parkinson's Disease and Related Disorders, and the Jerry McDonald Huntington Disease Research Fund.
We would like to acknowledge the National Brain and Tissue Resource for Parkinson's Disease and Related Disorders at Banner Sun Health Research Institute (NS072026), Sun City, Arizona, the Harvard Brain Tissue Resource Center McLean Hospital, Belmont, Massachusetts, and the Human Brain and Spinal Fluid Resource Center VA, West Los Angeles Healthcare Center, California for providing the brain samples used in these studies.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnmol.2017.00430/full#supplementary-material
Andersson, M., Konradi, C., and Cenci, M. A. (2001). cAMP response element-binding protein is required for dopamine-dependent gene expression in the intact but not the dopamine-denervated striatum. J. Neurosci. 21, 9930–9943.
Arroyo, D. S., Soria, J. A., Gaviglio, E. A., Rodriguez-Galan, M. C., and Iribarren, P. (2011). Toll-like receptors are key players in neurodegeneration. Int. Immunopharmacol. 11, 1415–1421. doi: 10.1016/j.intimp.2011.05.006
Block, R. C., Dorsey, E. R., Beck, C. A., Brenna, J. T., and Shoulson, I. (2010). Altered cholesterol and fatty acid metabolism in Huntington disease. J. Clin. Lipidol. 4, 17–23. doi: 10.1016/j.jacl.2009.11.003
Braak, H., Del Tredici, K., Rüb, U., de Vos, R. A., Jansen Steur, E. N., and Braak, E. (2003). Staging of brain pathology related to sporadic Parkinson's disease. Neurobiol. Aging 24, 197–211. doi: 10.1016/S0197-4580(02)00065-9
Cam, M., Bid, H. K., Xiao, L., Zambetti, G. P., Houghton, P. J., Cam, H., et al. (2014). p53/TAp63 and AKT regulate mammalian target of rapamycin complex 1 (mTORC1) signaling through two independent parallel pathways in the presence of DNA damage. J. Biol. Chem. 289, 4083–4094. doi: 10.1074/jbc.M113.530303
Canal, M., Romaní-Aumedes, J., Martín-Flores, N., Pérez-Fernández, V., and Malagelada, C. (2014). RTP801/REDD1: a stress coping regulator that turns into a troublemaker in neurodegenerative disorders. Front. Cell. Neurosci. 8:313. doi: 10.3389/fncel.2014.00313
Capurro, A., Bodea, L.-G., Schaefer, P., Luthi-Carter, R., and Perreau, V. M. (2014). Computational deconvolution of genome wide expression data from Parkinson's and Huntington's disease brain tissues using population-specific expression analysis. Front. Neurosci. 8:441. doi: 10.3389/fnins.2014.00441
Choi, S. H., Labadorf, A. T., Myers, R. H., Lunetta, K. L., Dupuis, J., and DeStefano, A. L. (2017). Evaluation of logistic regression models and effect of covariates for case-control study in RNA-Seq analysis. BMC Bioinformatics 18:91. doi: 10.1186/s12859-017-1498-y
Crotti, A., Benner, C., Kerman, B. E., Gosselin, D., Lagier-Tourenne, C., Zuccato, C., et al. (2014). Mutant huntingtin promotes autonomous microglia activation via myeloid lineage-determining factors. Nat. Neurosci. 17, 513–521. doi: 10.1038/nn.3668
Dao, P., Numanagić, I., Lin, Y. Y., Hach, F., Karakoc, E., Donmez, N., et al. (2014). ORMAN: optimal resolution of ambiguous RNA-Seq multimappings in the presence of novel isoforms. Bioinformatics 30, 644–651. doi: 10.1093/bioinformatics/btt591
Degli-Esposti, M. A., Dougall, W. C., Smolak, P. J., Waugh, J. Y., Smith, C. A., Goodwin, R. G., et al. (1997). The novel receptor TRAIL-R4 induces NF-κB and protects against TRAIL-mediated apoptosis, yet retains an incomplete death domain. Immunity 7, 813–820. doi: 10.1016/S1074-7613(00)80399-4
Dennis, M. D., McGhee, N. K., Jefferson, L. S., and Kimball, S. R. (2013). Regulated in DNA damage and development 1 (REDD1) promotes cell survival during serum deprivation by sustaining repression of signaling through the mechanistic target of rapamycin in complex 1 (mTORC1). Cell. Signal. 25, 2709–2716. doi: 10.1016/j.cellsig.2013.08.038
Devi, L., and Ohno, M. (2014). PERK mediates eIF2α phosphorylation responsible for BACE1 elevation, CREB dysfunction and neurodegeneration in a mouse model of alzheimer's disease. Neurobiol. Aging 35, 2272–2281. doi: 10.1016/j.neurobiolaging.2014.04.031
Dobbs, R. J., Charlett, A., Purkiss, A. G., Dobbs, S. M., Weller, C., and Peterson, D. W. (1999). Association of circulating TNF-α and IL-6 with ageing and Parkinsonism. Acta Neurol. Scand. 100, 34–41. doi: 10.1111/j.1600-0404.1999.tb00721.x
Dumitriu, A., Latourelle, J. C., Hadzi, T. C., Pankratz, N., Garza, D., Miller, J. P., et al. (2012). Gene expression profiles in Parkinson disease prefrontal cortex implicate FOXO1 and genes under its transcriptional regulation. PLoS Genet. 8:e1002794. doi: 10.1371/journal.pgen.1002794
Elstner, M., Morris, C. M., Heim, K., Bender, A., Mehta, D., Jaros, E., et al. (2011). Expression analysis of dopaminergic neurons in Parkinson's disease and aging links transcriptional dysregulation of energy metabolism to cell death. Acta Neuropathol. 122, 75–86. doi: 10.1007/s00401-011-0828-9
Feng, Y., Hartig, S. M., Bechill, J. E., Blanchard, E. G., Caudell, E., Corey, S. J., et al. (2010). The cdc42-interacting protein-4 (CIP4) gene knock-out mouse reveals delayed and decreased endocytosis. J. Biol. Chem. 285, 4348–4354. doi: 10.1074/jbc.M109.041038
Flood, P. M., Qian, L., Peterson, L. J., Zhang, F., Shi, J.-S., Gao, H.-M., et al. (2011). Transcriptional factor NF-κB as a target for therapy in Parkinson's disease. Parkinsons Dis. 2011:216298. doi: 10.4061/2011/216298
Ghosh, A., Roy, A., Liu, X., Kordower, J. H., Mufson, E. J., Hartley, D. M., et al. (2007). Selective inhibition of NF-κB activation prevents dopaminergic neuronal loss in a mouse model of Parkinson's disease. Proc. Natl. Acad. Sci. U.S.A. 104, 18754–18759. doi: 10.1073/pnas.0704908104
Habicher, J., Haitina, T., Eriksson, I., Holmborn, K., Dierker, T., Ahlberg, P. E., et al. (2015). Chondroitin/dermatan sulfate modification enzymes in zebrafish development. PLoS ONE 10:e0121957. doi: 10.1371/journal.pone.0121957
Hadzi, T. C., Hendricks, A. E., Latourelle, J. C., Lunetta, K. L., Cupples, L. A., Gillis, T., et al. (2012). Assessment of cortical and striatal involvement in 523 Huntington disease brains. Neurology 79, 1708–1715. doi: 10.1212/WNL.0b013e31826e9a5d
Halliday, G. M., Del Tredici, K., and Braak, H. (2006). Critical appraisal of brain pathology staging related to presymptomatic and symptomatic cases of sporadic Parkinson's disease. J. Neural Transm. Suppl. 70, 99–103. doi: 10.1007/978-3-211-45295-0_16
Harrow, J., Frankish, A., Gonzalez, J. M., Tapanari, E., Diekhans, M., Kokocinski, F., et al. (2012). GENCODE: the reference human genome annotation for the ENCODE project. Genome Res. 22, 1760–1774. doi: 10.1101/gr.135350.111
Holbert, S., Dedeoglu, A., Humbert, S., Saudou, F., Ferrante, R. J., Néri, C., et al. (2003). Cdc42-interacting protein 4 binds to huntingtin: neuropathologic and biological evidence for a role in Huntington's disease. Proc. Natl. Acad. Sci. U.S.A. 100, 2712–2717. doi: 10.1073/pnas.0437967100
Hoss, A. G., Kartha, V. K., Dong, X., Latourelle, J. C., Dumitriu, A., Hadzi, T. C., et al. (2014). MicroRNAs located in the hox gene clusters are implicated in Huntington's disease pathogenesis. PLoS Genet. 10:e1004188. doi: 10.1371/journal.pgen.1004188
Hoss, A. G., Labadorf, A., Latourelle, J. C., Kartha, V. K., Hadzi, T. C., Gusella, J. F., et al. (2015). miR-10b-5p expression in Huntington's disease brain relates to age of onset and the extent of striatal involvement. BMC Med. Genomics 8:10. doi: 10.1186/s12920-015-0083-3
Kinsella, R. J., Kähäri, A., Haider, S., Zamora, J., Proctor, G., Spudich, G., et al. (2011). Ensembl BioMarts: a hub for data retrieval across taxonomic space. Database 2011:bar030. doi: 10.1093/database/bar030
Kwan, W., Träger, U., Davalos, D., Chou, A., Bouchard, J., Andre, R., et al. (2012). Mutant huntingtin impairs immune cell migration in Huntington disease. J. Clin. Invest. 122, 4737–4747. doi: 10.1172/JCI64484
Labadorf, A., Hoss, A. G., Lagomarsino, V., Latourelle, J. C., Hadzi, T. C., Bregu, J., et al. (2015). RNA sequence analysis of human Huntington disease brain reveals an extensive increase in inflammatory and developmental gene expression. PLoS ONE 10:e0143563. doi: 10.1371/journal.pone.0143563
López-Gómez, C., Fernández, O., García-León, J. A., Pinto-Medel, M. J., Oliver-Martos, B., Ortega-Pinazo, J., et al. (2011). TRAIL/TRAIL receptor system and susceptibility to multiple sclerosis. PLoS ONE 6:e21766. doi: 10.1371/journal.pone.0021766
Malagelada, C., López-Toledano, M. A., Willett, R. T., Jin, Z. H., Shelanski, M. L., Greene, L. A., et al. (2011). RTP801/REDD1 regulates the timing of cortical neurogenesis and neuron migration. J. Neurosci. 31, 3186–3196. doi: 10.1523/JNEUROSCI.4011-10.2011
Mantamadiotis, T., Lemberger, T., Bleckmann, S. C., Kern, H., Kretz, O., Martin Villalba, A., et al. (2002). Disruption of CREB function in brain leads to neurodegeneration. Nat. Genet. 31, 47–54. doi: 10.1038/ng882
Marcora, E., and Kennedy, M. B. (2010). The Huntington's disease mutation impairs huntingtin's role in the transport of NF-κB from the synapse to the nucleus. Hum. Mol. Genet. 19, 4373–4384. doi: 10.1093/hmg/ddq358
Martín-Flores, N., Romaní-Aumedes, J., Rué, L., Canal, M., Sanders, P., Straccia, M., et al. (2015). RTP801 is involved in mutant Huntingtin-Induced cell death. Mol. Neurobiol. 53, 2857–2868. doi: 10.1007/s12035-015-9166-6
Mizukoshi, E., Fushimi, K., Arai, K., Yamashita, T., Honda, M., and Kaneko, S. (2012). Expression of chondroitin-glucuronate c5-epimerase and cellular immune responses in patients with hepatocellular carcinoma. Liver Int. 32, 1516–1526. doi: 10.1111/j.1478-3231.2012.02853.x
Moreira, P. I., Carvalho, C., Zhu, X., Smith, M. A., and Perry, G. (2010). Mitochondrial dysfunction is a trigger of alzheimer's disease pathophysiology. Biochim. Biophys. Acta 1802, 2–10. doi: 10.1016/j.bbadis.2009.10.006
Müller, S. K., Bender, A., Laub, C., Högen, T., Schlaudraff, F., Liss, B., et al. (2013). Lewy body pathology is associated with mitochondrial DNA damage in Parkinson's disease. Neurobiol. Aging 34, 2231–2233. doi: 10.1016/j.neurobiolaging.2013.03.016
Ota, K. T., Liu, R. J., Voleti, B., Maldonado-Aviles, J. G., Duric, V., Iwata, M., et al. (2014). REDD1 is essential for stress-induced synaptic loss and depressive behavior. Nat. Med. 20, 531–535. doi: 10.1038/nm.3513
Pizzorusso, T., Ratto, G. M., Putignano, E., and Maffei, L. (2000). Brain-derived neurotrophic factor causes cAMP response element-binding protein phosphorylation in absence of calcium increases in slices and cultured neurons from rat visual cortex. J. Neurosci. 20, 2809–2816.
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Romaní-Aumedes, J., Canal, M., Martín-Flores, N., Sun, X., Pérez-Fernández, V., Wewering, S., et al. (2014). Parkin loss of function contributes to RTP801 elevation and neurodegeneration in Parkinson's disease. Cell Death Dis. 5:e1364. doi: 10.1038/cddis.2014.333
Russo, C. V., Salvatore, E., Saccà, F., Tucci, T., Rinaldi, C., Sorrentino, P., et al. (2013). Insulin sensitivity and early-phase insulin secretion in normoglycemic Huntington's disease patients. J. Huntingtons Dis. 2, 501–507. doi: 10.3233/JHD-130078
Sanders, L. H., McCoy, J., Hu, X., Mastroberardino, P. G., Dickinson, B. C., Chang, C. J., et al. (2014). Mitochondrial DNA damage: molecular marker of vulnerable nigral neurons in Parkinson's disease. Neurobiol. Dis. 70, 214–223. doi: 10.1016/j.nbd.2014.06.014
Scaramuzzino, C., Casci, I., Parodi, S., Lievens, P. M., Polanco, M. J., Milioto, C., et al. (2015). Protein arginine methyltransferase 6 enhances polyglutamine-expanded androgen receptor function and toxicity in spinal and bulbar muscular atrophy. Neuron 85, 88–100. doi: 10.1016/j.neuron.2014.12.031
Sharma, M., Maraganore, D. M., Ioannidis, J. P., Riess, O., Aasly, J. O., Annesi, G., et al. (2011). Role of sepiapterin reductase gene at the PARK3 locus in Parkinson's disease. Neurobiol. Aging 32, 2108.e1–5. doi: 10.1016/j.neurobiolaging.2011.05.024
Stachtea, X. N., Tykesson, E., van Kuppevelt, T. H., Feinstein, R., Malmström, A., Reijmers, R. M., et al. (2015). Dermatan Sulfate-Free mice display embryological defects and are neonatal lethal despite normal lymphoid and Non-Lymphoid organogenesis. PLoS ONE 10:e0140279. doi: 10.1371/journal.pone.0140279
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102
Tobin, J. E., Cui, J., Wilk, J. B., Latourelle, J. C., Laramie, J. M., McKee, A. C., et al. (2007). Sepiapterin reductase expression is increased in Parkinson's disease brain tissue. Brain Res. 1139, 42–47. doi: 10.1016/j.brainres.2007.01.001
Träger, U., Andre, R., Lahiri, N., Magnusson-Lind, A., Weiss, A., Grueninger, S., et al. (2014). HTT-lowering reverses Huntington's disease immune dysfunction caused by NFκB pathway dysregulation. Brain 137, 819–833. doi: 10.1093/brain/awt355
Vadysirisack, D. D., Baenke, F., Ory, B., Lei, K., and Ellisen, L. W. (2011). Feedback control of p53 translation by REDD1 and mTORC1 limits the p53-dependent DNA damage response. Mol. Cell. Biol. 31, 4356–4365. doi: 10.1128/MCB.05541-11
Yu, G., Wang, L.-G., Yan, G.-R., and He, Q.-Y. (2015). DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics 31, 608–609. doi: 10.1093/bioinformatics/btu684
Zuccato, C., Marullo, M., Conforti, P., MacDonald, M. E., Tartari, M., Cattaneo, E., et al. (2008). RESEARCH ARTICLE: systematic assessment of BDNF and its receptor levels in human cortices affected by huntington's disease. Brain Pathol. 18, 225–238. doi: 10.1111/j.1750-3639.2007.00111.x
Keywords: Huntington disease, Parkinson disease, mRNA-seq, neuroinflammation, bioinformatics and computational biology, neurodegenerative diseases, immune system
Citation: Labadorf A, Choi SH and Myers RH (2018) Evidence for a Pan-Neurodegenerative Disease Response in Huntington's and Parkinson's Disease Expression Profiles. Front. Mol. Neurosci. 10:430. doi: 10.3389/fnmol.2017.00430
Received: 14 June 2017; Accepted: 12 December 2017;
Published: 11 January 2018.
Edited by:Hermona Soreq, Hebrew University of Jerusalem, Israel
Reviewed by:Alessandro Guffanti, Genomnia srl, Italy
Maciej Maurycy Lalowski, University of Helsinki, Finland
Copyright © 2018 Labadorf, Choi and Myers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Adam Labadorf, firstname.lastname@example.org