Transcriptomic Profiling of Extracellular RNAs Present in Cerebrospinal Fluid Identifies Differentially Expressed Transcripts in Parkinson’s Disease

Human CSF samples

Frequency matching was used to recruit 27 PD patients and 30 age- and sex-matched healthy controls at the Pacific Northwest Udall Center of Excellence for Parkinson’s Disease Research (PANUC) and the Alzheimer’s Disease Research Center (ADRC) at the University of Washington (UW) and Veterans Affairs (VA) Puget Sound Health Care System (Seattle, Washington, USA) as previously described [8–10]. All subjects underwent the following evaluations: medical history, physical and neurological examinations, laboratory tests, and neuropsychological assessments. The inclusion and exclusion criteria for study subjects have been previously described [8–11]. Briefly, PD patients met UK PD Society Brain Bank clinical diagnostic criteria for PD [12], except that having “more than one affected relative” was not considered an exclusion criterion. Control subjects were the patients’ spouses or community volunteers in good health. They had no signs or symptoms suggesting cognitive impairment or neurological disease; all had a Mini Mental Status Exam (MMSE) score between 28 and 30, a Clinical Dementia Rating (CDR) score of 0, and New York University paragraph recall scores (immediate and delayed) of >6. All procedures were performed in strict compliance with guidelines for human experimentation and all participants underwent detailed informed consent procedures. Subjects provided consent in writing in accordance with procedures approved by the institutional review boards at the UW and VA Puget Sound Health Care System.

CSF was collected by lumbar puncture as previously described [9]. Similar sample collection protocols and quality control procedures were followed at all participating centers, in particular, the use of polypropylene collection and storage tubes, rapid separation into single use aliquots, and freezing of CSF samples, to minimize potential site variability. De-identified samples were aliquoted and stored at –80°C, coded as group “A” or “B”, and then transferred to the University of Miami for RNA analysis. The UW CSF data dictionary is shown in Supplementary Table 1.

Library preparation of the CSF transcriptome and sequencing

