Identification of Risk and Age-at-Onset Genes on Chromosome 1p in Parkinson Disease

Patients and Families

Affected individuals and family members were collected by the Morris K. Udall Parkinson Disease Research Center of Excellence (PDRCE) located within the Duke Center for Human Genetics (DCHG) and by the 13 centers of the Parkinson Disease Genetics Collaboration (Scott et al. 2001). A standard clinical evaluation involves neurological examination, including evaluation with the Unified Parkinson’s Disease Rating Scale (Goetz et al. 2003). A rigorous clinical assessment was performed by all participating clinicians to provide a clear diagnosis of PD and to exclude any individuals who displayed atypical features of parkinsonism (Hubble et al. 1999; Scott et al. 2001). Individuals categorized as “affected” showed at least two of the cardinal signs of PD (resting tremor, bradykinesia, and rigidity). AAO for affected individuals was defined as the age at which an affected individual first noticed one of the cardinal signs of PD. Participants categorized as “unaffected” demonstrated no signs of the disease, and participants categorized as “unclear” showed only one cardinal sign and/or atypical features. All participants signed informed consent forms prior to blood and data collection. Institutional review boards at each participating center approved the study protocols and consent forms.

The data set consists of multiplex (n=267) and singleton (n=361) white families. We defined singleton and multiplex families on the basis of the total number of parent-child triads and discordant sib pairs (DSPs) in the family that can contribute to the association test. Singleton families have only one group (either a triad or a DSP) contributing to the association test—that is, they have one affected individual—and either the parent (affected or unaffected) or the unaffected sibling is sampled in addition to the affected individual. Multiplex families have at least two groups (triads or DSPs) contributing to the association test—that is, they have at least two affected siblings sampled in the family. Families with Parkin mutation carriers were excluded from this study. The multiplex data set includes 609 affected individuals (average AAO ± SD = 61.0 ± 11.6 years; range 14–90 years; 58.8% males) and 666 unaffected individuals (42.8% males). The singleton families include 391 affected individuals (average AAO ± SD = 55.5 ± 13.0 years; range 15–85 years; 69% males) and 356 unaffected individuals (42.7% males).

DNA Extraction and Genotyping

DNA samples were prepared and stored by the DCHG DNA bank core. Genomic DNA was extracted from whole blood by use of the PureGene system (Gentra Systems Autopure LS). A total of 284 SNPs (table 1) were genotyped with the use of Applied Biosystems (ABI) Assays-on-Demand (AoDs) or Assays-by-Design (AbDs) or with the use of primers and probes designed using ABI Primer Express 2.0 software. The SNPs were chosen first on the basis of their location (e.g., SNPs were initially chosen with an average distance of 100 kb between them) and then on the basis of frequency to capture a wide range of frequencies among all selected SNPs. The TaqMan allelic discrimination assay was used to genotype all SNPs. PCR amplification was performed in 5-μl reactions (comprising 2.6 ng of dried DNA, 1× TaqMan universal PCR master mix [ABI], and 1× genotyping mix for AoDs and AbDs or, for self-designed assays, 900 nM of each primer and 200 nM of each probe). PCR was performed using the GeneAmp PCR system 9700 thermocycler (ABI), with a 40-cycle program (1 cycle at 95°C for 10 min, followed by 40 cycles at 95°C for 15 s and at T_m for 1 min, where T_m is 60°C for AoDs and AbDs and ranges from 58°C to 64°C for self-designed assays). The fluorescence generated during the PCR amplification was read using the ABI Prism 7900HT sequence detection system and was analyzed with SDS software (ABI).

Stringent quality-control measures were taken to ensure data consistency. Internal controls consisted of 24 duplicated individuals per 384-well plate. In addition, two samples from CEPH were plated eight times per plate to assure plate-to-plate consistency. All genotypers were blinded to these internal controls. Quality-control samples were compared at the DCHG Data Coordinating Center. Data were stored and managed by the PEDIGENE system (Haynes et al. 1995). To pass quality control, genotyping plates must have retained a 100% match for quality-control samples and must have at least 95% overall efficiency.

Candidate Genes Derived from the Genomic Convergence Approach

Two independent gene expression studies on human midbrain tissues from patients with PD and normal controls, by use of microarray and serial analysis of gene expression (SAGE) technologies, were conducted as part of current projects of the Duke PDRCE (Hauser et al. 2003; Noureddine et al. 2005). By combining these two studies, we found six genes that were significantly differentially expressed between patients with PD and controls and that mapped to the chromosome 1p AAO linkage region (table 2). In the present study, we tested SNPs in these six genes for association with risk and AAO in PD.

