Skip to main content
  • Research article
  • Open access
  • Published:

Gene expression within the periaqueductal gray is linked to vocal behavior and early-onset parkinsonism in Pink1 knockout rats



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 [1]. 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 [16] 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 [22], and appears to be involved in the speech patterns of individuals with PD (where PAG connectivity is correlated to speech loudness) [23].

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 [24]. Moreover, WT female conspecifics show decreased motivation to approach male Pink1−/− vocalization stimuli, suggesting that these deficits impair the social communication function [27]. More recent data suggest that these vocal deficits are not rescued by pharmacological dopamine replacement (levodopa) [28], but can be, as in humans, modulated by a vocal-exercise therapy [29]. 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 [15]. Further, lesions to the PAG result in muteness [17], and its neurotransmitters and projections are connected to the peripheral structures (vocal fold) of vocalization [30]. At 8 months, Pink1−/− male rats exhibit significant aggregation of insoluble α-synuclein, a pathological marker of PD [31]. Moreover, PAG protein aggregation has also been implicated in vocal dysfunction in a mouse α-synuclein over expression model [32]. 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 [33]) and has contributed to the knowledge of genetic variants underlying the disease [33]. Many of these studies are focused on the midbrain substantia nigra [34, 35], frontal lobe [36], and more recently blood-brain [37] and gut biomarkers [38]. 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: [1] loss of Pink1 in the PAG will influence expression of genes that interact with Pink1 [2]; loss of Pink1 will emphasize other genes that relate to PD; and [3] 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 [7] compared to females [4].

Fig. 1
figure 1

Ordered KEGG pathway rank in males and females. Significant metabolic pathways in (A) male and (B) female periaqueductal gray RNA-sequencing datasets (contrast Pink1−/−: WT). RNO = KEGG pathway entry for Rattus norvegicus. Significantly enriched KEGG pathways are ordered from most to least significant. The number of genes in the specified KEGG pathway are denoted by the size of a square (males) or circle (females) and the specific number of expressed genes for each sex, respectively. The color filling each circle corresponds to the specific KEGG pathway description

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.

Table 1 Example of directional overlap with human PD datasets. Grey text is associated with non-significant p values
Fig. 2
figure 2

Protein-protein interaction networks. STRING model of protein interactions via genes that are commonly up and down regulated in male and female Pink1−/− rats and in humans with PD. Nodes are proteins and edges indicate protein interactions

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: [1] male (Supplementary Table 1), [2] female (Supplementary Table 2), and a [3] 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.

Table 2 List of gene candidates appearing in modules with Pink1 in more than one condition (male, female, combined-sexes)

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).

Fig. 3
figure 3

Characterization of the combined-sexes MEPink WGCNA module. The interaction network of the differentially expressed genes was visualized using Cytoscape 3.3.0 software. Data were plotted using weight (level of significance) as a factor. Only annotated genes are plotted from the combined sexes dataset (37 nodes, 203 edges). Pink1 (yellow) was the top gene with the highest number of significant connections. Lines represent significant correlations between two genes. Red lines represent significant correlations between Pink1 and other module genes; black lines are connections between module genes (note: node colors do not correlate to the WGCNA module color)

Female ultrasonic vocalization behavior and WGCNA

Previously recorded female vocalization data [26] 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.

Fig. 4
figure 4

Module-trait relationships. Detected statistical associations between expression profiles for each of the WGCNA identified modules (y-axis) and female frequency modulated vocal behavior (x-axis). Boxes contain correlation coefficients (p values) where Red = strong positive correlation; Green = strong negative correlation). The genes comprising each module are listed in Supplementary Tables 3 and 4. *** indicate modules that correlated significantly with vocalization variables (see Table 3)

Table 3 Significant female vocalization modules

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.

Fig. 5
figure 5

Gene-vocalization links. MeCyan gene module is the most positively correlated module with all female rat acoustic properties. ENRICHR gene enrichment software was used to examine gene pathways within this module, as presented here. Data were plotted using weight as a factor using Cytoscape 3.3.0 software. The most significant connections are plotted using weight as a factor; lines represent significant correlations between genes. Colors are indicative of the same enriched pathway: Dark purple (monoamine transmembrane transporter activity; dopamine transmembrane transporter activity; indolalkylamine metabolic process): Slc6a4, Ddc, Tph2; Light yellow (nervous system development, axon guidance/axonogenesis): Gss, Gfra1, Sim1, Scg1; Brown (RNA splicing; mRNA processing; peptidyl-prolyl cis-trans isomerase activity): Ppil3, Fkbp14, Wbp11, Magohb; Light Purple (neuron migration, axon guidance; generation of neurons and depolarization): Sema3a, Ndnf, Runx2, Smad7, Gja5, Rangrf; Light Pink (neuropeptide, cell signaling; transcription): Pax8, Pax5, Nts

