Whole-exome sequencing of 228 patients with sporadic Parkinson’s disease

Parkinson’s disease (PD) is the most common neurodegenerative movement disorder, affecting 1% of the population over 65 years characterized clinically by both motor and non-motor symptoms accompanied by the preferential loss of dopamine neurons in the substantia nigra pars compacta. Here, we sequenced the exomes of 244 Parkinson’s patients selected from the Oxford Parkinson’s Disease Centre Discovery Cohort and, after quality control, 228 exomes were available for analyses. The PD patient exomes were compared to 884 control exomes selected from the UK10K datasets. No single non-synonymous (NS) single nucleotide variant (SNV) nor any gene carrying a higher burden of NS SNVs was significantly associated with PD status after multiple-testing correction. However, significant enrichments of genes whose proteins have roles in the extracellular matrix were amongst the top 300 genes with the most significantly associated NS SNVs, while regions associated with PD by a recent Genome Wide Association (GWA) study were enriched in genes containing PD-associated NS SNVs. By examining genes within GWA regions possessing rare PD-associated SNVs, we identified RAD51B. The protein-product of RAD51B interacts with that of its paralogue RAD51, which is associated with congenital mirror movements phenotypes, a phenotype also comorbid with PD.

Scientific RepoRts | 7:41188 | DOI: 10.1038/srep41188 when very large cohorts (> 10,000 samples) are investigated will there likely be sufficient power for such individual rare disease-associated variants to be revealed 5 .
Variants that occur more frequently in the patient cohort than in controls could still be of interest for predicting PD susceptibility despite individually not being genome-wide significant. For example, individual genes, or else multiple genes encoding functionally interacting proteins, might be found that harbour an unusually large number of such sub-genome-wide significant variants. These genes would then add to our understanding of the molecular deficiencies contributing to neurodegeneration in PD.
Here, we analyse exonic variants found in 228 individuals with sporadic PD, all recruited from a relatively homogeneous UK population. We apply multiple approaches at the level of the variant, the gene and the pathway to analyse whether -considered individually or together -these variants can contribute to PD diagnosis.

