Genome-wide Analysis of Copy Number Variation in Latin American Parkinson’s Disease Patients

Introduction

Parkinson’s disease (PD) is the second most common neurodegenerative disorder, and the fastest growing cause of disability due to a neurological disorder in the world ^1,2. PD is a multifactorial syndrome that is thought to be caused by the complex interaction of genetics, environmental factors, and aging ³.

Evidence for the genetic basis of PD has increased substantially over the past decades ^4–6. The first causal gene for PD, SNCA, was discovered in 1997 ⁷, and its protein product (α-synuclein) was further shown to be a major component of Lewy bodies, the pathological hallmark of PD. Dominant pathogenic single nucleotide variants (SNVs) in SNCA ^8–11, as well as copy number variants (CNVs), such as duplication or triplication of the entire gene with a clear dose effect, have been reported ^12–14. The discovery of SNCA was followed by that of PARK2 ¹⁵, where both pathogenic SNVs and CNVs are associated usually with autosomal recessive, early-onset form of the disease ¹⁶. Almost exclusively, genetic discoveries in PD have focused on SNVs, and studies on CNVs have been infrequent ^17–19. CNVs in PARK2, SNCA, PINK1, DJ1, and ATP13A2 (from more to less frequent) have been reported using a candidate gene approach ^20–22, while no CNVs have been shown for LRRK2. To date, only two studies have investigated the role of CNVs in PD at a genome-wide level, including exclusively European and Ashkenazi Jewish individuals ^17,19, with a sample size of 1672 and 432, respectively.

PD is a global disease affecting all ethnicities. Unfortunately, the majority of studies do not include individuals of non-European ancestry, creating a large gap in knowledge. This is especially true for Hispanics/Latinos (Box 1). Despite the fact that they are the largest and fastest growing ethnic minority in the US²³, Hispanics/Latinos are critically underrepresented in most genetic studies ²⁴. This is probably due to their complex admixed ancestry with influences primarily from European, Amerindian and African populations. In the US, the incidence and prevalence rates of PD among Hispanics are at least as high, if not higher than in non-Hispanic Whites, while the rates are lower for Asians and Blacks^25,26. Yet, little is known about the genetics of PD in Hispanics/Latinos, especially the frequency and characteristics of CNVs. No genome-wide studies in this population have been performed to date.

To address the lack of diversity in PD genetic studies and to understand the genetic architecture of PD in Latinos, we created the Latin American Research Consortium on the Genetics of PD (LARGE-PD) ²⁷. For this study, we used genome-wide genotypes of 1,497 individuals from LARGE-PD. The aim of this study was to elucidate genomic structural changes, as well as assess the CNV burden in this cohort of Latino PD patients and controls.

Methods

As part of our ongoing collaborative effort within LARGE-PD ²⁷, we examined data from a total of 1,497 individuals (807 PD patients and 690 controls) recruited from nine different sites across the following five different countries: Peru (n = 721), Colombia (n = 351), Brazil (n = 227), Uruguay (n = 191), and Chile (n = 13). All patients were evaluated by a movement disorder specialist at each of the sites and met the UK PD Society Brain Bank clinical diagnostic criteria²⁸. Controls were selected from ancestry matched individuals that did not have symptoms compatible with neurodegenerative disorders. All PD patients and controls provided signed informed consent according to the local ethical requirements of each site. All individuals were genotyped on Illumina’s Multi-Ethnic Global Array (MEGA) (Illumina, San Diego, CA, USA). A total of 1,779,819 markers were available before quality control (QC).

We performed an initial round of QC using PLINK 1.90 ²⁹, based on single nucleotide polymorphism (SNP) genotype data for all samples and following established protocols described in Niestroj et al ³⁰. Samples with a call rate < 0.96 or a discordant sex status were excluded. We filtered autosomal SNPs for low genotyping rate (> 0.98), case-control difference in minor-allele frequency (> 0.05), and deviation from Hardy-Weinberg equilibrium (HWE, P-value ≤ 0.001) before pruning SNPs for linkage disequilibrium (--indep-pairwise 200 100 0.2) using PLINK ²⁹. 1000 Genomes population ³¹ was used as a reference for visual clustering of the Principal Component Analysis (PCA) to assess for population stratification.