Fig. 6
figure 6

Neurotensin connections. The Nts node and its connections (from the MECyan module). Data were plotted using weight as a factor using Cytoscape 3.3.0 software. The most significant connections plotted using weight as a factor; red lines represent significant correlations between genes. Genes in yellow represent genes of interest, other connecting genes are in blue


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 [26]. 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 [44]. 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 [49]), with a strong correlation to PD. For example, abnormal iron concentrations in the basal ganglia are hypothesized to induce PD symptoms [50]. 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 [51] 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 [52]. Hfe polymorphisms also affect cellular glutamate levels in cell lines [53]; 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 [54]), 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 [58] 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 [59]. 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 [60]. 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 [57], and blocking glutamate receptors suppresses vocalization (as well as locomotion) in mice [61]. 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 [62]. In rats, neurotensin agonists decreases stress induced 22-kHz vocalizations [63], 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 [65].

Female vocalization testing and analysis

A study of female-specific behavior was initially published in Marquis et al., 2019 [26]. 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 [31]. Tissue samples were transferred to microcentrifuge tubes and stored at − 80 °C until processing for RNA.

RNA preparation

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 [67]. Technical quality was determined using several parameters. Briefly, trimming software skewer [68] 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 [69]. The p-value cutoff was set to 0.05 for significance. A Benjamini-Hochberg correction was applied to control the False Discovery Rate (FDR) [70]. 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 [71]; raw data were uploaded to the Gene Expression Omnibus (; 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 [72].

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 [73].

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 ( [74]. 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: [1] male rats, [2] female rats, and [3] 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; 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 [72].

Protein-protein interactions

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 [33]. This human-rat list of 160 genes was entered into STRINGdp (v 2.0; Search Tool for the Retrieval of Interacting Genes/Proteins, 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:

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;



Parkinson’s disease


periaqueductal gray












  1. 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.

    Article  PubMed  Google Scholar 

  2. 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.

    Article  PubMed  Google Scholar 

  3. 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.

    Article  PubMed  Google Scholar 

  4. Hawkes CH, Del Tredici K, Braak H. A timeline for Parkinson's disease. Parkinsonism Relat Disord. 2010;16(2):79–84.

    Article  PubMed  Google Scholar 

  5. 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.

    Article  Google Scholar 

  6. 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.

    Article  PubMed  CAS  Google Scholar 

  7. 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.

    Article  PubMed  CAS  Google Scholar 

  8. 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.

    Article  CAS  Google Scholar 

  9. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  10. 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.

    Article  PubMed  Google Scholar 

  11. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  12. 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.

    Article  PubMed  CAS  Google Scholar 

  13. 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.

    Article  PubMed  CAS  Google Scholar 

  14. Larson CR. The midbrain periaqueductal gray: a brainstem structure involved in vocalization. J Speech Lang Hear Res. 1985;28(2):241–9.

    Article  CAS  Google Scholar 

  15. 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.

    Article  PubMed  CAS  Google Scholar 

  16. Goodson JL. The vertebrate social behavior network: evolutionary themes and variations. Horm Behav. 2005;48(1):11–22.

    Article  PubMed  PubMed Central  Google Scholar 

  17. 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.

    Article  PubMed  CAS  Google Scholar 

  18. Davis PJ, Zhang SP, Winkworth A, Bandler R. Neural control of vocalization: respiratory and emotional influences. J Voice. 1996;10(1):23–38.

    Article  PubMed  CAS  Google Scholar 

  19. 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.

    Chapter  Google Scholar 

  20. Shiba K, Satoh I, Kobayashi N, Hayashi F. Multifunctional laryngeal Motoneurons: an intracellular study in the cat. J Neurosci. 1999;19(7):2717–27.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. 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.

    Article  PubMed  CAS  Google Scholar 

  22. 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.

    Article  PubMed  CAS  Google Scholar 

  23. 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.

    Article  PubMed  CAS  Google Scholar 

  24. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. 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.

  26. Marquis JM, Lettenberger SE, Kelm-Nelson CA. Early-onset Parkinsonian behaviors in female Pink1−/− rats. Behavioral Brain Research. 2020;13(377):112175.

    Article  CAS  Google Scholar 

  27. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 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.

  29. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. 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.

    Article  PubMed  CAS  Google Scholar 

  31. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  33. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  34. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. 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.

    Google Scholar 

  37. 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.

    Article  PubMed  CAS  Google Scholar 

  38. 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.

    CAS  Google Scholar 

  39. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  41. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  43. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. 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.

    Article  PubMed  CAS  Google Scholar 

  45. 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.

    Article  PubMed  CAS  Google Scholar 

  46. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. 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.

    Article  PubMed  CAS  Google Scholar 

  48. 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.

    Article  PubMed  CAS  Google Scholar 

  49. Eisenstein RS. Interaction of the hemochromatosis gene product Hfe with transferrin receptor modulates cellular Iron metabolism. Nutr Rev. 1998;56(12):356–8.

    Article  PubMed  CAS  Google Scholar 

  50. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Liddell JR, White AR. Nexus between mitochondrial function, iron, copper and glutathione in Parkinson's disease. Neurochem Int. 2018;117:126–38.

    Article  PubMed  CAS  Google Scholar 

  52. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Mitchell RM, Lee SY, Simmons Z, Connor JR. HFE polymorphisms affect cellular glutamate regulation. Neurobiol Aging. 2011;32(6):1114–23.

    Article  PubMed  CAS  Google Scholar 

  54. Jurgens U. The role of the periaqueductal grey in vocal behaviour. Behav Brain Res. 1994;62(2):107–17.

    Article  PubMed  CAS  Google Scholar 

  55. 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.

    Article  PubMed  CAS  Google Scholar 

  56. 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.

    Article  PubMed  Google Scholar 

  57. Jurgens U, Richter K. Glutamate-induced vocalization in the squirrel monkey. Brain Res. 1986;373(1–2):349–58.

    Article  PubMed  CAS  Google Scholar 

  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.

    Article  PubMed  CAS  Google Scholar 

  59. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  60. Nickols HH, Conn PJ. Development of allosteric modulators of GPCRs for treatment of CNS disorders. Neurobiol Dis. 2014;61:55–71.

    Article  PubMed  CAS  Google Scholar 

  61. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  62. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. 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.

    Article  CAS  Google Scholar 

  64. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. Guide for the Care and Use of Laboratory Animals. Washington DC: National Academy of Sciences.; 2011.

  66. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  67. Aken BL, Ayling S, Barrell D, Clarke L, Curwen V, Fairley S, et al. The Ensembl gene annotation system. Database. 2016;2016.

  68. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26.

  70. Reiner A, Yekutieli D, Benjamini Y. Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003;19(3):368–75.

    Article  PubMed  CAS  Google Scholar 

  71. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Subhash S, Kanduri C. GeneSCF: a real-time based functional enrichment tool with support for multiple organisms. BMC Bioinformatics. 2016;17(1):365.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references


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.

Author information

Authors and Affiliations



CKN conceived of the study, designed, analyzed, and interpreted the data. SG performed the statistical analysis, interpreted the data, and was a major contributor in writing the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Cynthia A. Kelm-Nelson.

Ethics declarations

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

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Gene modules identified by WGCNA in male rats.

Additional file 2: Table S2.

Gene modules identified by WGCNA for female rats.

Additional file 3: Table S3.

Gene modules identified by WGCNA for all rats.

Additional file 4: Table S4.

Relationships between WGCNA identified modules and female vocal behaivor.

Additional file 5: Table S5.

Relationships between WGCNA identified modules and female vocal behaivor.

Additional file 6: Table S6.

Male RNA-Seq data used to analyze differential gene expression using EdgeR.

Additional file 7: Table S7.

Female RNA-Seq data used to analyze differential gene expression using EdgeR.

Additional file 8: Table S8.

Results of a differential gene expression in female rats.

Additional file 9: Table S9.

Results of a differential gene expression in female rats.

Additional file 10: Table S10.

Top 500 common up- and down-regulated genes in male and female Pink1−/− rats.

Additional file 11: Table S11.

STRING Gene List.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: