- Research article
- Open Access
Gene expression within the periaqueductal gray is linked to vocal behavior and early-onset parkinsonism in Pink1 knockout rats
BMC Genomics volume 21, Article number: 625 (2020)
Parkinson’s disease (PD) is a degenerative disease with early-stage pathology hypothesized to manifest in brainstem regions. Vocal deficits, including soft, monotone speech, result in significant clinical and quality of life issues and are present in 90% of PD patients; yet the underlying pathology mediating these significant voice deficits is unknown. The Pink1−/− rat is a valid model of early-onset PD that presents with analogous vocal communication deficits. Previous work shows abnormal α-synuclein protein aggregation in the periaqueductal gray (PAG), a brain region critical and necessary to the modulation of mammalian vocal behavior. In this study, we used high-throughput RNA sequencing to examine gene expression within the PAG of both male and female Pink1−/− rats as compared to age-matched wildtype controls. We used a bioinformatic approach to (1) test the hypothesis that loss of Pink1 in the PAG will influence the differential expression of genes that interact with Pink1, (2) highlight other key genes that relate to this type of Mendelian PD, and (3) catalog molecular targets that may be important for the production of rat vocalizations.
Knockout of the Pink1 gene resulted in differentially expressed genes for both male and female rats that also mapped to human PD datasets. Pathway analysis highlighted several significant metabolic pathways. Weighted gene co-expression network analysis (WGCNA) was used to identify gene nodes and their interactions in (A) males, (B) females, and (C) combined-sexes datasets. For each analysis, within the module containing the Pink1 gene, Pink1 itself was the central node with the highest number of interactions with other genes including solute carriers, glutamate metabotropic receptors, and genes associated with protein localization. Strong connections between Pink1 and Krt2 and Hfe were found in both males and female datasets. In females a number of modules were significantly correlated with vocalization traits.
Overall, this work supports the premise that gene expression changes in the PAG may contribute to the vocal deficits observed in this PD rat model. Additionally, this dataset identifies genes that represent new therapeutic targets for PD voice disorders.
Parkinson’s disease (PD) is the second most common degenerative disorder affecting nearly 10 million people worldwide . The hallmark pathology of PD is death of dopaminergic neurons in the substantia nigra; however, pathology outside of dopamine loss often precedes clinical presentation of limb motor symptoms [2,3,4]. Vocal communication deficits, including hypokinetic dysarthria [5,6,7,8], are common and during the course of disease progression, over 90% of individuals present with these deficits. These signs negatively influence overall health, social interactions, employment, and quality of life [6, 7, 9, 10]. Despite this considerable clinical issue, the underlying central nervous system pathology that contributes to vocalization deficits in PD is poorly understood and understudied.
The standard pharmacological and surgical treatments aimed at restoring dopamine loss are generally ineffective for voice dysfunction which further suggests that vocal deficits are linked to other, unknown CNS pathologies [11,12,13]. Research suggests that PD in the central nervous system progresses in a caudal-to-rostral direction, with the brainstem impacted earlier than cortical regions [3, 4]. The periaqueductal gray (PAG) is a phylogenetically conserved brainstem nucleus implicated in the modulation and production of mammalian vocalizations [14,15,16]. Specifically, descending projections from the PAG modulate the motoneurons in the brainstem that control vocal fold adduction and respiration coordination [17, 18]. Additionally, the PAG is interconnected with the social behavior network  and the ascending cortical and limbic projections from the PAG are hypothesized to influence the motivation or emotional intent of species-specific vocalizations [19,20,21]. In humans, functional imaging studies demonstrate that the brainstem PAG is also involved in the circuitry and control of speech , and appears to be involved in the speech patterns of individuals with PD (where PAG connectivity is correlated to speech loudness) .
Genetic rodent models of PD are useful because they mimic aspects of idiopathic PD and translate with an early and progressive disease manifestation. The Pink1−/− rat, related to the PARK6 phenotype of human familial PD, demonstrates early motor deficits including changes to ultrasonic vocal production [24, 25]. Specifically, data show that male and female Pink1−/− rats exhibit early (2 months of age) ultrasonic vocalization changes to intensity (loudness) compared to non-affected wildtype (WT) control rats [24, 26]. Additional acoustic deficits include changes to frequency range in male rats at 8 months of age . Moreover, WT female conspecifics show decreased motivation to approach male Pink1−/− vocalization stimuli, suggesting that these deficits impair the social communication function . More recent data suggest that these vocal deficits are not rescued by pharmacological dopamine replacement (levodopa) , but can be, as in humans, modulated by a vocal-exercise therapy . In general, these deficits recapitulate findings observed in humans with PD, irregular vocal qualities that impact the ability to effectively communicate. In rats, the stimulation of the PAG increases production of vocalizations . Further, lesions to the PAG result in muteness , and its neurotransmitters and projections are connected to the peripheral structures (vocal fold) of vocalization . At 8 months, Pink1−/− male rats exhibit significant aggregation of insoluble α-synuclein, a pathological marker of PD . Moreover, PAG protein aggregation has also been implicated in vocal dysfunction in a mouse α-synuclein over expression model . Together, these data contribute to the working hypothesis that pathological changes in the PAG may account for aspects vocal dysfunction observed in the Pink1−/− PD rat.
Large scale gene expression analysis in post-mortem PD tissue and normal controls has led to the generation of detailed databases of the gene expression underlying the disease (see metanalysis ) and has contributed to the knowledge of genetic variants underlying the disease . Many of these studies are focused on the midbrain substantia nigra [34, 35], frontal lobe , and more recently blood-brain  and gut biomarkers . A primary advantage of a genetic model is the homogeneity of the sample that cannot be assessed with human post-mortem tissue. Because the etiology of PD is diverse and the pathology of vocal deficits is significantly understudied, in this study we used RNA sequencing to characterize differentially expressed genes within the PAG in the Pink1−/− rat compared to age-matched wildtype controls at 8 months of age. We tested the specific hypotheses:  loss of Pink1 in the PAG will influence expression of genes that interact with Pink1 ; loss of Pink1 will emphasize other genes that relate to PD; and  behavioral and bioinformatic approaches to data analysis will identify molecular targets important for rat vocalization. To further validate the Pink1 rat model, we compared our dataset to existing human databases to identify possible therapeutic targets for PD voice deficits.
Differential gene expression and KEGG pathway analysis between genotypes and sex
We first evaluated differential gene expression in the PAG between Pink1−/− and WT rats; Pink1 was the most significant downregulated gene in both males and females. Using the p value as a cutoff threshold (p < 0.05), there were 675 genes (520 annotated) within males and 1155 genes (973 annotated) within females found to be differentially expressed. Comparison of the top 500 downregulated genes showed a significant overlap between sexes (76 genes, hypergeometric p value< 0.0001). The gene overlap between the 500 upregulated gene list was also statistically significant (63 genes, hypergeometric p value < 0.0001).
Sex-specific enrichment analysis was used to compare the gene sets to existing KEGG network datasets (Fig. 1). Briefly, in males, pathways (genes) of interest include pentose and glucuronate interconversions (Akr1b10, Ugt1a9), and glycine, serine and threonine metabolism (Psat1, MaoB) (Fig. 1A). In females, glutathione metabolism (Gstm3, Sms), PPAR signaling (Ubc, Slc27a6), and metabolism (Tpmt, Gstm3) pathways were identified (Fig. 1B). Overlapping pathways include drug metabolism and metabolism of xenobiotics by cytochrome P450. There were more highlighted KEGG pathways in males  compared to females .
Enrichment analysis (using Enrichr) of the downregulated genes in females identified significant matches to genes downregulated in the CNS with PD in humans as well as PARK2 and PARK7 knockdown in human neuronal or tumor cell lines, suggesting the Pink1−/− rat model recapitulates genetic and idiopathic PD in humans (Table 1). Likewise, upregulated genes in Pink1−/−females match upregulated genes in the CNS in humans with PD. Similar matching of the Pink1−/− male brain gene expression with that found in humans with PD was identified. A weighted graph was created using the identified protein-protein interactions (Fig. 2) obtained from the STRING database. Initially, a human-rat list of 160 genes was entered into STRING and the default cutoff of gene-gene connectivity of 0.4 was used to exclude gene with little or no connectivity. The remaining genes were ranked based on the number of interactions with other unique genes and the top 83 genes with high interactions are plotted within the figure.
Downregulated genes of interest in females include genes involved in RNA binding and transcription (Hsp90aa1; Pih1d1; Hspe1), metabolism/mitochondria functions (Hsd17b; Adipor2; Cox5a; Cox6b1; Ugp2, Uqcrfs1), iron (Tf; Hspe1), vesicle signaling (Tspan8; Kif5c), ubiquitin (Rps27a; Psmd14; Hspd1; Dnaja1; BIrc2), ribosomal proteins (Rpl30, Rpl21, Rps24), and protein activity (Tomm20, Tomm7, Efr3a). In males, downregulated genes of interest that matched human Parkinson disease sequencing datasets included ubiquitin (Gpr37, Trim9), axon guidance (Sema3c), and syntaxin (Stxbp1).Notable upregulated genes that overlapped between the Pink1−/− female rat and human datasets included genes involved in neuron signaling (vesicle (Cd59) sodium channel activity/membrane depolarization (Scn1b), calcium signaling (Ednrb; Mylk), synapse organization (Lrp4)). There were associations with amyloid-beta binding (Clstn1, Cryab, Atp1a3) a neurological feature of protein aggregation in Alzheimer’s disease.
Gene interactions with the Pink1 gene using WGCNA
A primary goal of this work was to determine which genes interacted with Pink1. The annotated gene lists were run with WGCNA and data is presented in three ways:  male (Supplementary Table 1),  female (Supplementary Table 2), and a  combined-sexes approach (Supplementary Table 3). A number of modules (male = 29 modules; female = 26, combined sexes = 9) were generated; these tables can be sorted by module and trait (Supplementary Tables 4 and 5).
The following WGCNA modules included the Pink1 gene: for males the red module, females the lightcyan module, and for the combined-sexes the pink module. The sex-specific analysis of red and lightcyan modules that included Pink1 showed interactions with a different subset of genes, suggesting a possible sex-specific difference. For all three modules, Pink1, Krt2, Tert were the only genes that overlapped. A list of genes that appeared in at least two of the modules can be found in Table 2.
To determine the genes and their functions that interact with Pink1, the 45 genes that were in combined sexes pink module (displayed in Fig. 3) were put into the gene enrichment analysis tool to evaluate this gene list against preexisting data sets. Using Enrichr to evaluate gene ontology and biological processes, several areas of enrichment were identified including establishment of protein localization (Tert, Pih1d1) and glutamate receptor signaling (Grm2, Grm6).
Female ultrasonic vocalization behavior and WGCNA
Previously recorded female vocalization data  from WT and Pink1−/− rats was used in a WGCNA to identify vocalization-associated modules (Fig. 4). Using the 26 modules for female rats, the statistical links with acoustic and non-acoustic variables were identified (Table 3). Vocal modules with the top significant trait relationships included cyan, pink, midnightblue, darkgreen and grey60.
Because the cyan module was significantly linked with all vocal acoustic parameters of the ultrasonic calls, including duration, bandwidth, intensity, and peak frequency, as well as the percent of complex calls we focused enrichment analysis on this specific module (all other modules are sortable within Supplementary Tables 3 and 4). Performing gene enrichment analysis on the gene list from the cyan module, the top biological pathways were the regulation of membrane depolarization (GO:0003254) and indolalkylamine metabolic process (GO:0006586). Other pathways of interest include monoamine transport (GO:0015844); generation of neurons (GO:0048699); dopamine transmembrane transporter activity (GO:0005329); nervous system development (GO:0007399), axon guidance (GO:0007411); mRNA splicing (GO:00003898); and neuropeptide hormone activity (GO:0005184). Figure 5 represents top pathways (genes) of interest within the cyan module. Further evaluation of the Nts node suggests several biologically relevant connections (Fig. 6) including Slc6a4, Slc10a4, Ndnf, Pax5, Pax8, Sema3a, Gja5, and Gpr160.
In general, the understanding of Mendelian inherited forms of PD is limited, but necessary in order to provide insight into the genetic nature of the disease. At 8 months of age, the Pink1−/− genetic rat model of PD exhibits observable limb and cranial sensorimotor dysfunction including deficits in ultrasonic vocal communication that co-occur with aggregated α-synuclein in the PAG. The aim of the present study was to identify differentially expressed genes using high throughput RNA-sequencing in the brainstems of Pink1−/− male and female rats as compared to their respective age-matched wildtype (WT) controls and use bioinformatic approaches to further validate this PD model and highlight specific gene targets that may modulate vocal behaviors in the rat. This study successfully generated datasets of genes including lists of differentially expressed genes in both males and females; there was statistically significant overlap between both sexes. The upregulated and downregulated genes showed similarity to human Parkinson disease studies as well as PARK gene models which provides further validity at the gene-level. WGCNA bioinformatic approaches highlighted several gene pathways that may be important for vocalization. Together these findings highlight new directions for targeting vocal biology pathways.
Phenotypically, motor signs of PD manifest differently in males compared to female Pink1−/− rats . For example, females do not show the classical limb motor deficits in tapered balance beam or cylinder tests compared to progressive slowness and increased number of errors in males. Both females and males show decreases in vocal loudness; however, males also have changes to other vocal variables including reductions in peak frequency and bandwidth. In this study, the comparison of sequencing runs in males and females yielded strong significant relationships between sexes suggesting that the removal of the Pink1 gene in both sexes yields similar genetic expression in the PAG. What contributes to the observed differences in behavioral phenotypes is unknown and should be addressed in future work. In normal cases, Pink1 has a protective role against mitochondrial quality, dysfunction, and mitophagy by activating a mitochondrial damage-response signaling pathway. The Pink1: Parkin: Ubiquitin interaction has been well studied . Here, we have recapitulated the important relationship between Pink1 and Ubc in both male and female gene datasets as Pink1 and Ubc were the top genes within the same-sex module pink. Activation of this pathway results in selective autophagy. Thus, loss of function mutations in this process ultimately lead to increases in oxidative buildup and mitochondrial damage which has been previously reported [42, 45,46,47,48].
Interestingly, comparison of the modules that include Pink1 yielded different co-expressed genes; the only overlapping genes were Pink1, Kert2, and Hfe. The interaction between Pink1 and Hfe suggests a mitochondrial: iron link, that has not previously been investigated in this particular in vivo model. Hfe expression modulates cellular iron absorption (reviewed in ), with a strong correlation to PD. For example, abnormal iron concentrations in the basal ganglia are hypothesized to induce PD symptoms . Previous reports in the Pink1−/− rat suggest the presence of mitophagy including striatal mitochondrial proteomic alterations, and challenges within the mitochondria respiratory system. Iron accumulation within the neuron may be a mechanism of neurodegeneration and additional changes with age (i.e. disease progression). Pathway analysis from the female data showed significant genotype contrasts for the glutathione pathway  which is consistent with the strong links between iron accumulation, mitochondrial dysfunction (removal of impaired mitochondria), and resulting oxidative inflammation in the central and peripheral nervous system . Hfe polymorphisms also affect cellular glutamate levels in cell lines ; glutamate receptor alterations observed in this study are discussed below. Our data suggest the Pink1−/− rat is a representative model of these types of complex interactions.
Activation of the PAG has significant effects on neurotransmitter release (reviewed in ), including GABAergic and glutamatergic signaling which suggests a role for multiple neurotransmitter systems in vocal biology [55,56,57]. Previously, we have shown that glutamate decarboxylase 1 (Gad1) in the PAG is significantly reduced in Pink1−/− rats  and our most recent work suggests that modulation of the GABA neurotransmitter represents therapeutic targets for rescuing vocalization. For example, within the VTA of Pink1−/− rats Gad gene expression is upregulated with an intensive behavioral vocalization intervention . In this dataset, there are several glutamatergic pathways (genes) that were identified with gene enrichment analysis from the combined sex pink module. Two genes, Grm2 (codes for L-glutamate) and Grm6 (glutamate metabotropic receptor 6) are of particular interest as metabotropic glutamate receptors have been hypothesized to be promising anti-Parkinsonian drugs . And, glutamate remains an important excitatory neurotransmitter for the production of vocalizations. For example, in the squirrel monkey, glutamate induces vocalizations when injected into the periaqueductal gray , and blocking glutamate receptors suppresses vocalization (as well as locomotion) in mice . Thus, compounds that promote glutamate receptor expression and/or upregulate glutamate production should be analyzed in future studies on Pink1−/− vocal behaviors and as vocal modulators.
In addition to GABA/glutamatergic systems, the gene enrichment analysis of the female vocalization module cyan (this module was significantly correlated to all USV parameters) suggest an important relationship between monoamine solute carriers and neuromodulators in female rat vocalization. From the enrichment analysis (visualized in Fig. 5), the vesicular monoamine transporter Slc6a4 (transports serotonin) and Ddc (dopamine/serotonin) both directly interact with neurotensin (Nts). Nts is a neuromodulator that has been linked to social behavior and song production in birds, its mRNA expression is correlated to song production . In rats, neurotensin agonists decreases stress induced 22-kHz vocalizations , but its effects on 50-kHz social vocalizations is unknown. Moreover, neurotensin shows significant relationships with other genes; for example, there are interactions with PAX transcription factor genes and neuron derived neurotrophic factor. These gene targets provide insight into the PAG-specific relationships with vocalization parameters and are promising therapeutic targets for observed decrease in vocal loudness in the Pink1−/− female rat model and may be translatable to the male acoustic dysfunction. In general, this work presents researchers with data mining opportunities, specifically to identify relationships between acoustic traits and any subset of these genes.
The data suggest that loss of Pink1 influences gene pathways within the brainstem PAG. Additionally, this study provides links and validation of this rat model as being useful for studying the disease as our data sets also match to human idiopathic Parkinson and PARK2 gene databases. These results indicate a key role for specific genes in vocal behavior and suggest potential drug targets for PD vocal deficits.
Rats, housing, and acclimation
A total of 16 rats were used in this study; 8 male (n = 4 Pink1−/−, n = 4 WT) and 8 female (n = 4 Pink1−/−, n = 4 WT) Long Evans rats (SAGE™ Research Labs, Boyertown, PA, USA) were aged to 8 months [25, 64]. Rats were housed in groups of two (within like genotypes/sex) in standard polycarbonate cages (290 mm × 533 mm × 210 mm) with corn cob bedding on a reversed 12:12 h light: dark cycle. All test procedures occurred during the dark period of the cycle. Food and water were available ad libitum. All procedures were approved by the University of Wisconsin-Madison Animal Care and Use Committee and were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory animals .
Female vocalization testing and analysis
A study of female-specific behavior was initially published in Marquis et al., 2019 . Methods for data collection have been previously described [24, 58, 66]. Frequency-modulated acoustic (duration (s), bandwidth (Hz), intensity (loudness, dB), and peak frequency (Hz)) and non-acoustic variables (number of calls, percent of complex calls, call rate (calls/90 s)) were used in the statistical correlation analysis, discussed below.
Tissue harvest and processing
Rats were deeply anesthetized with isoflurane and rapidly decapitated. The brains were dissected and immediately frozen and stored at − 80 °C. Brains were sliced coronally on a freezing cryostat at 250 μm thickness at − 15 °C and mounted on glass slides. Anatomically equivalent sections were used from each rat and 2 mm tissue punches (left and right) were collected within the PAG (Bregma ~ − 7.8 mm) using the Brain Punch Set (FST 18035–02, Foster City, CA, USA) as in previously reported methods . Tissue samples were transferred to microcentrifuge tubes and stored at − 80 °C until processing for RNA.
Sample order was randomized throughout the molecular portion of the study. Brain tissue was homogenized with an electric sonic dismembrator (Fisher Scientific, Hampton, NH, USA) and was extracted using the Bio-Rad Aurum Total RNA Fatty and Fibrous Tissue Kit (Catalog No. 732–6830; Bio-Rad, Hercules, CA, USA). The total RNA was measured using a Nanodrop system (Thermo Scientific, Wilmington, DE, USA) and yielded significant concentrations of RNA. Additionally, the 28S:18S rRNA was quantified with an Agilent RNA 6000 Pico kit (Eukaryote Total RNA Pico, Agilent Technologies, Santa Clara, CA) and verified using a Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). All samples had an A260/A280 ratio that fell within the 1.92–2.02 range and showed satisfactory marker and ribosomal peaks. RNA Integrity Numbers were above 7.2.
Library construction and RNA sequencing
All RNA-sequencing procedures followed guidelines by ENCODE and were performed by the University of Wisconsin-Madison Biotechnology Center’s Next Generation Sequencing Facility. The Illumina® Total RNA-Seq TruSeq platform (Illumina Inc., San Diego, CA, USA) was used to profile differential expression of genes in the PAG between Pink1−/− and WT rats. Males (WT and Pink1−/−) and females (WT and Pink1−/−) were processed in separate batches, but preparation and processing of tissue as described above was identical. Briefly, the Stranded Total RNA Library Prep Kit was used to remove rRNA, and a sequencing library was generated. To prepare libraries, 500 ng was used as an input, rRNA reduction was done using H/M/R reagents (TruSeq Stranded Total RNA kit), samples were processed with RNA clean beads and 70% ethanol. The 3′ and 5′ adapters were ligated to small RNAs, followed by reverse transcription to obtain single-stranded complementary DNA. PCR-amplification used a universal primer and a primer containing a unique index sequence (Illumina, Inc). The amplified complementary DNA libraries were gel purified and used to construct RNA sequencing libraries. Libraries were quantified using Qubit DNA HS kit, diluted 1:100, assayed on Agilent DNA1000 chip. Libraries showed no adapter dimer contamination. Sequencing was performed on an HiSeq 2000 high-throughput sequencing system within a single run (Illumina, Inc). Adaptor sequences, contamination and low-quality reads were removed. Reads were mapped to the annotated rat (Rattus norvegicus) genome in Ensembl . Technical quality was determined using several parameters. Briefly, trimming software skewer  was used to preprocess raw fastq files. Combined cycle base quality, per cycle base frequencies and average base quality, relative 3-Kmer diversity, Phred quality distribution, mean quality distribution, average read length and read occurrence distribution were used as additional quality control measures on the read data. Additionally, the biological coefficient of variation was estimated to be approximately 0.12. As such, a number of differentially expressed genes were identified, including Pink1−/− (Raw data (RSEM), is displayed in Supplementary Tables 6 and 7).
Differential gene expression analysis and KEGG pathway enrichment
Analysis was performed with the glm using the EdgeR Bioconductor Package, v. 3.9 . The p-value cutoff was set to 0.05 for significance. A Benjamini-Hochberg correction was applied to control the False Discovery Rate (FDR) . Statistics and FDR for both male and female datasets (EdgeR Results) are provided in Supplementary Tables 8 and 9, respectively. RSEM approach for normalizing RNA seq data was used ; raw data were uploaded to the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150939; GSE150939). Genes were ranked according to p value and sorted by up- or down-regulation. The average LogCPM for all genes was similar for males (3.69) and females (3.92). The top 500 genes for males and females (both up- and down -regulated; see Supplementary Table 10) were compared to identify similar expression changes between the sexes in response to knockout (hypergeometric test). The cutoff of 500 was used as this is expected to capture significant biological changes in both directions and this cutoff is the same as used by GEO2Enrichr function to extract information from gene expression studies .
KEGG pathway enrichment was used to identify biological themes in the collection of differentially expressed genes. To conduct this test, all genes with a p-value < 0.05 were selected and for each pathway input genes that are part of the specific pathway were counted. The list of genes in every pathway was tested for over- or under-representation with respect to the input list of DE genes. The number of “background” genes was determined by counting the number of protein coding genes [Biotype] in the annotation model (described in .
Weighted gene co-expression network (WGCNA) analysis
WGCNA was used to construct gene co-expression networks and gene modules from the gene expression datasets. Briefly, data were log 2 + 1 transformed, low expression genes were removed (specifically, using the filter function in EdgeR genes that had no expression in any of the individuals were removed) and WGCNA was run (on males: 13253 genes; females: 13277 genes; combined-sexes: 13253 genes) using R software (https://www.r-project.org/) . Using a weighted network of genes and expression correlates (nodes and edges), correlations were raised to a soft thresholding power β of 12. Searchable networks were created for:  male rats,  female rats, and  combined-sexes (Supplementary Tables 1, 2, 3, 4 and 5).
Unsupervised hierarchical clustering for WGCNA included: the minimum module size of 30 genes, the signed mode, the deepSplit parameter set to 2, the mergeCutHeight parameter set to 0.15, and a threshold setting for merging modules of 0.25. Module eigengene values were also evaluated in terms of their genotype and vocal traits. Modules were visualized using Cytoscape (v3.7.1 48; https://cytoscape.org/) and the network file was exported and manually trimmed as to consist of genes of interest and the specific gene-to-gene correlations.
Gene enrichment analysis
In order to identify genes that are over-represented in the data and associated with particular functions and relevance to PD, gene enrichment analysis was performed on both the differentially expressed gene dataset and the gene modules produced by the WGCNA analysis. Gene enrichment was performed with EnrichR .
Genes dysregulated in Pink1−/− rats were compared with genes dysregulated in the same direction in human PD using multiple datasets (GSE7621, GSE19587, GSE17204) and the meta-analysis list in . This human-rat list of 160 genes was entered into STRINGdp (v 2.0; Search Tool for the Retrieval of Interacting Genes/Proteins, http://string.embl.de/) to identify genes with high protein-protein interactions (Supplementary Table 11).
Availability of data and materials
The datasets generated during the current study are available in the Gene Expression Omnibus repository (GSE150939). Website: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150939
The datasets analyzed during the current study are available in the Gene Expression Omnibus repositories (GSE7621, GSE19587, GSE17204) and as described in Table 1 as well as the annotated rat genome (Rat Rnor_6.0 assembly; https://www.ebi.ac.uk/ena/browser/view/GCA_000001895.4).
de Lau LML, Giesbergen PCLM, de Rijk MC, Hofman A, Koudstaal PJ, Breteler MMB. Incidence of parkinsonism and Parkinson disease in a general population: the Rotterdam study. Neurology. 2004;63(7):1240–4.
Braak H, Ghebremedhin E, Rüb U, Bratzke H, Tredici K. Stages in the development of Parkinson’s disease-related pathology. Cell Tissue Res. 2004;318(1):121–34.
Braak H, Tredici KD, Rüb U, de Vos RAI, Jansen Steur ENH, Braak E. Staging of brain pathology related to sporadic Parkinson’s disease. Neurobiol Aging. 2003;24(2):197–211.
Hawkes CH, Del Tredici K, Braak H. A timeline for Parkinson's disease. Parkinsonism Relat Disord. 2010;16(2):79–84.
Harel BT, Cannizzaro MS, Cohen H, Reilly N, Snyder PJ. Acoustic characteristics of Parkinsonian speech: a potential biomarker of early disease progression and treatment. J Neurolinguist. 2004;17(6):439–53.
Hartelius L, Svensson P. Speech and swallowing symptoms associated with Parkinson’s disease and multiple sclerosis: a survey. Folia Phoniatrica et Logopaedica. 1994;46(1):9–17.
Ho AK, Iansek R, Marigliani C, Bradshaw JL, Gates S. Speech impairment in a large sample of patients with Parkinson's disease. Behav Neurol. 1999;11(3):131–7.
Holmes RJ, Oates JM, Phyland DJ, Hughes AJ. Voice characteristics in the progression of Parkinson's disease. International journal of language & communication disorders / Royal College of Speech & Language Therapists. 2000;35(3):407–18.
Plowman-Prine EK, Sapienza CM, Okun MS, Pollock SL, Jacobson C, Wu SS, et al. The relationship between quality of life and swallowing in Parkinson's disease. Mov Disord. 2009;24(9):1352–8.
Marras C, McDermott MP, Rochon PA, Tanner CM, Naglie G, Lang AE. Predictors of deterioration in health-related quality of life in Parkinson's disease: results from the DATATOP trial. Mov Disord. 2008;23(5):653–9.
Ciucci MR, Grant LM, Rajamanickam ESP, Hilby BL, Blue KV, Jones CA, et al. Early identification and treatment of communication and swallowing deficits in Parkinson disease. Semin Speech Lang. 2013;34(03):185–202.
Plowman-Prine EK, Okun MS, Sapienza CM, Shrivastav R, Fernandez HH, Foote KD, et al. Perceptual characteristics of Parkinsonian speech: a comparison of the pharmacological effects of levodopa across speech and non-speech motor systems. NeuroRehabilitation. 2009;24(2):131–44.
Sanabria J, Ruiz PG, Gutierrez R, Marquez F, Escobar P, Gentil M, et al. The effect of levodopa on vocal function in Parkinson's disease. Clin Neuropharmacol. 2001;24(2):99–102.
Larson CR. The midbrain periaqueductal gray: a brainstem structure involved in vocalization. J Speech Lang Hear Res. 1985;28(2):241–9.
Larson CR, Kistler MK. The relationship of periaqueductal gray neurons to vocalization and laryngeal EMG in the behaving monkey. Exp Brain Res. 1986;63(3):596–606.
Goodson JL. The vertebrate social behavior network: evolutionary themes and variations. Horm Behav. 2005;48(1):11–22.
Zhang SP, Davis PJ, Bandler R, Carrive P. Brain stem integration of vocalization: role of the midbrain periaqueductal gray. J Neurophysiol. 1994;72(3):1337–56.
Davis PJ, Zhang SP, Winkworth A, Bandler R. Neural control of vocalization: respiratory and emotional influences. J Voice. 1996;10(1):23–38.
Larson CR. Activity of PAG neurons during conditioned vocalization in the macaque monkey. In: Depaulis A, Bandler R, editors. The midbrain periaqueductal gray matter: functional, anatomical, and neurochemical organization. Boston, MA: Springer US; 1991. p. 23–40.
Shiba K, Satoh I, Kobayashi N, Hayashi F. Multifunctional laryngeal Motoneurons: an intracellular study in the cat. J Neurosci. 1999;19(7):2717–27.
Fardin V, Oliveras J-L, Besson J-M. A reinvestigation of the analgesic effects induced by stimulation of the periaqueductal gray matter in the rat. I. the production of behavioral side effects together with analgesia. Brain Res. 1984;306(1):105–23.
Schulz GM, Varga M, Jeffires K, Ludlow CL, Braun AR. Functional Neuroanatomy of human vocalization: an H215O PET study. Cereb Cortex. 2005;15(12):1835–47.
Rektorova I, Mikl M, Barrett J, Marecek R, Rektor I, Paus T. Functional neuroanatomy of vocalization in patients with Parkinson's disease. J Neurol Sci. 2012;313(1):7–12.
Grant LM, Kelm-Nelson CK, Hilby BL, Blue KV, Rajamanickam ESP, Pultorak J, et al. Evidence for early and progressive ultrasonic vocalization and oromotor deficits in a PINK1 knockout rat model of Parkinson disease. J Neurosci Res. 2015;93(11):1713–27.
Dave KD, De Silva S, Sheth NP, Ramboz S, Beck MJ, Quang C, et al. Phenotypic characterization of recessive gene knockout rat models of Parkinson's disease. Neurobiol Dis 2014;70(0):190–203.
Marquis JM, Lettenberger SE, Kelm-Nelson CA. Early-onset Parkinsonian behaviors in female Pink1−/− rats. Behavioral Brain Research. 2020;13(377):112175.
Pultorak J, Kelm-Nelson CK, Holt LR, Blue KV, Ciucci MR, Johnson AM. Decreased approach behavior and nucleus accumbens immediate early gene expression in response to Parkinsonian ultrasonic vocalizations in rats. Soc Neurosci. 2015;11(4):365–79.
Kelm-Nelson CAaC, M.R. Levodopa improves a subset of motor function associated with nigrostriatal deficits in a Pink1 −/− rat model of Parkinson disease. In publication.
Kelm-Nelson CA, Yang KM, Ciucci MR. Exercise effects on early vocal ultrasonic communication dysfunction in a PINK1 knockout model of Parkinson's disease. J Parkinsons Dis. 2015;5(4):749–63.
Yajima Y, Hayashi Y, Yoshi N. The midbrain central gray substance as a highly sensitive neural structure for the production of ultrasonic vocalization in the rat. Brain Res. 1980;198(2):446–52.
Kelm-Nelson CA, Stevenson SA, Ciucci MR. Atp13a2 expression in the periaqueductal gray is decreased in the Pink1 −/− rat model of Parkinson disease. Neurosci Lett. 2016;621:75–82.
Grant LM, Richter FR, Miller JE, White SA, Fox CM, Zhu C, et al. Vocalization deficits in mice over-expressing alpha-synuclein, a model of pre-manifest Parkinson's disease. Behav Neurosci. 2014;128(2):110–21.
Kelly J, Moyeed R, Carroll C, Albani D, Li X. Gene expression meta-analysis of Parkinson’s disease and its relationship with Alzheimer’s disease. Molecular Brain. 2019;12(1):16.
Lesnick TG, Papapetropoulos S, Mash DC, Ffrench-Mullen J, Shehadeh L, de Andrade M, et al. A genomic pathway approach to a complex disease: axon guidance and Parkinson disease. PLoS Genet. 2007;3(6):e98.
Schulze M, Sommer A, Plotz S, Farrell M, Winner B, Grosch J, et al. Sporadic Parkinson's disease derived neuronal cells show disease-specific mRNA and small RNA signatures with abundant deregulation of piRNAs. Acta neuropathologica communications. 2018;6(1):58.
Dumitriu A, Golji J, Labadorf AT, Gao B, Beach TG, Myers RH, et al. Integrative analyses of proteomics and RNA transcriptomics implicate mitochondrial processes, protein folding pathways and GWAS loci in Parkinson disease. BMC Med Genet. 2016;9(1):5.
Chatterjee P, Roy D. Comparative analysis of RNA-Seq data from brain and blood samples of Parkinson's disease. Biochem Biophys Res Commun. 2017;484(3):557–64.
Heintz-Buschart A, Pandey U, Wicke T, Sixel-Döring F, Janzen A, Sittig-Wiegand E, et al. The nasal and gut microbiome in Parkinson's disease and idiopathic rapid eye movement sleep behavior disorder. MovDisord. 2018;33(1):88–98.
Clements CM, McNally RS, Conti BJ, Mak TW, Ting JPY. DJ-1, a cancer- and Parkinson's disease-associated protein, stabilizes the antioxidant transcriptional master regulator Nrf2. Proc Natl Acad Sci U S A. 2006;103(41):15091–6.
Lewandowski NM, Ju S, Verbitsky M, Ross B, Geddie ML, Rockenstein E, et al. Polyamine pathway contributes to the pathogenesis of Parkinson disease. Proc Natl Acad Sci U S A. 2010;107(39):16970–5.
Gong Y, Zack TI, Morris LG, Lin K, Hukkelhoven E, Raheja R, et al. Pan-cancer genetic analysis identifies PARK2 as a master regulator of G1/S cyclins. Nat Genet. 2014;46(6):588–94.
Gautier CA, Kitada T, Shen J. Loss of PINK1 causes mitochondrial functional defects and increased sensitivity to oxidative stress. Proc Natl Acad Sci U S A. 2008;105(32):11364–9.
Foti R, Zucchelli S, Biagioli M, Roncaglia P, Vilotti S, Calligaris R, et al. Parkinson disease-associated DJ-1 is required for the expression of the glial cell line-derived neurotrophic factor receptor RET in human neuroblastoma cells. J Biol Chem. 2010;285(24):18565–74.
Koyano F, Okatsu K, Kosako H, Tamura Y, Go E, Kimura M, et al. Ubiquitin is phosphorylated by PINK1 to activate parkin. Nature. 2014;510(7503):162–6.
Heeman B, Van den Haute C, Aelvoet S-A, Valsecchi F, Rodenburg RJ, Reumers V, et al. Depletion of PINK1 affects mitochondrial metabolism, calcium homeostasis and energy maintenance. J Cell Sci. 2011;124(7):1115–25.
Liu W, Vives-Bauza C, Acín-Peréz R, Yamamoto A, Tan Y, Li Y, et al. PINK1 defect causes mitochondrial dysfunction, proteasomal deficit and α-Synuclein aggregation in cell culture models of Parkinson's disease. PLoS One. 2009;4(2):e4597.
Piccoli C, Sardanelli A, Scrima R, Ripoli M, Quarato G, D'Aprile A, et al. Mitochondrial respiratory dysfunction in familiar parkinsonism associated with PINK1 mutation. Neurochem Res. 2008;33(12):2565–74.
Ferris CF, Morrison TR, Iriah S, Malmberg S, Kulkarni P, Hartner JC, et al. Evidence of neurobiological changes in the Presymptomatic PINK1 knockout rat. J Parkinsons Dis. 2018;8(2):281–301.
Eisenstein RS. Interaction of the hemochromatosis gene product Hfe with transferrin receptor modulates cellular Iron metabolism. Nutr Rev. 1998;56(12):356–8.
Costello DJ, Walsh SL, Harrington HJ, Walsh CH. Concurrent hereditary haemochromatosis and idiopathic Parkinson's disease: a case report series. J Neurol Neurosurg Psychiatry. 2004;75(4):631–3.
Liddell JR, White AR. Nexus between mitochondrial function, iron, copper and glutathione in Parkinson's disease. Neurochem Int. 2018;117:126–38.
Urrutia PJ, Mena NP, Nunez MT. The interplay between iron accumulation, mitochondrial dysfunction, and inflammation during the execution step of neurodegenerative disorders. Front Pharmacol. 2014;5:38.
Mitchell RM, Lee SY, Simmons Z, Connor JR. HFE polymorphisms affect cellular glutamate regulation. Neurobiol Aging. 2011;32(6):1114–23.
Jurgens U. The role of the periaqueductal grey in vocal behaviour. Behav Brain Res. 1994;62(2):107–17.
Jurgens U, Lu CL. Interactions between glutamate, GABA, acetylcholine and histamine in the periaqueductal gray's control of vocalization in the squirrel monkey. Neurosci Lett. 1993;152(1–2):5–8.
Jürgens U, Lu CL. The effects of Periaqueductally injected transmitter antagonists on forebrain-elicited vocalization in the squirrel monkey. Eur J Neurosci. 1993;5(6):735–41.
Jurgens U, Richter K. Glutamate-induced vocalization in the squirrel monkey. Brain Res. 1986;373(1–2):349–58.
Kelm-Nelson CA, Trevino MA, Ciucci MR. Quantitative analysis of Catecholamines in the Pink1 −/− rat model of early-onset Parkinson's disease. Neuroscience. 2018;379:126–41.
Stevenson SA, Ciucci MR, Kelm-Nelson CA. Intervention changes acoustic peak frequency and mesolimbic neurochemistry in the Pink1−/− rat model of Parkinson disease. PLoS One. 2019;14(8):e0220734.
Nickols HH, Conn PJ. Development of allosteric modulators of GPCRs for treatment of CNS disorders. Neurobiol Dis. 2014;61:55–71.
Roccaro-Waldmeyer DM, Girard F, Milani D, Vannoni E, Prétôt L, Wolfer DP, et al. Eliminating the VGlut2-dependent Glutamatergic transmission of Parvalbumin-expressing neurons leads to deficits in locomotion and vocalization, decreased pain sensitivity, and increased dominance. Front Behav Neurosci. 2018;12:146.
Merullo DP, Asogwa CN, Sanchez-Valpuesta M, Hayase S, Pattnaik BR, Wada K, et al. Neurotensin and neurotensin receptor 1 mRNA expression in song-control regions changes during development in male zebra finches. Dev Neurobiol. 2018;78(7):671–86.
Prus AJ, Hillhouse TM, LaCrosse AL. Acute, but not repeated, administration of the neurotensin NTS1 receptor agonist PD149163 decreases conditioned footshock-induced ultrasonic vocalizations in rats. Prog Neuro-Psychopharmacol Biol Psychiatry. 2014;49:78–84.
Baptista MAS, Dave KD, Sheth NP, De Silva SN, Carlson KM, Aziz YN, et al. A strategy for the generation, characterization and distribution of animal models by the Michael J. Fox Foundation for Parkinson’s research. Dis Model Mech. 2013;6(6):1316–24.
Guide for the Care and Use of Laboratory Animals. Washington DC: National Academy of Sciences.; 2011.
Kelm-Nelson CA, Brauer AFL, Ciucci MR. Vocal training, levodopa, and environment effects on ultrasonic vocalizations in a rat neurotoxin model of Parkinson disease. Behav Brain Res. 2016;307:54–64.
Aken BL, Ayling S, Barrell D, Clarke L, Curwen V, Fairley S, et al. The Ensembl gene annotation system. Database. 2016;2016.
Jiang H, Lei R, Ding S-W, Zhu S. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinformatics. 2014;15(1):182.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26.
Reiner A, Yekutieli D, Benjamini Y. Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003;19(3):368–75.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.
Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128.
Subhash S, Kanduri C. GeneSCF: a real-time based functional enrichment tool with support for multiple organisms. BMC Bioinformatics. 2016;17(1):365.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.
The authors thank the University of Wisconsin-Madison Biotechnology Center Gene Expression Center & DNA Sequencing Facility for providing library preparation and next generation sequencing services. We would like to acknowledge the UW Bioinformatics Resource Center for their assistance with data analysis.
This work was supported by the National Institutes of Health (R21 DC016135 (Kelm-Nelson)). Funding bodies played no role in the design of the study or analysis or interpretation of data or in writing the manuscript.
Ethics approval and consent to participate
All procedures were approved by the University of Wisconsin-Madison Animal Care and Use Committee (IACUC; protocol M005177) and were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory animals.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Gene modules identified by WGCNA in male rats.
Gene modules identified by WGCNA for female rats.
Gene modules identified by WGCNA for all rats.
Relationships between WGCNA identified modules and female vocal behaivor.
Relationships between WGCNA identified modules and female vocal behaivor.
Male RNA-Seq data used to analyze differential gene expression using EdgeR.
Female RNA-Seq data used to analyze differential gene expression using EdgeR.
Results of a differential gene expression in female rats.
Results of a differential gene expression in female rats.
Top 500 common up- and down-regulated genes in male and female Pink1−/− rats.
STRING Gene List.
About this article
Cite this article
Kelm-Nelson, C.A., Gammie, S. Gene expression within the periaqueductal gray is linked to vocal behavior and early-onset parkinsonism in Pink1 knockout rats. BMC Genomics 21, 625 (2020). https://doi.org/10.1186/s12864-020-07037-4
- Parkinson’s disease
- RNA sequencing
- Periaqueductal gray