For CNV calling, we focused only on autosomal CNVs due to the higher quality of CNV calls from non-sex chromosomes. A custom population B-allele frequency (BAF) file was generated as a reference before calling CNVs. Then, we created GC wave-adjusted Log R Ratio (LRR) intensity files for all samples and employed PennCNV ³² software to detect CNVs in our dataset.

We assessed cryptic relatedness using KING ³³ software, and excluded individuals who were closely related (up to second degree) to another participant in our cohort by using the unrelated algorithm in KING. We performed an intensity-based QC to remove samples with low-quality data as previously described in Huang et al. ³⁴ Following this step, all samples had a LRR standard deviation of < 0.27, absolute value of waviness factor < 0.03, and a BAF drift < 0.0014.

Called CNVs were removed from the dataset if they spanned < 20 markers, were < 20 Kb in length, and had a SNP density < 0.0001 (amount of markers/length of CNV). SNP density was not considered for CNVs spanning ≥ 20 markers and ≥ 1Mb in length, as larger CNVs are not likely to be artifacts. To ensure that only high-quality CNVs passed our filters, we implemented a quality score calculation for each CNV following the methods of Macé et al. ³⁵, in which various CNV metrics are combined to estimate the probability of a called CNV to be a consensus call. Quality scores ranged from 0 (lowest) to 1 (highest) for duplications and similarly from 0 to −1 for deletions, and CNVs with quality scores between −0.5 to 0.5 were filtered out. A subset of final QC-passed CNVs were also inspected visually by five different investigators with expertise in the interpretation of BAF and LRR plots. CNVs were annotated for gene content using Ensembl ³⁶ including gene name and the corresponding exonic coordinates in hg19 assembly using bedtools 2.27.0 ³⁷.

We calculated CNV burden for PD using different categories to evaluate the relative contribution on PD risk: (1) the carrier status of overall CNV burden, including CNVs in non-genic regions (2) the carrier status of CNVs intersecting ‘any gene’ but none of the PD genes, (3) the carrier status of CNVs intersecting a list of “known PD genes”, and (4) the carrier status of large CNVs (≥ 1Mb in length). P values were adjusted with the false discovery rate (FDR) method to correct for multiple testing. For the overall CNV burden category, deletions and duplications were also analyzed separately. We selected 19 genes for the “known PD genes” category that were grouped as follows: six are well-established causal genes for PD (LRRK2, PARK2, PARK7, PINK1, SNCA, and VPS35), one is a susceptibility factor for PD with a large effect size (GBA), and 12 either result in a parkinsonian syndrome that sometimes overlaps PDptio or are putative causal genes for PD that have not be adequately validated (ATP13A2, DNAJC6, DNAJC13, EIF4G1, FBXO7, GIGYF2, HTRA2, PLA2G6, RAB39B, SYNJ1, TMEM230, and VPS13C) ^6,38,39.

To assess for the difference in CNV burden between PD patients and controls, we fitted a logistic regression model using the “glm” function of the stats package ⁴⁰ in R 3.6.0 ⁴¹. Cox-proportional hazards regression analyses and Kaplan-Meier curves were calculated using the survival package ⁴². For all burden analyses, odds ratios (OR), 95% confidence intervals (CIs), and significance were calculated. ORs were calculated by the exponential of the logistic regression coefficient. For Cox-proportional hazards regression, hazard ratios (HR) were calculated to allow for censored observations. Potential confounding variables were used as covariates and included age, sex, and the first five ancestry principal components for all regression models.

Results

We had available data from a total of 1,497 individuals in LARGE-PD. We excluded 39 individuals due to relatedness, and 79 due to failing our intensity-based QC steps. Thus after QC, our final cohort included 1,379 individuals (747 PD patients and 632 controls) from Peru (n = 677), Colombia (n = 320), Brazil (n = 192), Uruguay (n = 177), and Chile (n = 13). There were more males in PD patients compared to controls (53.2% vs 33.1%, P < 0.001). Sample demographics are shown in Table 1. To visualize the ethnic composition of our cohort, we performed PCA using 1000 Genome populations ²⁸ as a reference. Our samples overlapped strongly with the projection of Admixed Amerindian samples (AMR) (Supp. Fig. 1).