Iterative Association Mapping

We developed a second approach, iterative association mapping, to identify candidate genes in a linkage region. The overall concept is to reduce the number of SNPs genotyped while maximizing the chance of discovering a significant association. SNPs are first chosen at 100-kb intervals and then tested for association with traits of interest, which, in this case, are risk and AAO in PD. If no significant association is detected, the marker-to-marker distance is decreased by one-half each time (to 50 kb, 25 kb, etc.), until a significant association result is found. When a significant association is detected, additional SNPs in the surrounding region are then tested on the basis of known linkage disequilibrium (LD) patterns or physical iteration in the surrounding region of the associated SNP if no previous LD patterns are available.

Statistical Analyses

All SNPs were tested for Hardy-Weinberg equilibrium (HWE) and LD in the affected group (comprising one affected individual from each family) and the unaffected group (comprising one unaffected individual from each family). An exact test implemented in the Genetic Data Analysis program was used to test HWE, for which 3,200 permutations were performed to estimate the empirical P value for each marker (Zaykin et al. 1995). The Graphical Overview of Linkage Disequilibrium (GOLD) package was used to calculate LD (as measured by the Pearson correlation coefficient, r², and Lewontin’s standardized disequilibrium coefficient, D′) between pairs of SNPs (Abecasis and Cookson 2000). Both r² and D′ range from 0 (no LD) to 1 (perfect LD). However, there is no clear definition to use for interpretation of intermediate LD values. Here, we chose an arbitrary cutoff and considered two markers to be in strong LD if r²>0.60 or D^′>0.90.

AAO was treated as a quantitative trait. We used both the orthogonal model (OM) (Abecasis et al. 2000) and the Monks-Kaplan method (MKM) (Monks and Kaplan 2000) implemented in the Linkage Disequilibrium Analyses for Quantitative and Discrete Traits (QTDT) program (see Web Resources) to test the association between markers and AAO. The MKM not only provides an association signal but also detects the direction of association—that is, positive association for allele A is declared when the majority of allele A carriers have an AAO higher than the average AAO. In addition to obtaining nominal P values, we also performed 10,000 permutation tests to obtain an empirical P value for each marker on the basis of the MKM. The global significance level was derived from permutation tests.

We performed haplotype analysis for genes with significant markers. Prior to the haplotype analysis, we identified tagging SNPs (tagSNPs) for each gene by using the ldSelect program (see Web Resources) (Carlson et al. 2004). The ldSelect program generates groups of markers in LD on the basis of a given threshold for r². These groups are referred to as “LD bins.” A tagSNP is then selected from each LD bin. To perform the haplotype association analysis for AAO on the tagSNPs, we first used the FBAT –o option (Laird et al. 2000) to estimate the optimal offset of the AAO for each tagSNP. We then performed the HBAT –e option (Horvath et al. 2004) on the adjusted AAO data (subtracting from AAO the average optimal offset estimate), to test the association between haplotypes and AAO. When the number of tagSNPs is large, the computational time is substantial and the haplotype frequencies tend to be small, which makes interpretation difficult, even if significant P values are found. Therefore, we limited our haplotype computations to five tagSNPs. For genes with more than five tagSNPs, we analyzed all possible combinations of five tagSNPs.

The pedigree disequilibrium test (PDT) (Martin et al. 2000, 2003) was used to determine the association between markers and PD risk. Two PDT statistics were used: the PDT-sum statistic, for allelic effects, and the GenoPDT statistic, for genotypic effects. We also performed haplotype analysis on the risk genes detected by the PDT. The approach of selecting tagSNPs was as described above. We used the HBAT –e option to test the haplotype association between a set of tagSNPs and PD.

Several criteria were used to determine the final levels of significance in the presence of multiple comparisons. First, a significance level of P⩽.05 was used for evaluating the initial set of markers with 100-kb spacing. Second, a cluster approach (described below) was used to generate a significance level for further iterations. This approach requires that two or more markers, which have an r² correlation <0.6, be significant within a cluster of SNPs. Finally, at least one marker in the candidate gene or region needs to meet the global significance level derived from the permutation test.

Assume a total of N markers with low LD (r²<0.6) across the region of interest and x markers located in each cluster, which leads to y clusters (y=N/x). We hypothesized that a cluster would be significant only if two markers within the cluster are significant. We can formulate the probability (α_c) that 1 of the y clusters is significant as a function of the probability that a marker is significant (where α is the significance level of a marker):