Results
Exomes of 244 Parkinson's disease (PD) patients selected from the Oxford Parkinson's Disease Centre Discovery Cohort were sequenced. After quality control, and after one individual with relatively strong identity-by-descent (> 0.1875) was discarded, 228 PD patient exomes were considered for analysis and compared with 884 control exomes from UK10K cohort and matching our control in term of ancestry (Method). 94,369 sites with base calls in at least 99% of samples were retained for analysis. Eight individuals known to carry variants in GBA (N370S and L444P) or LRRK2 (R1441C and G2019S) were included as positive controls for mutation detection through exon resequencing. All known variants were detected, although the single GBA L444P variant call failed to exceed the VQSLOD threshold. We called a total of 94,369 single nucleotide variants (SNVs), of which 39,569 had a minor allele frequency of less than 1% (in controls), affecting 10,252 genes. Of 94,369 SNVs, 21,235 were non-synonymous (NS) affecting 9444 protein coding genes, including 313 that caused the gain or loss of a stop codon and 6,541 that were considered by Polyphen-2 6 as being damaging ( Supplementary Fig. 1).
No SNV/gene reached genome-wide significance for PD association. Fisher's exact test and logistic regression using between 1 and 10 covariates were performed using PLINK 7 (http://pngu.mgh.harvard.edu/purcell/plink/). The best-fitting quantile-quantile (Q-Q) plot was obtained with the use of two covariates in logistic regression and the control sets listed in the Methods (Figs 1A and 2). As expected, given the low number of cases and the multiple testing correction required, no substitution achieved significance of association with PD status (p-value upper threshold 2.4 × 10 −6 ; Fig. 2). We next used SKAT 8 to perform a burden test which found no single gene showing a significant association of multiple SNVs with PD risk. (Fig. 1B and Supplementary Fig. 2). After performing burden analyses, we again failed to find a significant enrichment in the number of NS per individual between cases and controls (NS cases = 6639 vs NS controls = 6650, p-value = 0.34).
PD association in previous reported PD candidate genes. We then examined globally whether known genes involved in the familial or sporadic forms of PD were enriched in PD associated NS SNV (Method). For 15 genes involved in familial PD 9 (Supplementary Table 1), we did not find a significant enrichment of PD-associated NS SNVs (10/15 genes had at least one NS SNV; p-value enrichment based on a logistic regression test = 0.054 and p-value enrichment based on the SKAT test = 1). However, for 329 genes harboured within 26 GWA intervals (Methods), a significant 1.2-fold enrichment in PD-associated NS SNVs was achieved (163/329 genes had at least one NS SNV, p-value enrichment based on logistic regression test = 0.036).
Next, we looked amongst 15 genes involved in familial and 329 genes harboured in 26 GWA intervals (Fig. 1, Tables 1 and 2) for NS SNVs with nominal p-value associations of < 0.05 with PD. We identified two PD associated NS SNVs with MAF < 1% within genes harboured in PD risk GWA intervals, each with an odds ratio > 3.5. Specifically, these were (1) rs34094401, a missense mutation (L172W) in the RAD51B gene encoding a protein essential for DNA repair (18 cases versus 20 controls are heterozygous for this variant) lying in the PD GWA interval defined by the lead SNP rs1555399, and (2) rs41309351, a missense mutation (R480W) in the CPXM1 gene encoding for a carboxylpeptidase (16 cases versus 12 controls are heterozygous for the variant) lying in the PD GWA interval defined by the lead SNP rs55785911. Truncating and missense variants in the RAD51B-interacting paralog RAD51 have been reported to cause congenital mirror movement syndrome 10,11 . By examining the frequencies for rs34094401 in an independent cohort including 536 PD patients and 260 controls matching PD patients for age and genotyped on HumanCoreExome chip from Illumina 12 , we observed a similar but non-significant effect at the RAD51B locus (p-value = 0.14): MAF cases = 0.01306 (14 cases heterozygous) vs MAF controls = 0.005769 (3 controls heterozygous). We noted that the control minor allele frequency of these two SNPS were consistent with those reported by the Exome Aggregation Consortium (ExAC) in the European (Non-Finnish) population 13  Functional analyses of PD-associated SNVs. Although we did not find genome-wide significant results for PD association, genes whose proteins contribute to a common cellular process or pathway could contain an enrichment of SNVs that are sub-genome-wide significant for PD association. For this analysis, we examined whether the top N genes with most PD-associated NS SNVs (N = 100, 200, 300) were enriched in particular Gene Ontology terms 14 . Only when considering the top 300 genes were any significant enrichments found, specifically extracellular matrix part (GO:0044420; 1.8-fold enrichment; q-value = 0.02) and extracellular matrix disassembly (GO:0022617; 1.8-fold enrichment; q-value = 0.02) (Supplementary Table 2). Although only significant when considering such a large number of genes, the genes annotated with a role in extracellular matrix disassembly include those encoding the metalloproteinases MMP-7 and MMP-8, and the reduced expression of metalloproteinases has been noted in PD post-mortem brain tissue 15 while a polymorphism in MMP-9 has been associated with PD and amyotrophic lateral sclerosis 16 .