The initial number of CNV calls was 249,101 including 176,462 deletions and 72,639 duplications. After all QC steps, including filtering by consensus quality score as described in the methods, the final number of high-quality CNVs was 8,412, including 5,155 CNVs in PD patients and 3,257 CNVs in controls. CNV analysis showed 1,274 of the samples (92.4%) carrying at least one QC-passed CNV. The length of the CNVs in the overall cohort ranged from 20 kilobases (Kb) to 3.4 megabases (Mb), with a median size of 52.4 Kb. CNV characteristics are shown in Supp. Table 1.

We applied logistic regression to compare the CNV burden in PD patients and controls on all categories defined earlier, adjusting P values for multiple testing (see methods for details). We found no significant difference in overall CNV burden (OR: 1.19 [0.78 - 1.8], P = 0.64), CNVs in any gene (OR: 1.07 [0.81 - 1.4], P = 0.77), and large CNVs that are ≥ 1Mb in length (OR: 1.46 [0.82 - 2.65], P = 0.4) (Fig. 1). Interestingly, nine PD patients and six controls carried a duplication that was ≥ 1Mb on chromosome 11 (P = 0.8), which did not overlap with any known gene region. We also calculated overall CNV burden independently in duplications (OR: 1.36 [1.05 - 1.77], P = 0.06) and deletions (OR: 0.98 [0.74 - 1.29], P = 0.89), and neither were statistically significant.

• Download figure
• Open in new tab

• Download figure
• Open in new tab

Figure 1:

(A) Forest plot showing the CNV burden compared between PD patients and controls. Odds ratios (ORs) and P values were calculated using a logistic regression for CNVs corrected with age, sex, and first 5 components of PCA. P values were adjusted with FDR for multiple testing. ORs > 1 indicates an increased risk for PD per unit of CNV burden. (B) Table showing number of CNV carriers in any of the 19 known PD genes. (C) Visualization of CNVs on PARK2.

We then explored CNVs on genomic regions that were previously associated with typical PD and other parkinsonian phenotypes, and we found that PD patients were significantly enriched with CNVs overlapping these genes (OR: 3.97 [1.69 - 10.5], P = 0.018) (Fig. 1). This finding was largely driven by CNVs on PARK2 in 20 patients, followed by two patients with a CNV on SNCA, compared to six controls carrying a CNV on PARK2 and none on SNCA. In addition one control had a CNV on PLA2G6 (Fig. 1).

To assess whether PD patients carrying a PARK2 CNV in our cohort had an earlier age at onset (AAO), we performed a Mann-Whitney test and compared patients that carry a CNV overlapping PARK2 to those carrying a CNV overlapping any other gene. The median AAO for patients with a CNV overlapping PARK2 was 31 years old, while that for patients with a CNV overlapping any other gene was 57 years old (P = 7.46 x 10^-7). To further investigate the CNV burden in early-onset PD (EOPD) patients, we performed a subset analysis in which cases with an AAO < 50 years old were compared to controls. Again, we observed a significant enrichment in CNVs on known PD genes in patients with EOPD (OR: 4.91 [1.92 - 13.68], P = 0.006). This result was also driven by CNVs on PARK2.

Kaplan-Meier estimates of AAO showed that individuals carrying a CNV on a known PD gene had significantly earlier onset of symptoms compared to individuals with other or no CNVs (log-rank test, P < 0.001) (Fig. 2). Using a Cox proportional-hazards regression analysis with age, sex, and the first five ancestry principal components included as covariates, we found that the effect of carrying a CNV on a known PD gene on the hazard of AAO was highly significant (HR: 2.42 [1.57 - 3.71], P = 5.70 x 10^-5). We also assessed AAO in PD patients only, comparing CNV carriers on a known PD gene to PD patients with other or no CNVs, and found that having at least one CNV results in earlier onset of symptoms. (HR: 1.92 [1.22 - 3.02], P < 0.001) (Supp. Fig. 2).

