A Genome Wide Association Study Links Glutamate Receptor Pathway to Sporadic Creutzfeldt-Jakob Disease Risk

We performed a genome-wide association (GWA) study in 434 sporadic Creutzfeldt-Jakob disease (sCJD) patients and 1939 controls from the United Kingdom, Germany and The Netherlands. The findings were replicated in an independent sample of 1109 sCJD and 2264 controls provided by a multinational consortium. From the initial GWA analysis we selected 23 SNPs for further genotyping in 1109 sCJD cases from seven different countries. Five SNPs were significantly associated with sCJD after correction for multiple testing. Subsequently these five SNPs were genotyped in 2264 controls. The pooled analysis, including 1543 sCJD cases and 4203 controls, yielded two genome wide significant results: rs6107516 (p-value=7.62x10-9) a variant tagging the prion protein gene (PRNP); and rs6951643 (p-value=1.66x10-8) tagging the Glutamate Receptor Metabotropic 8 gene (GRM8). Next we analysed the data stratifying by country of origin combining samples from the pooled analysis with genotypes from the 1000 Genomes Project and imputed genotypes from the Rotterdam Study (Total n=12967). The meta-analysis of the results showed that rs6107516 (p-value=3.00x10-8) and rs6951643 (p-value=3.91x10-5) remained as the two most significantly associated SNPs. Rs6951643 is located in an intronic region of GRM8, a gene that was additionally tagged by a cluster of 12 SNPs within our top100 ranked results. GRM8 encodes for mGluR8, a protein which belongs to the metabotropic glutamate receptor family, recently shown to be involved in the transduction of cellular signals triggered by the prion protein. Pathway enrichment analyses performed with both Ingenuity Pathway Analysis and ALIGATOR postulates glutamate receptor signalling as one of the main pathways associated with sCJD. In summary, we have detected GRM8 as a novel, non-PRNP, genome-wide significant marker associated with heightened disease risk, providing additional evidence supporting a role of glutamate receptors in sCJD pathogenesis.


Introduction
Sporadic Creutzfeldt-Jakob disease (sCJD), although rare, with a yearly incidence of one to two cases per million, is the most common form of human prion disease. This group of disorders is characterized by spongiform changes in the brain, as well as accumulation of misfolded, often protease-resistant, conformers (PrP Sc ) of the normal prion protein (PrP C ). The PrP C gene (PRNP) plays a central role in prion disease susceptibility. Expression of PrP C is indispensable for disease transmission [1] and the polymorphism coding for methionine (M) or valine (V) at codon 129 (PRNP M129V) has been linked to disease risk [2]. Homozygosity at PRNP M129V has been consistently associated to sCJD, being one of the strongest common genetic risk factors reported for neurodegenerative diseases. The remarkable disease-determining effect of this PRNP polymorphism is observed in variant CJD, a subtype acquired from dietary exposure to bovine spongiform encephalopathy [3], where all definite and probable clinical cases studied to date have been PRNP129MM [4].
Similar to other diseases, several genetic association studies with candidate genes have been performed on sCJD susceptibility [5]. Only one previous genome wide association study (GWAS) of sCJD risk has been published to date, showing that the PRNP locus was strongly associated with disease risk, specifically with rs1799990 (PRNP M129V) [6].
To further scrutinize genomic variations related to sCJD risk we have performed a threestage GWAS encompassing a total of 1,543 sCJD cases and 4,203 controls, as well as a metaanalysis encompassing data from the 1000 Genomes Project and imputations from the Rotterdam Study.