Discussion
In this study, we compared the exomes of 228 PD cases with 884 controls exomes drawn from the UK10K 17 study. We performed association tests both at the level of single-nucleotide variants and at the gene level but found that no variant was significantly associated with PD after applying a multiple-testing correction. While power for discovery is clearly the major issue, there may be additional considerations affecting our study. As the heritability of PD is not high (34-40%), the identification of causative variants or risk alleles for PD may be more difficult than for some other, common neurological disorders 18,19 . Furthermore, as our controls are not aged-matched, we do not know whether control individuals will go on to develop PD. Indeed, for example the LRRK2-G2019S  mutation has a penetrance of 28% at age 59 years, 51% at age 69 years, and 74% at age 79 years 20 . Burden tests or collapsing methods, such as SKAT, might be expected to improve power over single-marker tests, but they still require thousands of samples to reach acceptable statistical power 8,21. We considered different subsets of the UK10K 17 controls for this study and the best-fitting Q-Q plot was obtained using a set of six UK10K data sets 17 (see Methods). Notably for other studies seeking exome controls, the use of the two largest UK10K cohorts, COHORT_TWINS (~1,700 samples) and COHORT_ALSPAC (740 samples), produced inflated p-values that could not be controlled by filtering or covariates. Some of the minor allele frequencies from these two cohorts that were responsible for this effect were very different from the allele frequencies in other controls, suggesting a population stratification problem. While DNA sequencing enables the direct detection of causal mutations and may supersede microarray-based genotyping methods, some of the limitations of GWAS still apply to these studies. Thousands of cases and controls are required to robustly associate rare variants with a common disease and population stratification has to be controlled for 5,22 .
NS SNVs most associated with PD were enriched both in genes whose proteins have roles associated with the extracellular matrix and in genes lying within PD-associated GWA intervals. Furthermore, we identified a nominally-associated SNV in the RAD51B gene which lies within a region of the genome associated with PD through GWA. Variants in the paralog RAD51, whose protein product interacts with that of RAD51B, have been associated with congenital mirror movement disorders (OMIM 614508) and mirror movement phenotypes are also associated with PD 23 . Thus, RAD51B represents an interesting candidate gene for further examination within this GWA locus.   25 , following which patients were diagnosed as having PD, with an estimated clinical probability for this diagnosis being made by the research neurologist. Atypical parkinsonian features were systemically screened for using the NINDS Parkinson's tool, and patients with secondary parkinsonism due to head trauma or medication use, or features of atypical parkinsonism syndromes, were excluded from the study. Subjects were also excluded from subsequent analysis if they had a < 80% baseline probability of manifesting idiopathic PD. Clinical assessments were performed while the subject was taking their normal medications, and in the 'on' state. The study was undertaken with the understanding and written consent of each subject, with the approval of the NHS South Central Berkshire Research Ethics Committee (study reference 10/ H0505/71), and in compliance with national legislation and the Declaration of Helsinki. 244 consecutive PD subjects from the Discovery cohort were recruited to this exome study, of whom 224 (90%) were selected on the following basis: Control Exome samples. As our PD exome patients have ben recruited in UK, we used public whole-exome sequencing data of individuals matching our PD cohort in term of ancestry and thus with UK origin. The following exome sample datasets were used from the UK 10K dataset: NEURO_IOP_COLLIER (152 samples), NEURO_ABERDEEN (317 samples), NEURO_EDINBURGH (213 samples), RARE_NEUROMUSCULAR (72 samples), RARE_THYROID (18 samples) and NEURO_IMGSAC (112 samples), making a total of 884 control exomes. The two largest UK10K 17 cohorts, COHORT_TWINS and COHORT_ALSPAC, were not used as they contained minor allele frequencies that differed from those of other control cohorts, indicating a possible population stratification problem. Furthermore, to avoid inflation of low p-values due a stratification issue with particular subsets of the UK10K controls, best-fitting Q-Q plot was obtained using the set of six control data sets.
Exome sequencing. 120 μ l of 25 ng/μ l of DNA was fragmented to an average size of 150 bp (75-300 bp), and subjected to Illumina DNA sequencing library creation using Bravo automated liquid handling. Adapter ligated libraries were amplified via 6 cycles of PCR. 500 ng/ul of each amplified library was hybridized to RNA baits (SureSelect Human All Exon 50 Mb) and post sequence capture processed following the Agilent SureSelect protocol. Enriched libraries were post capture indexed using a unique DNA barcode via 12 cycles of PCR. Equimolar pools of 8 libraries were sequenced using the Illumina HiSeq machine using the 75 base paired end protocol. Data processing. Sequencing reads from each lane were aligned to the human reference genome used in the 1000 Genomes phase 2 (NCBI build 37 + decoy v5) 26 using BWA version 0.5.9-r16 and the parameters '-q 15 -t 6' . The GATK 27 'IndelRealigner' was used to realign reads near known indels [ftp://ftp.1000genomes.ebi.ac.uk//vol1/ ftp/technical/reference/phase2_mapping_resources/]. The BAM files were then re-sorted and quality scores were recalibrated using GATK 'TableRecalibration' 27 ignoring all SNP positions in dbSNP v135 28 . Finally, SAMtools 'calmd' was used to recalculate MD/NM tags in the BAM files. All lanes from the same library were then merged into a single BAM file using Picard tools 29 and PCR duplicates were marked using Picard 'MarkDuplicates' . Finally, the library BAM files were merged into a single BAM per individual.
Variant calling. SNP and indel calls were made using samtools 30  Additional filtering. Additionally, the following filters were applied to individual genotype calls based on their depth or genotype quality: IndelsMinSampleDP,Description = "Genotypes set to. for samples with DP < 4 (FORMAT/DP)" IndelsMaxSampleDP,Description = "Genotypes set to. for samples with DP > 2000 (FORMAT/DP)" SNPsMinSampleGQ,Description = "Genotypes set to. for samples with GQ < 20 (FORMAT/GQ)" IndelsMinSampleGQ,Description = "Genotypes set to. for samples with GQ < 60 (FORMAT/GQ)".
Annotation. The calls were annotated with dbSNP 137 rsIDs and population allele frequencies from the 1000 Association analysis. VCF files were converted to PED format (format of file used by PLINK (Purcell et al. 2007)) with VCFtools 87 33 [-remove-indels -remove-filtered-all]. We performed single variant association analyses by using logistic regression test incorporating two first ancestry covariate. We evaluated evidence for association at the genetic level by using the sequence kernel association test (SKAT) 8 that allows for risk and protective effects in the same gene". We used asymptotic p-values and incorporated ancestry covariate (two first Principal Components). SKAT assumes that rare variants are more likely to be causal variants with large effect sizes and assigns weight to the variant according to minor allele frequency (we used the default frequency weight (1 and 25 as the two beta distribution shape parameters, to up-weight lower frequency mutations)). As power was clearly an issue in this study, we did not repeat these analyses by functional category of variants.
PD GWA intervals. The 26 PD GWA intervals were defined: (i) by identifying the most distant pair of SNPs with either (i) an r 2 > 0.5 with the lead SNP reported in http://pdgene.org 3,34 using the haplotypes of 1000 Genome Project or (ii) an r 2 > 0.5 computed by considering the European haplotypes only, and then adding an additional 250 kb on either side of the interval to include genes that may be regulated by GWA-associated regulatory variants.

Gene Sets enriched in genes with PD-associated Non-synonymous SNVs.
To evaluate if a set of genes was enriched in PD associated NS SNVs, we computed the sum of -log 10 of p-value association of logistic regression test of NS SNVs within these genes and compared this measure with this got with random set genes matching each gene of test gene set for the number of NS SNVs. We performed an analogue test by replacing the p-value association of logistic regression by p-value associated with SKAT test.
Tissue enrichment analysis. FPKM= > Feature Scaling= > Euclidean Normalization (such that the square root of the sum of the square of expression levels is one). Compare the sum of expression value for PD risk gene with a set of random gene (match for CDS length).
Overlap with PD GWAS intervals. 26 genome-wide significant PD intervals were acquired from http:// pdgene.org 3,34 . Intervals were extended at each end by 500kb. The number of the top N genes (10 ≤ N ≤ 100) most significantly associated by SKAT with PD located within either set of intervals was compared with the number of randomly sampled length-matched genes in these intervals.