By restricting the significance level of a cluster to be α_c, we can compute the probability that a marker is significant. In other words, the probability that two markers within a cluster are significant at the level of α will result in the probability α_c that one cluster is significant. Clearly, α decreases when the number of significant markers within a cluster decreases or when α_c, the significance level of a cluster, decreases. The calculation of the global significance level is described above.

The multiplex families used in this study include 167 families that were studied in the previous AAO linkage study (hereafter called “the linkage data set”) (Li et al. 2002). We performed SOLAR (Almasy and Blangero 1998) PEDLOD analysis with our previously found chromosome 1 peak marker (D1S2134), to obtain family-specific LOD scores for the 167 families. We then stratified the linkage data set into positive- and negative-linkage subsets on the basis of the family-specific LOD scores. The genes significantly associated with AAO in the overall data set were also tested for association with AAO by use of the MKM in the positive- and negative-linkage subsets. We did not use the OM, because it requires a normal distribution for the quantitative trait of interest, which is a problem for these small stratified data sets.

mRNA Analysis of USP24

Total RNA was isolated from human midbrain tissue and was reverse transcribed using poly-dT primers to generate a cDNA library. Primers used to amplify fragments of the USP24 transcript were designed using Primer3 (see Web Resources). (The sequences are available on request.) We generated several PCR products of the expected size from the cDNA library and sequenced them. The exon-intron structure of the complete USP24 transcript was deduced from genomic alignment of the overlapping RT-PCR fragments.

Identification of the Linkage Subsets of Families

The SOLAR PEDLOD analysis of D1S2134 identified 83 families with positive LOD scores (i.e., with positive linkage) and 84 with negative LOD scores (i.e., with negative linkage) from the linkage data set (Li et al. 2002). Throughout the present study, we performed association analyses with the overall PD data set as well as with these two stratified linkage subsets.

Genomic Convergence

We identified two differentially expressed genes from a previous microarray study (Hauser et al. 2005) and four from a SAGE study (Noureddine et al. 2005) that mapped to our chromosome 1p AAO linkage region (table 2). We analyzed SNPs (table 1) in each of these six genes, using the PD multiplex data set. The LD patterns of these six genes (pairwise r² values) are shown in tables 3–8.

The exclusion of a gene as a candidate from an association study is not always straightforward. The degree of confidence with which one excludes a gene from association is based on the depth of the search. One measure is at the level of LD defined by the current HapMap data set. Because we began genotyping our data set prior to the availability of the HapMap data set, and because we genotyped as many SNPs with as wide a variety of frequencies as possible from what was available in public (NCBI) and private (ABI) databases, some of our markers are not in the HapMap data set. To evaluate whether we have sufficiently covered each gene, we compared our SNP coverage of each gene with that available in the current HapMap data. The number of LD bins identified on the basis of HapMap SNPs with an MAF >10% is as follows: one LD bin each for ATP6V0B, UQCRH, and C1orf8; two for TTC4; three for RNF11; and 12 for PPAP2B. Overall, our SNPs included the HapMap tagSNPs in all genes except RNF11 and PPAP2B; we missed one HapMap tagSNP in RNF11 and covered only two HapMap tagSNPs (of seven genotyped SNPs) in PPAP2B.

None of these genes show significant association with PD risk, and only SNP 193 in C1orf8 showed significant association with AAO in PD (fig. 2). The association of SNP 193 was not verified in the positive-linkage subset.

ELAVL4

The embryonic-lethal, abnormal vision, Drosophila-like 4 gene (ELAVL4) encodes for a neuron-specific RNA-binding protein. This gene was studied as a biological candidate gene through an ongoing project of the Duke PDRCE (Antic and Keene 1997). Two polymorphisms (SNPs 136 and 143) were previously found to be significantly associated with AAO in PD (Noureddine et al. 2005). However, these markers were not found to reach significant P values in the positive-linkage subset in the present study.

Iterative Association Mapping and LD

The initial association map consisted of 200 SNPs (1 SNP genotyped, on average, every 100 kb) in the 1-LOD score region (40.4–59.2 Mb on NCBI build 34) (fig. 1). With additional genotyping in the regions of interest, the average SNP density in our final association map was 1 marker every 66 kb, with a total of 284 SNPs genotyped. The MAFs of the SNPs varied from 0.03 to 0.50 (median and average of 0.29). All but 20 (7%) of the SNPs were in HWE in both the affected and unaffected samples at a P=.05 level (table 1). The genotype distributions of these 20 SNPs were reexamined by a technician in the laboratory and were tested for HWE again. The results remained the same. Given a 5% random chance of obtaining markers not in HWE, the 7% frequency detected in our project is within a reasonable range. Furthermore, it is important to note that the MKM and the PDT do not require HWE.