• Download figure
• Open in new tab

Figure 2:

Kaplan-Meier estimates of individuals (PD patients and controls) carrying a CNV on a PD gene and individuals with other or no CNVs. Controls are censored observations since it is only known that they did not develop PD up to the age of their last visit. Probability: probability of not having symptoms associated with PD. Age at Visit or Age at Onset: time to onset of PD symptoms for cases and time to last visit for controls. Highlight around the curves shows 95% confidence intervals.

Discussion

Here, we present a genome-wide characterization of CNVs in a cohort of Latino PD patients and controls from LARGE-PD ²⁷. We analyzed genotypes of 1,497 individuals on the same platform and analyzed all samples with the same CNV calling and quality control pipeline. We used ancestry matched controls for the interpretation of CNVs detected in PD patients. This is particularly important considering that the data for Latino population frequency of CNVs is limited, especially in neurologically healthy adults ^43,44. We assessed the CNV burden for different categories and observed an increased burden of CNVs overlapping known PD genes in PD patients. We identified 22 patients that carried CNVs overlapping two established PD genes (PARK2 & SNCA), and found that 14 of these patients had a disease AAO < 50 years. The median AAO for patients with a CNV overlapping PARK2 was almost 20 years earlier than that of other patients, in agreement with the literature ^45–47.

PARK2 mutations are the most common genetic cause for EOPD ^5,15,48, but an important caveat is that this information is mostly derived from studies in populations of European ancestry. The frequency of all PARK2 pathogenic variants (CNVs and SNVs) in EOPD patients with European ancestry ranges from 49% in familial cases, to 15-18% in isolated patients, while the frequency of carrying a CNV in PARK2 in isolated EOPD patients is approximately 10%^45,49,50. Some studies suggest that alterations in PARK2 are more frequent in Hispanic populations. One study showed a 2.7 fold increase for carrying any PARK2 alteration and a 2.8 fold increase for carrying a heterozygous mutation in Hispanic EOPD patients (N = 77) compared to White non-Hispanics, and the frequency of PARK2 CNVs in Hispanic PD patients in this cohort was 6.4% ⁴⁸. In another study examining Mexican-mestizo EOPD patients (N = 63), the frequency of PARK2 CNVs was found to be 50%, and 18% of these patients were heterozygous ⁵¹. In our cohort, we had 250 patients with EOPD, and 5.6% (N = 14) of these patients carried a CNV on PARK2.

The role of homozygous and compound heterozygous variants, including CNVs on PARK2 is well known, especially in EOPD ^5,15,45. However, there is also increasing evidence that PARK2 heterozygosity is a risk factor for PD and is associated with a decreased AAO ^45,46. In our cohort, there was a significant association between the AAO and PARK2 carrier status. Still, the role of heterozygous PARK2 CNVs in altering PD susceptibility remains controversial ^16,52. In order to correctly characterize PD patients, an integrated SNV-CNV analysis is needed, given the importance of both allele types for comprehensive genetic diagnosis in PD ⁵³. Pankratz et al. showed that the frequency of carrying a single PARK2 CNV was higher in PD patients compared to controls, while it was similar for carrying a single point mutation ⁵⁴ Heterozygous PARK2 CNV carrier status may still play a role in the development of PD despite its recessive inheritance ^46,55, through a haploinsufficiency effect ⁵⁶.