Results
Demographic and clinical features of the sCJD case populations are shown in S1 Table. The Q-Q plots for autosomal and X chromosome SNPs are given in S1 Fig; during discovery (stage one) the genomic inflation factor λ was 1.053 for autosomal and 1.057 for X chromosome SNPs.
S2 Table shows the top 100 SNPs associated with sCJD at discovery stage, sorted by allelic differences p-value. A total of 23 SNPs were taken forward to replication. In stage two, we successfully genotyped 22 of the 23 SNPs with the Sequenom iPLEx GOLD platform in an independent population of 1,109 samples of sCJD. Table 1 shows that five SNPs remained significant  Table 2 shows differences in allelic frequencies in the discovery, replication and pooled population sets including 1,543 sCJD cases and 4,203 controls. The A allele of rs6951643 was associated with a 1.27 fold increased risk of sCJD (95%CI = 1.17-1.38), and this effect was consistently observed across all tested populations (Table 3). Rs6107516 (p-value = 3.00x10 -8 ) and rs6951643 (p-value = 3.91x10 -5 ) remained as the two most statistically significant associations in a meta-analysis of the discovery data adjusted for PCAs and analysis of replication samples stratified by country. This meta-analysis included discovery and replication samples plus data from the 1000 Genomes Project and imputed genotypes from the Rotterdam Study (Total n = 12967), ( were 12 tagging GRM8 (S2 Table). Pathway analysis using IPA showed that glutamate receptor signalling was the most over-represented canonical pathway within our top results (p-value = 8.01x10 -5 ) (Fig 1). Analysis with ALIGATOR software showed (S3 Table) that the GO category glutamate receptor activity was the 4th most represented amongst our top results. Three genes from the GO category "glutamate receptor pathway" (S2 Table) were within our top 100 SNPs (GRM8, GRIN2B and GRIA1). Fine mapping in search for functional variants linked to rs6951643 was attempted sequencing the 11 exons of GRM8 gene and the corresponding intronic flanking regions in 96 sCJD patients. We found 8 intronic SNPs (rs73231278, rs62468898, rs1008274, rs111546739, rs17685327, rs6951643, rs2074012, rs35648111) and one exonic non-coding SNP (rs34182595). Functional prediction analysis of these genetic variants by FuncPred or RESCUE-ESE showed no indication of regulatory implications. In the Gtex eQTL database we found no eQTL associated to any of the GRM8 SNPs.
There is no co-expression of GRM8 and PRNP in blood eQTLs in Genenetwork. However, based on PITA prediction we found that rs34182595, which is an insertion/deletion SNP in exon 11 in partial LD with rs6951643 (D' = 0.41; r 2 = 0.16), is located in a target site for micro-RNAs (miR-103 and miR-107), a finding that has been previously reported [7].
We then performed an immunohistochemical study in brain samples from 48 sCJD patients. Immunoreactivity for mGluR8 was observed in neurons and microglial cells (both in the white and grey matter) While in neurons we did not observe obvious differences we found a trend towards higher combined score of mGluR8 immunoreactivity in microglia related to the A allele at rs6951643 (Fig 2); however, these differences were not statistically significant (ordinal regression A-carriers versus A-non-carriers adjusting by c129 and levels of microglia in temporal cortex; p-value = 0.093). S3 Fig shows semi-quantitative assessments of mGluR8 expression in microglia (0 to 3) in the temporal region across the three rs6951643 genotypes. We found no association between rs6951643 genotypes and patient's age of disease onset or disease duration.

Discussion
In this study we report a non-PRNP genetic risk variant for sCJD, GRM8. Moreover, building on findings from previous studies [8], pathway analyses yielded glutamate receptor signalling as one of the main pathways linked to sCJD pathogenesis.
A previous GWAS of prion diseases employed 1,259 sCJD samples, in addition to other disease subtypes and 6,015 shared controls [6]. In the sCJD sub-group and in a meta-analysis including all prion cases only variants in PRNP were found to be significantly associated with disease risk. Several SNPs outside PRNP were identified, but in contrast to our multi-national disease cohort the association with sCJD was not homogeneous across the different geographical groups. Our main novel finding (rs6951643) was not statistically significant in that analysis. One possible explanation for this discrepancy is the fact that both studies have limited power to detect small effects due to the overall restricted case numbers included for study and therefore a SNP with a relatively modest OR of 1.27 might not be detected by all analyses. The experience and insights drawn from GWAS performed in more prevalent disorders, such as Alzheimer's disease, shows that increasing sample size and performing meta-analysis might be essential in order to discover and confirm, relatively small-effect genetic risk factors. However, it is worth re-iterating that in our study, the association with rs6951643 was consistent across the geographically diverse sCJD population tested, showing an over-representation of the A allele in cases, and across the different genotyping methods used.
A limitation of our primary analysis was the fact that we were not able to adjust by country of origin at the replication stage. In order to overcome this caveat we performed a metaanalysis including control data from the 1000 Genomes Project and imputed data from the Rotterdam Study. In the meta-analysis of the stratified analysis by country of origin the association of both SNPs became less significant suggesting that some signal might be caused by population stratification. Despite the fact that both variants were nominally significant in the metaanalysis of the replication data, the PRNP variant was the only genome wide significantly associated to sCJD in the pooled meta-analysis. However, this meta-analysis is not perfect either, as ideally the cases and controls should originate from the same population and in our analysis we attempted to match on country of origin (see S4 Table) and this approach is statistically less powerful because of the uneven distribution of cases and controls, as might reflect the fact that the PRNP SNP p-value is also less significant in the meta-analysis than in the primary analysis. Still we acknowledge that the results are a call for caution and the association between GRM8 genetic variants and sCJD should be confirmed in larger analysis with extensive control for population stratification. GRM8 encodes for mGluR8, a protein that belongs to the metabotropic glutamate receptor family, a family that has recently been linked to the transduction of physiological and cytotoxic signals mediated by PrP C . MGluR1 and mGluR5 have been shown to interact with PrP C , with such associations appearing important for promoting neurite outgrowth [9,10]. Additionally, a recent study demonstrated that mGluR5 coupled with PrP C mediated the cellular toxicity of soluble β-amyloid oligomers [10]. Further, in the APPswe/PS1dE9 Alzheimer's disease mouse model, altered PrP C processing and a selective increase in cortical mGluR1 expression has been reported. The authors hypothesized that complex processing of PrP C in connection with mGluR1 over-expression is triggered by β-amyloid peptides [11]. In the setting of accumulation of PrP Sc , our immunohistochemical assessment of a limited sample of sCJD patients found that carriers of the risk allele at rs6951643 tended to have higher mGluR8 expression in microglial cells compared to non-carriers. Although our functional prediction analysis for rs6951643 was uninformative, it is of interest that a nearby variant in partial LD (rs34182595) is located in a micro-RNA target site and could influence gene expression.
In conclusion, our study has detected a GRM8 genetic variant as a suggestive marker for sCJD risk mapping outside the PRNP region. Our findings provide evidence supporting a role for glutamate receptor signalling pathways in sCJD susceptibility. The involvement of glutamate receptor pathway in sCJD should be addressed in future studies in order to provide new insights in its pathogenesis. Our results underscore the importance of increasing sample sizes in future studies in order to augment the likelihood of detecting additional non-PRNP genetic risk factors for sCJD.

Ethics statement
The present study was conducted according to the revised Declaration of Helsinki and Good Clinical Practice guidelines. A signed informed consent to participate in genetic research was obtained from all participants or patients' relatives. The study was approved by Comité de

Populations and study design
All patients were ascertained by National CJD Surveillance Centers. Only definite or probable sCJD cases, according to accepted classification criteria, were included [12]. Table 5 summarizes the study design and populations included after quality control. All cases and controls were of Caucasian origin. Legal representatives gave written informed consent, and all samples were taken in accordance with the Helsinki declaration.

Genotyping and quality control
At the discovery stage (stage one) patient's samples were characterized using the Affymetrix 500k array at the USC node of the Spanish National Genotyping Center. Out of the initially available 554 sCJD samples there was not enough DNA in 61 (36 samples from the UK and 25 from Germany), leaving 493 for genotyping. Genotypes were called first by the Dynamic Model (DM) and those samples with an overall call rate >93% were subsequently called by a Bayesian Robust linear model with Mahalanobis distance (BRLMM). Samples with BRLMM call rates <95% and with call and discrimination rates of the Modified Partitioning Around Medoids algorithm (MDR-MCR) >10% were excluded (31 patient samples out of 493 did not reach these thresholds, leaving in total 462 available for the analysis). Genotype and quality control of the WTCCC controls have been described elsewhere [13]. The other set of controls (Rotterdam Study) were genotyped at the Genetic Laboratory, Department of Internal Medicine at Erasmus Medical Center (Rotterdam) also using the Afymmetrix 500K array following same protocols as those used for cases. After quality control, 28 patients and 25 controls were excluded resulting in a final discovery population of 434 sCJD patients and 1,939 controls. Replication was attempted in independent sCJD patients (n = 1,201) (stage two) and control series  (n = 2,264) (stage three) with the Sequenom iPLEx GOLD platform. SNPs for replication were selected as follows: after excluding outliers, we chose the 10 SNPs with lowest p-values (when a cluster of SNPs was in LD tagging the same gene we selected only the one with lowest p-value) and the remainder from the top 100 when a) they were in linkage disequilibrium with other top 100 SNPs or b) agreed in direction with the previous sCJD GWAS [6], or c) were connected to pathways of interest based on our previous study [14] (phosphatidylinositol) or based on the pathway analyses performed with the current data (glutamate receptors).

Statistical analysis of genetic data
The statistical analyses were conducted using GenABEL [15]. For the individual SNP analysis, we excluded those SNPs: 1) with call rates <98% within either group (n = 219,182); 2) with minor allele frequency<0.01 (n = 1); 3) with controls not in Hardy-Weinberg equilibrium (False Discovery Rate for unacceptably high individual heterozygosity<1%) (n = 1,300). The quality control further included a check of the sex chromosomes against the reported sex and unexpected sample duplicates (Identity by State, IBS>95%). All quality control was performed with the 'check.marker' function of GenABEL. After quality control, 279,389 out of 499,872 SNPs were included in the analyses. In the discovery analyses (stage one), we conducted a one degree of freedom additive score test between cases and controls coding the presence of the minor genotype 0 for non-carriers, 1 for heterozygous and 2 for homozygous carriers. We controlled subpopulation structure using principal components analysis (PCA). In brief, we selected a random set of 10,000 SNPs and calculated a genomic kinship matrix using pair wise identity by state statistic (function 'ibs'). We then performed PCA analysis (function 'cmdscale'), and adjusted our analysis by the three main IBS matrix principal components (function 'mlreg'). We calculated the genomic inflation factor lambda (λ) for both autosomal and X chromosome SNPs. In the replication phase (stages two and three), we compared allele frequencies in cases and controls by a Chi-squared test implemented in Haploview [16]. We used a Bonferroni correction to adjust for multiple testing, setting the threshold for significance to a p-value of 0.0023 (0.05/22 SNPs successfully genotyped). We combined the gene discovery and replication series of patients for all validated SNPs and used as the criterion for genome wide significance a p-value of 5x10 -8 . We attempted to adjust for country of origin after stage three. We analysed the data from stage two and three with genotypes of the 1000 Genomes Project [17] and imputed genotypes from the Rotterdam Study [18]. We extracted the SNPs from the 1000 Genomes reference set excluding the children and other family members (Version: phase I v3). Additionally we imputed the same SNPs from the Rotterdam Study (imputations of the same phase I v3of the 1000 Genomes). Imputed genotypes had a very high quality score (R 2 > = 0.99). Next we paired our cases and controls by country of origin (S4 Table) and calculated effect estimates of the SNPs per country. The effect estimates of the SNPs from the discovery stage (adjusted by PCAs) and the effect estimates from the replication were meta-analysed using a inverse variance weighted meta-analysis(METAL version released 2011-03-25). [19]. We performed pathway analyses using ALIGATOR, a method for studying groups of genes by testing for over-representation of members of those groups within lists of genes containing significantly associated SNPs from GWA studies [20]. This analysis used 410 SNPs with a p-value for significance of <0.001 out of 115,565 within-gene SNPs from a total of 279,389 available. We found these SNPs were associated with 94 genes from a total of 13,092 genes with GO annotation covered in this study. Additionally we used Ingenuity Pathway Analysis (IPA) www. ingenuity.com) to determine the functional pathways in the genes tagged by our top ranked SNPs. We again selected those genes tagged by SNPs with p-value <0.001 and selected the canonical pathway analysis implemented in IPA software.