Figure 3 depicts the overall LD structure in this region (for SNPs that did not deviate from HWE, since LD calculations assume equilibrium) in the unaffected group. A similar LD pattern was observed in the affected group. LD is mostly restricted to intragenic areas, with no extensive LD for long stretches of DNA, or across distant loci for the majority of polymorphisms. Only SNPs with a low MAF (recent SNPs) show high values of D′ with most neighboring SNPs.

To obtain a P value for the cluster analysis, 210 markers were identified whose r² was <0.6 for LD. Using these 210 markers and assuming 7 markers within each cluster, we derived a significance level of .01 for each marker. In addition, we obtained a global significance level of .001. Figures 2 and 4 depict the nominal P values of all 284 markers in −log₁₀ scale. Among the first 200 SNPs studied (from the 100-kb map), evidence for association with AAO was found by either the OM or the MKM in the genes for the eukaryotic translation initiation factor EIF2B3 (for SNP 63, P=.009 by the OM and .0004 by the MKM), the testis-specific protein kinase 2 (TESK2 [for SNP 76, P=.008 by the MKM]), hypothetical protein FLJ14442 (for SNP 117, P=.01 by the MKM), and the ubiquitin-specific protease 24 (USP24 [for SNP 220, P=.004 by the OM]) (figs. 2 and 5). All these markers have empirical P values by the permutation test that are slightly lower than the nominal P values. For example, the empirical P value for SNP 63 in EIF2B3 was .0002. Evidence of association with risk for PD by use of the multiplex data set was found only in the HIV type 1 enhancer-binding protein 3 gene (HIVEP3), for SNPs 13 (P=.008) and 19 (P=.004) (fig. 4). We proceeded to increase the SNP density in these genes.

TESK2 and FLJ14442

Additional SNPs (SNPs 72, 74, and 75 in TESK2 and SNPs 116, 118, 120, 122, and 124 in FLJ14442) were genotyped, to a final average density of 1 marker per 29 kb for TESK2 and 1 marker per 51 kb for FLJ14442. Although we detected two sets of cluster markers for AAO association, no markers were significant after correction for multiple testing, nor did they show evidence of association in the positive-linkage subset.

EIF2B3

Ten additional SNPs (SNPs 57–62 and 64–67) were genotyped in the EIF2B3 gene (size 136 kb), leading to a final average density of 1 marker per 12 kb (fig. 2). Interestingly, several markers that were close to significance in the overall data set became significantly associated with AAO in the positive-linkage subset (table 9), despite the fact that the subset is only one-third (83 families) of the total sample. Therefore, at least two clusters of markers in low LD (r²<0.6) (SNPs 59–61 and 62–64) are strongly associated with AAO in this gene. More interestingly, SNPs 62–64 are still significant after correction for multiple testing (P<.001).

Five tagSNPs (SNPs 59, 60, and 64–66) were found in EIF2B3. Haplotype analysis with these five tagSNPs and the use of the overall PD data set produced two haplotypes significantly associated with AAO: C-C-G-T-G (haplotype frequency 17.2%; P=.002) and A-C-A-T-G (haplotype frequency 15.2%; P=.002) (table 10). These two haplotypes showed P values comparable to what we detected for SNP 64 alone (P=.01 by the OM and .0001 by the MKM).

USP24 and AK127075

In total, we genotyped 14 SNPs (SNPs 218–231) with ∼17-kb spacing in the region from USP24 to the cDNA FLJ45132 clone BRAWH3037979 (GenBank accession number AK127075), a region in which seven SNPs (SNPs 220–222, 224, 227, 230, and 231) are significantly associated with AAO (P<.01) (fig. 2). The most significant marker was SNP 227, with P=.0006 by the OM and .007 by the MKM.