Library preparation was optimized to achieve successful sequencing runs from low input RNA by scaling up reverse transcription and PCR reaction volumes without increasing PCR amplification cycles. Extracted RNA was between 550 t0 3600 picograms as shown in Supplementary Figure 1. We utilized NEB-based directional RNA-Seq (Ipswich, MA, USA) sample preparation kit to prepare all libraries. End repair of cDNA libraries was performed and 5’ and 3’ adaptors were ligated, sequentially. Eluted cDNA was enriched and remaining RNA products were digested with USER excision and PCR amplification. During enrichment, each sample was indexed with one unique index primer from a set of 12 different index primers for Illumina multiplex sequencing. RNA was reverse transcribed to synthesize first strand cDNA, followed by a purification step, using Agencourt RNA Clean XP beads (Beckman Coulter, Brea, CA, USA). Second strand cDNA synthesis was performed and libraries were amplified (12 cycles), followed by purification using 1.8 X Agencourt AMPure beads (Beckman Coulter, Brea, CA, USA). Finally, three uniquely indexed libraries were pooled and sequenced on a HigSeq2000 Illumina sequencer. A quality control assessment was performed, using the Bowtie software package by estimating genome coverage, percent alignment and nucleotide quality. Raw sequencing data was transformed to FastQ format and stored on a dedicated server. All raw sequencing data and corresponding codes are uploaded onto the Parkinson’s Disease Biomarker Program database (PDBP) Data Management Resource (PDBP DMR) (http://pdbp.ninds.nih.gov) and are freely available to registered users.

Sequencing read processing and genome mapping

Genome mapping was performed using the Spliced Transcripts Alignment to Reference (STAR) software [13], based on the latest available versions of STAR 2.1.3 and the most recent human transcriptome (reference hg19) downloaded from the UCSC database (Santa Cruz, CA, USA). Transcript assembly was processed and counts per million (CPM) values for each aligned sequence (known and novel) were determined. To ensure that reads were indeed novel, we applied an initial filter based on expression levels to eliminate low abundance single-exon transcripts to avoid mislabeling of genomic DNA contaminated reads.

RNA sequencing analysis and identification of differentially expressed transcripts

A trimming algorithm based on Trim Galore software (Babraham Institute, Babraham, United Kingdom) was used to trim reads for adapter removal. EdgeR [14], a Bioconductor software package for differential expression analysis of replicated count data, was used to correct for multiple testing (false discovery rate (FDR) cutoff of <0.1) and for the identification of differentially expressed transcripts based on CPM values [14].

Functional annotation and pathway analysis of differentially expressed transcripts (DETs)

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DETs was performed using the Database for Annotation Visualization and Integrated Discovery (DAVID) online analytical tools [15–17]. Protein-protein interaction networks of DETs were constructed via the search tool for the retrieval of interacting genes/proteins (STRING; http://www.stringdb.org/) and the similarity between the PPI networks and the KEGG pathways were calculated via EnrichNet [18, 19].

Quantitative real time PCR (qRT-PCR) validation of RNA-seq data

The cDNA synthesis protocol is described in Supplementary Figure 2. The majority of RNA in the CSF is short in length, most likely the result of RNA fragmentation. We synthesized cDNA and utilized SYBR green-based qRT-PCR reactions with transcript-specific forward primers and a universal PCR (reverse) primer (Supplementary Table 3) that anneals to the 5’ end of the tagged cDNA (Supplementary Figure 3). Small CSF sample volumes limited qRT-PCR gene validation to a subset of candidate sequencing tags; therefore 21 transcripts (Supplementary Table 3) were selected for validation based on their relevance to known PD genes in the PPI analysis.

We performed Sanger sequencing to validate the specificity of qRT-PCR products. Raw data was normalized to the geometric mean of the 18 S gene that met the criteria of having similar expression between sample groups. Relative expression was compared between PD and control groups using the student’s t-test. All measurements were performed in triplicate and data were analyzed using the ΔCt method. Statistical analyses were performed using SPSS 21 for Mac OSX (IBM, Armonk, NY, USA).

High-throughput sequencing of the CSF transcriptome

We analyzed RNA transcripts from the CSF of 27 PD patients and 30 controls. CSF RNA-seq data were analyzed using STAR to align reads to the human genome, generating an average of 84×10⁶ paired-end reads that aligned to the reference human genome as well as to annotated genes. The average number of trimmed reads was 81×10⁶ per sample, of which 82.7% aligned to the human reference genome (UCSC_hg19) (Supplementary Table 2). Our analysis focused on RNA fragments that were reproducibly present in at least 65% of samples in each group and CPM >1. After filtering, a total of 3521 fragmented transcripts were detected: protein-coding (2862 transcripts; 81.3% of total transcripts), non-coding RNAs (458 transcripts; 13% of total transcripts) and pseudo genes (201 transcripts; 5.7% of total transcripts, Supplementary Table 4).

Differentially expressed transcripts (DETs) in the CSF of PD patients

Differential expression analysis was performed using p <  0.05 to identify a total of 467 DETs between PD patients and controls. Using p-value cutoffs of p <  0.001, 0.01 <  p <  0.001 and 0.01 <  p <  0.05 as shown in the volcano plot narrowed down the list of DETs to 207, 103 and 157, respectively (Fig. 1).

Differential expression analysis of CPM values using p <  0.05, FDR <0.1 and |log2 FC|>1 as cutoff criteria produced a total of 201 transcripts differentially expressed between PD patients and controls. Among these DETs, 92 transcripts were upregulated and 109 transcripts were downregulated in PD patients. Of the 201 DETs, 142 (70.6%) were related to protein-coding transcripts (Supplementary Table 5). Interestingly, non-coding RNA and pseudo genes are more abundant (29.4% vs 18.1% ; p = 0.001, Relative Risk = 1.8, 95% CI: 1.34 –2.41) in DETs as compared to non-DETs. Protein-coding transcripts are present in a lower proportion in DETs as compared to non-DETs (70.6% vs 81.9% , Fig. 2).

Functional annotation clustering of differentially expressed transcripts

To determine the relationship between differentially expressed transcripts identified in RNA-seq to relevant PD processes and pathways, gene ontology (GO) and pathway enrichment analyses were performed using DAVID online tools. We considered significantly enriched genes as those with the following criteria: p-value <0.05, functional categories with the highest classification stringency and enrichmentscore >1.

Based on the cutoff criteria of |log2 FC|>1, p <  0.05 and FDR <0.1, we screened 201 DETs. Up-regulated transcripts were mostly enriched for chromatin regulators, while down-regulated transcripts were significantly enriched for protein-tyrosine phosphatases, endoplasmic reticulum membrane proteins, and genes that regulate phosphorylation, dephosphorylation and phospholipase activity (Table 1). Details of functional annotation clustering in up and down-regulated DETs are shown in Supplementary Tables 6 and 7.

Construction of protein-protein interaction (PPI) networks

We mapped 201 DETs to the STRING database (the hub protein was selected according to the node degree) and screened for significant interactions (>0.7). These relationships were integrated to construct interaction networks among interacting proteins (Fig. 3) and the similarity between PPI networks and KEGG pathways were calculated via EnrichNet [18, 19]. This analysis detected the following proteins: EZH2, DNMT1, NDUFV2, MAPK11, CCR3, LPAR1, SSTR5, PTPRC, and BMP7 that interact with known PD genes PARK7, LRRK2, and SNCA. Although these three PD transcripts were not in our DETs, because of their direct interaction with our DETs, were added to the module (Fig. 3) as these protein-coding genes play an important role in the pathological process of Parkinson’s disease.

Validation of CSF sequencing

Target transcripts included 19 genes from the PPI network (Fig. 3) and two non-coding RNAs that are differentially expressed. The differential expression observed in RNA-seq between the two groups was validated in 13 out of 17 transcripts (76%) by qRT-PCR (Table 2). Our comparison of RNA-Seq and qRT-PCR showed good agreement in both the direction and magnitude of fold change for candidate transcripts. Log2 fold changes of these transcripts in qRT-PCR correlated with RNA-seq (p = 0.01, r² = 0.65). AC010127.3 is an antisense RNA to SCN9a that is down-regulated in qRT-PCR while SCN9a is significantly up-regulated. UC001lva.4 (p = 0.01, Log2FC = –1.6) and AC079630 (p = 0.001, Log2FC = –6.72) are two lncRNAs on the LRRK2 locus that are significantly down-regulated in qRT-PCR while the LRRK2 gene itself is not significantly altered in expression. NDUFV2 and LRPPRC (mitochondrial biogenesis and maintenance factors), are both significantly down-regulated.