Sequencing
We sequenced the 11 exons and intronic flanking areas of the GRM8 gene in 96 sCJD patients using specific primers (S5 Table). The amplification reactions were carried out with 25 ng of genomic DNA and 0.5 units of Taq DNA Polymerase in a volume of 25 μl. The final concentrations of other reactants were: 1x Taq DNA Polymerase Buffer, 0.1 mM dNTPs and 0.1 μM of each primer. The final concentration of MgCl 2 was 2 mM for the amplification of exons 1-7 and 10; 4 mM for exons 8 and 9a; 1.5 mM for exon 9b and 1 mM for exon 11. Additionally, the amplification reactions of exons 1 and 8 were supplied with DMSO 10% and 5%, respectively. The PCR cycling conditions were as follows: initial denaturation at 96°C for 3 min followed by 30 cycles of 96°C for 30 s, annealing temperature (see S5 Table) for 30 s and extension temperature of 72°C for 1 min with a final extension at 72°C for 10 min. A 2 μl aliquot of the amplification reaction was sequenced using 0.1μM of the above primers.

Phenotypic correlations
We performed an immunohistochemical study in brain samples from 48 sCJD patients. The rs6951643 genotype distribution was as follows: 13 AA, 25 AG, and 9 GG. Sections from the hippocampal region CA1 sub-region, temporal cortex and white matter were immunostained for mGluR8. The intensity of mGluR8 (1:50; polyclonal rabbit antibody, Novus Biologicals, Cambridge, UK) immunostaining was evaluated using a scale of 0-3 (0: no; 1: weak; 2: moderate; 3: strong staining). The frequency of mGluR8 positive cells was scored semi-quantitatively using three categories: 1, <10%; 2, 10-50%; 3, >50% and was evaluated to provide information about the relative number of mGluR8 positive cells within the tissue. The product of these two values (intensity and frequency scores) was used to give the overall scores (total scores). In adjacent sections microglial activation using immunostaining for HLA-DP, DQ, DR (clone CR3/ 43, 1:00, monoclonal mouse), and spongiosis were also estimated semi-quantitatively. The analysis was performed (GGK) blinded to rs6951643 genotype. Ordinal regression was employed to test for association between rs6951643 genotypes and brain semi-quantitative expression of mGluR8. We adjusted by disease duration, PRNP M129V genotype, gliosis and spongiosis. We correlated clinical variables (age at death and patient's disease onset) with the presence of 0, 1 or 2 sCJD risk alleles. In order to assess rs6951643 influence on age at onset and disease duration we performed a time to event analysis using Cox regression adjusting by PRNPM129V genotype and other relevant factors like sex and country of origin.
Supporting Information

Acknowledgments
The author thanks Simon Mead for sharing his data in order to select SNPs for replication. We also thank the Wellcome Trust Case Control Consortium (WTCCC) for the use of their control data. UK. The national UK CJD Surveillance Unit is grateful to clinicians, patients, and family members throughout the UK for a remarkable level of cooperation with the CJD surveillance program. France We thank all members of the French National Surveillance Network for Creutzfeldt-Jakob disease and all physicians for case notification. Italy We are very grateful to neurologists and neuropathologists, patients, and family members throughout Italy for their collaboration. The Netherlands We thank Pascal Arp, Mila Jhamai, Marijn Verkerk, Lizbeth Herrera, and Marjolein Peters for their help in creating the GWAS database, and Karol Estrada and Maksim V. Struchalin for their support in creation and analysis of imputed data. The authors are grateful to the study participants, the staff from the Rotterdam Study, and the participating general practitioners and pharmacists. Spain. The authors thank Biobanco Valdecilla and IDIVAL for its support throughout this project thanks to J.M. Polo for their critical advice, thanks to Beatriz Sobrino and María Torres (CEGEN, Spanish National Genotying Center) for all their work in the sample genotyping, and Weyma Notel for her help in editing the manuscript. Australia. The ANCJDR thanks all patients, families, clinicians and allied health personnel for their support and cooperation to facilitate this study.