In silico, several lines of evidence suggested that the annotated USP24 gene in NCBI build 34 (as defined by the mRNA for KIAA1057 protein [GenBank accession number AB028980]) may actually be a truncated version of the full-length USP24 transcript. The 5′ end of the AB028980 transcript (exons 1–11) matches the 3′ end of the AK127075 mRNA (exons 25–35), and the human THC1877380 transcript from the TIGR Human Gene Index overlaps both genes. GenScan predicts the existence of the NT_032977.390 mRNA (composed of the AB028980 and AK127075 mRNAs and 12 additional exons at the 5′ end), and there is a cluster of human overlapping spliced ESTs (e.g., GenBank accession numbers BM458550, AW853346, and CD687922) that support the existence of a longer USP24 transcript. Furthermore, the mouse AK045043 significantly overlaps with this cluster of ESTs but has two additional distant exons at the 5′ end. The putative first exon is supported by the FirstEF program prediction, contains an ATG start codon with sequences conforming to a Kozak consensus ([A/G]CC ATG G), has a nearby CpG island, and is close to predicted promoter sequences; all of which strongly reinforce the idea that it encodes the first exon of the larger USP24 ORF. This gene produces a predicted mRNA of ∼8 kb.

To evaluate the existence of this larger USP24 transcript, termed “USP24_L,” we used strategically positioned primers to amplify overlapping transcript fragments from a human midbrain cDNA library. We obtained RT-PCR products of the expected sizes, and direct sequencing of these products confirmed the existence of USP24_L. Using the BLAST tool implemented at the University of California–Santa Cruz Web site, we aligned the experimentally amplified composite cDNA with the genomic sequence. The sequence of our USP24_L transcript carried more exons than the GenScan NT_032977.390 and XM_371254 predictions, some of which are supported by human or mouse ESTs. All splice junctions followed the canonical AG/GT rule. The composite cDNA is predicted to encode a protein of 2,590 aa (fig. 5), distributed over 69 exons and spanning 146 kb of genomic sequence (54904635–55050704 bp in chromosome 1). The LD block observed from SNP 216 through SNP 231 (fig. 3), which encompasses USP24_L and flanking regulatory sequences only, also supports the predicted size of USP24_L.

Since the SNPs significantly associated with AAO in this region completely span USP24_L, and since strong LD exists throughout USP24_L but not with neighboring genes, we concluded that the association originates from USP24_L itself. Three LD bins were found in this region on the basis of 14 SNPs genotyped in this study (SNPs 218–231). The seven SNPs significantly associated with AAO were, in fact, originating from two LD bins. The first LD bin is formed by SNPs 220, 221, 224, and 230 (maximum P=.007), and the second is formed by SNPs 222, 227, and 231 (maximum P=.003), which implies that there are two independent polymorphisms in USP24_L that have significant effect on AAO. Although none of the SNPs in USP24_L were significantly associated in either the positive- or negative-linkage subsets by the MKM, SNPs 221, 224, and 230 were close to significant (.05<P<.06) in the positive-linkage subset (table 9).

Three tagSNPs (SNPs 218, 219, and 227) were identified in USP24. Two haplotypes, C-T-T (frequency 62.6%; P=.003) and C-T-C (frequency 19.9%; P=.026), were found to be significantly associated with AAO (table 10). Overall, these haplotypes in USP24 did not provide any more information on the association with AAO than did SNP 227 alone.

HIVEP3

A total of nine markers in this gene were genotyped, with a final average density of 1 marker every 45 kb. The new SNPs failed to reveal any further significant association with risk for development of PD (fig. 4). However, SNP 12 was close to significant in both the allelic (P=.058) and the genotypic (P=.057) association tests, and SNP 18 (P=.059) was close to significant in the PDT, since it is in relatively high LD with SNP 19 (r²=0.75 in the unaffected group). To test for association of SNPs 13 and 19 in a second independent data set, we genotyped these two markers in the PD singleton data set. We did not find evidence of association of these SNPs in the singleton data set alone. However, both markers showed stronger significant association in the combined multiplex and singleton data set (P=.006 for SNP 13; P=.002 for SNP 19) than in the multiplex data set. Clearly, some singleton families also contribute to the association of these two markers.

We identified eight tagSNPs (SNPs 13–17 and 19–21) in HIVEP3. Haplotype analyses based on five tagSNPs revealed the best results by use of tagSNPs 13, 15, 17, 19, and 21, in which a rare A-G-T-G-C haplotype (frequency 2.1%) was significantly associated with risk for PD (P=.003) (table 10). HIVEP3 is a relatively large gene (408 kb), and very low levels of LD were observed among the SNPs genotyped (fig. 3). The lack of LD between SNPs 13 and 19 (r²=0; D^′=0.02) provides two independent lines of evidence for the involvement of this gene in controlling risk for development of PD.