Some limitations of our study are that we did not validate all CNVs with a different method or sequence all PARK2 CNV carriers. From our overall cohort (N = 1,379), a small portion were previously sequenced and/or CNV screened for other ongoing projects (see Supp. methods). Out of the 250 EOPD patients included in our analysis, 77 (30.8%) of them had been previously screened for CNVs with multiplex ligation-dependent probe amplification. This included seven of the 22 patients that were found in this current study carrying a CNV on PARK2 or SNCA. Results for all 77 samples matched, confirming negatives and all seven positive carriers, showing that our CNV analysis pipeline was accurate and detected 100% of the known CNV carriers (with no false positives) found with a different method. In addition, out of the 14 EOPD patients who carried a PARK2 CNV in our cohort, eight of them were also previously sequenced for PARK2. Based on data from the previous studies ^48,51, out of the 14 EOPD patients who carried a PARK2 CNV in our cohort, one would expect five to seven of them to be compound heterozygous for a PARK2 mutation of any kind. We found three of them to be compound heterozygous, carrying both a CNV and a pathogenic SNV in PARK2, and one of these patients had a novel acceptor splice site mutation, reported previously ⁵⁷. We also identified a fourth compound heterozygous sample through imputation (Supp. methods). Details for these four compound carriers are included in Supp. Table 2. This is very similar to the frequency reported in previous studies mentioned above^48,51.

The large genetic variation in Latinos due to admixture from several populations (mostly European, Amerindian and African) creates a challenge when analyzing this population. For this reason, we established a workflow with rigorous quality control. We also constructed a Latino reference file from scratch for CNV calls, as publicly available reference files were all based on Europeans. In this study, we analyzed all samples together in order to boost statistical power. However, separate calling of CNVs in subpopulations based on admixture analysis is likely to yield more refined results. Thus, larger sample sizes will be needed to make discoveries specific to subpopulations of Latinos. Admixture mapping to examine the chromosomal location of the CNVs could also provide more insights about the relationship between PD genetics and ethnicity^58,59.

To our knowledge, this is the first study that focuses on genome-wide CNVs in PD patients from Latin America. We believe that expanding the diversity of genetic studies for PD is necessary to understand the genetic profiles of these individuals and that our work will enrich current scientific knowledge about CNVs in this underrepresented population.

Glossary

AAO = age at onset; BAF = B-allele frequency; CI = confidence interval; CNV = copy number variant; EOPD = early-onset Parkinson’s disease; HR = hazards ratio; Kb = kilobases; LRR = log R ratio; Mb = megabases; OR = odds ratio; PCA = principal component analysis; PD = Parkinson’s disease; QC = quality control; SNP = single nucleotide polymorphism; SNV = single nucleotide variant.

Funding

This work was funded by a Stanley Fahn Junior Faculty Award from the Parkinson’s Foundation and supported by a research grant from the American Parkinson’s Disease Association, with resources and the use of facilities at the Veterans Affairs Puget Sound Health Care System.

Authors’ Roles

EIS and IFM designed the study. EIS, EPP, and LMN analyzed the data. IFM and DL supervised the study. EIS and IFM wrote the manuscript. All authors interpreted the data and revised the manuscript.

Supplementary

• Download figure
• Open in new tab

Supp. Figure 1:

(A) PCA of LARGE-PD individuals included in the study, color-coded by country of origin. (B) PCA of 1000 Genomes population for comparison, color-coded by different populations.

• Download figure
• Open in new tab

Supp. Figure 2:

Kaplan-Meier estimates of PD patients carrying a CNV on a known PD gene and patients with other or no CNVs. Probability: probability of not having symptoms associated with PD. Age at Onset: time to onset of PD symptoms. Highlight around the curves shows 95% confidence intervals.

Supp. Methods:

From our overall cohort (N = 1,379), 117 (8.4%) of the individuals were previously sequenced and/or CNV screened for other ongoing projects. Forty of the samples with self-reported family history were capture-sequenced for all coding regions using a custom neurodegenerative panel which included 21 PD related genes: ATP13A2, ATP1A3, DNAJC6, DNAJC13, EIF4G1, FBXO7, GCH1, GBA, GIGYF2, HTRA2, LRRK2, PARK2, PARK7, PINK1, PLA2G6, RAB39B, SLC6A3, SNCA, SNCB, SYNJ1, TAF1, TARDBP, TMEM230, VPS13C, and VPS35. Seventy-seven patients with EOPD were Sanger sequenced and CNV screened using multiplex ligation-dependent probe amplification for PARK2 and SNCA. We also imputed LARGE-PD dataset (807 PD patients and 690 controls) using the Trans-Omics for Precision Medicine (TOPMed) Project ⁶⁰ which identified a fourth compound heterozygous individual.