We enrolled 118 consecutive unrelated PD patients (68 males, 50 females) with 97 unrelated controls (28 males, 69 females), matched according to age and gender whenever possible. Thirty-five patients (16.28%) reported having a relative (≤ 2°) with any sort of movement disorder; of those, 21 (9.77%) had a formal PD diagnosis. All subjects resided and were originated from Costa Rica and were recruited at the Movement Disorders Unit of the Department of Neurology, Hospital San Juan de Dios, Caja Costarricense de Seguro Social. All patients fulfilled Gelb criteria for the clinical diagnosis of PD, while controls had no signs or personal history of any neurodegenerative disease and were mainly the spouses of the PD cases. We preferred using Gelb criteria over the United Kingdom Parkinson's Disease Society Brain Bank (UKPDSBB) as it provided different clinical diagnostic levels of certainty (possible and probable) and it has shown to have similar positive and negative predictive values, as well as sensitivity and global accuracy when compared to UKPDSBB (29). Albeit both diagnostic criteria sets have low specificity and are mainly focused on motor features, UKPDSBB criteria further err by challenging PD diagnosis in the presence of genetic risk factors (30). Our last patient was enrolled by 2011, which is 4 years earlier than when the Movement Disorder Society (MDS) task force proposed the new clinical diagnostic criteria for PD (MDS-PD criteria) (31); therefore, we were not able to use those for clinical diagnosis of patients enrolled in our study.

We gathered information concerning work and educational status as well as history of exposure to risk and protective factors of PD. We further obtained detailed information on PD history, comorbidities, and antiparkinsonian treatments. Additionally, motor disability of the patients was evaluated by means of the Unified Parkinson's Disease Rating Scale (UPDRS), Hoehn & Yahr (H&Y), and Schwab & England (S&E) scales. Cognitive status was assessed using the Montreal Cognitive Assessment (MoCA) test.

Molecular inversion probes were used to sequence coding and untranslated regions in familial PD and atypical parkinsonism-associated genes including GBA, SNCA, VPS35, LRRK2, GCH1, PRKN, PINK1, DJ-1, VPS13C, and ATP13A2 at McGill University with Illumina HiSeq 4000 as previously described (32). The full protocol can be found at https://github.com/gan-orlab/MIP_protocol. All sequences were aligned using Burrows-Wheeler Aligner (BWA) using the reference genome hg19 (33). Genome Alignment Tool Kit (GATK v3.8) was used to call variants and perform quality control and ANNOVAR was used to annotate each variant (34, 35). Exons 10 and 11 of GBA were sequenced using Sanger sequencing as previously described (36), and GBA variants in other exons were also validated using Sanger sequencing. We decided to focus on genes that are involved in typical PD, as our selected cohort is of typical PD (10, 37, 38).

All samples and variants were filtered based on standard quality control process as previously reported (39). In brief, variants were separated into common and rare by minor allele frequency (MAF) in the cohort. Rare variants (MAF < 0.01) with a minimum depth of coverage of > 30× were included in the analysis, along with common variants (MAF ≥ 0.01) with >15× coverage. We have established that for common variants, we get reliable reads at 15×; however, to get reliable reads for rare variants, we need > 30×; otherwise, there are many false positives (40). Variant calls with a genotype frequency of <25% of the reads or genotype quality of <30 were excluded. Samples and variants with more than 10% missingness were also excluded.

The atomic coordinates of the human Lrrk2 C-terminal domain structure (a.a. 1327–2527) were downloaded from the Protein Data Bank (ID 6VP6). The figure was generated using PyMol v.2.4.0.

We used Stata® (version 14) for the statistical analysis of sociodemographic and clinical variables. Normally distributed variables are reported as mean with its standard deviation (SD), whereas continuous but non-normally distributed variables are reported as median with the 25th and 75th percentile values (interquartile range, IQR). Normally distributed variables were compared with paired or unpaired t-tests, while non-normally distributed variables were compared with Mann–Whitney U-test or Wilcoxon match-paired signed-rank test. Frequencies were compared with χ² and Fisher's exact test. Tests were two-tailed, and significance was set at p < 0.05. We modeled through linear regression the association between demographic and clinical variables with the severity of the disease, as indexed by UPDRS and MoCA, as dependent variables in the models.

For genetic analysis, common and rare variants were analyzed separately. Association of common variants was tested using logistic regression adjusted for age and sex in PLINK v1.9. For rare variants' analysis, we examined the burden of rare variants in each gene using optimized sequence Kernel association test (SKAT-O) adjusted for age and sex (41). Rare variants were separated into different categories based on their potential pathogenicity to examine specific enrichment in different variant subgroups as described previously (40): (1) variants with Combined Annotation Dependent Depletion (CADD) score of ≥12.37 (representing the top 2% of potentially deleterious variants) (42); (2) regulatory variants predicted by ENCODE (43); (3) potentially functional variants including all non-synonymous variants, stop gain/loss variants, frameshift variants, and intronic splicing variants located within two base pairs of exon–intron junctions; (4) loss-of-function variants, which includes stop gain/loss, frameshift, and splicing variants; and (5) only non-synonymous variants. Bonferroni correction for multiple comparisons was applied as necessary.

This study was approved by the Ethics Committee of Hospital San Juan de Dios, Caja Costarricense de Seguro Social (CLOBI-HSJD #014-2015) and the University of Costa Rica (837-B5-304). Written informed consent was obtained from all participants.

At enrollment, PD probands had a mean age of 62.12 ± 13.51 years (range 25–86), and the mean age at onset was 54.62 ± 13.54 (range 16–83) years. Male PD patients comprised 57.63% of the sample. Despite the fact that a significantly larger proportion of the male PD patients reported current or previous jobs involving agricultural activities (19.40% male, 2.08% female; p = 0.01), the mean number of years of education of these men was significantly higher than women (10.74 ± 3.81 vs. 8.86 ± 4.01; p = 0.03). Table 1 details subjects' baseline characteristics along with the frequency of exposure to main risk and protective factors for PD. Most of the risk and protective factors were more prevalent in men. Tables 1, 2 detail the frequency of clinical manifestations as well as the standardized scale scores reported for PD cases. The most frequent initial symptoms included resting tremor (71.30%), rigidity (24.07%), and pain (10.19%). Most of the patients had asymmetric onset (94.12%) and a good response to levodopa (89.11%). Other frequently reported motor features comprised dystonia (46.08%), falls (39.22%), and dysphagia (36.27%). Common non-motor manifestations such as hyposmia, sleep disorders and depressive/anxious mood were seen in more than 50% of the cases. Overall median score of UPDRS “ON” was 36 (22–60), most of our patients were graded in the “2.5” and “3” categories of the H&Y scale with a median for S&E score of 80% (80–90%). The median value for the MoCA test was 22 (17–25). There were no statistically significant differences between sex, regarding these scores.

We were able to establish through multivariate linear regression modeling that an increased disease duration along with the presence of orthostasis, dysphagia, and mood disorders significantly correlated with increased scores in total ON UPDRS. Furthermore, we found an interaction between performing regular physical activity and duration of disease, where despite having increased years of evolution, patients that performed regular physical activity still scored less in the total ON UPDRS (see Supplementary Figure 1). Additionally, lower scores in MoCA testing significantly correlated with increased age, coffee consumption, and the presence of hallucinations, falls, and mood disorders (depression/anxiety), whereas increased years of education correlated with better MoCA scores (see Supplementary Figure 2).

The average coverage of the 10 genes analyzed in this study was >588× for all genes. The coverage per gene and the percentage of nucleotides covered at >15× and >30× for each gene are detailed in Supplementary Table 1. There were no differences in the coverage across the samples (patients and controls). Overall, after quality control, we identified 163 rare variants (Supplementary Table 2) and 158 common variants (Supplementary Table 3) across all genes and all samples that were included in the analysis. Specific protein and DNA changes are listed in Supplementary Tables 4, 5 for rare and common exonic variants, respectively.

Burden and SKAT-O analyses did not identify an association of any of the tested genes and PD (Table 3) after correction for multiple comparisons, as expected given the small sample size. We also did not identify any PD patients with potentially damaging homozygous or compound heterozygous variants in any of these genes. Rare variants in LRRK2 were nominally associated with PD, and 11 (9.2%) patients carried a rare non-synonymous variant, compared to four (4.1%) among the controls. Interestingly, six of these rare non-synonymous variants, all located between amino acids p.1620 and 1623 in the C-terminal-of-ROC (COR) domain of Lrrk2, were found in six patients and none in controls (Table 4).

Non-synonymous GBA variants were identified in three individuals: p.T369M was identified in a male patient with age at onset of 48 years, p.N370S was identified in a healthy female individual recruited at the age of 78 years, and p.L444P was identified in a healthy female individual recruited at the age of 64. While we cannot rule out that these healthy individuals will develop PD in the future, it is unlikely that GBA variants have a major role in PD among Costa Rican patients. One PD patient carried a pathogenic GCH1 variant, p.K224R, further emphasizing the role of this gene in PD.

In the analysis of common variants, none of the variants was associated with PD after correction for multiple comparisons (Supplementary Table 3), which set the corrected p-value for statistical significance at p < 0.00031. One non-synonymous variant in LRRK2, p.I723V, was found with allele frequency of 0.01 in patients and 0.09 in controls (OR = 0.11, 95% CI = 0.02–0.52, p = 0.005), yet this difference was not statistically significant after correction for multiple comparisons.

PD prevalence has been increasing over time with a global age-standardized prevalence rate increase of 21.7% from the years 1990 to 2016 (44). Furthermore, PD prevalence seems to be lower in Eastern compared to Western countries (45). Few studies have explored the prevalence of PD in Latin America providing values that are similar either to other developing countries (46) or to European cohorts (47, 48). PD also becomes more common with advancing age (44, 45). Our sample average age of PD at onset and at diagnosis was lower when compared to other cohorts (49–51), although it could suggest that PD presents earlier in Costa Rica, and more epidemiological studies are needed as it could also be related to recruitment bias.

The majority of our patients fulfilled Group A Gelb criteria while up to 60% also reported at least one of Group B symptoms, the most frequent being dystonia, falls, and dysphagia. The median for years of evolution of the disease for both men and women was 5; thus, we would expect to find Group B criteria in these patients along with the evolution of the disease. Few studies have explored ethnic variations in motor symptoms of PD, suggesting increased atypical features in Black and South Asian PD patients (52, 53); however, there is not enough evidence available along with a lack of standardized methodology to determine motor subtypes across studies and to further establish ethnic patterns of motor features (12). Common non-motor manifestations such as hyposmia, sleep disorders, and depressive/anxious mood were seen in more than half of our PD cases. Regardless of ethnicity, non-motor features are commonly present in PD with subtle differences described. Gastrointestinal non-motor features along with depression seem to be high in East Asian cohorts (54, 55). Likewise, Latino populations, such as Mexican (56), Peruvian (57), and ours, also reported high frequency of mood disorders including depression and anxiety, when compared to studies from UK and USA (58, 59). We also observed in our sample a frequency of sleep disorders and hyposmia that is higher than those reported in other cohorts (12).

Overall, our patients had a low education, which has been previously associated with a higher hazard of incident parkinsonism (60). A reduced education has also been suggested as a risk factor for cognitive impairment in PD (61). A history of non-potable water consumption along with exposure to pesticides and herbicides was reported in up to 40% of our patients. This type of exposure agrees with a mostly rural origin and the fact that 12.2% of the subjects reported involvement in agricultural activities as a main income source. We did not assess the frequency of protective and risk factors in the control group; hence, we are not able to establish any comparison with PD cases. Previous exposure to pesticides and herbicides is associated with the development of PD (62); yet, the identification of a given specific agent and the exact timing and dosing of exposure are almost impossible to establish through observational studies (63, 64). Nonetheless, key work detailing specific mechanisms that render patients vulnerable to pesticide-induced injury has been elegantly shown in animal models, further establishing biologic and toxicological pathways for specific chemicals to potentially cause PD (65). A similar situation is present regarding the exposure to welding and heavy metals. Manganese (66), copper, iron (67), and mercury (68) have been proposed as possible agents associated with the development of PD. In this study, 22.1 and 11.6% of the patients reported frequent exposure to welding and other heavy metals, respectively; however, the exact timing and dosing of exposure was not possible to assess.

Other literature has underscored the presence of protective factors for PD development, among which the most notable and with the strongest evidence include tobacco (69) and coffee consumption (7, 70–73). For both protective factors, there is also a dosing effect described, where the protective effect increases along with an increasing exposure (74, 75). Paradoxically, over 90% of our PD cases had been exposed to a protective factor in the past, most of them having a regular coffee intake (two to three cups per day for over 15 years), and yet they all developed PD.

Performing regular physical activity correlated with lower ON UPDRS scores in spite of increasing age. Physical activity has been established as a possible protective factor for incident Parkinsonism (76); our data would suggest that physical activity could determine reduced severity of disease, specifically concerning motor features. Although exercise has not been proven to slow the progression of akinesia, rigidity, and gait disturbances, it promotes a feeling of physical and mental well-being, and at the same time, it can alleviate rigidity-related pain and improve patients' motor (77) and non-motor symptoms (78).

Increasing age, coffee consumption, hallucinations, falls, and mood disorders along with reduced years of education significantly correlated with worse MoCA scores. Older age and duration of PD are determinant risk factors for incidence of dementia in PD (79). Furthermore, hallucinations have been established as risk factors for cognitive impairment (79, 80) along with gait disturbances (manifested by falls) (81) and depression (82). Reduced education years also have been proposed as a risk factor for cognitive impairment in patients with PD (61). Poor global cognition has been previously associated with a higher risk of incident parkinsonism (60). Coffee consumption has been suggested to reduce risk of dementia (83) with a dosing effect (84, 85); however, there have been inconsistent findings regarding the effects of coffee consumption on specific cognitive domains. It has been suggested to be in association with improved executive performance but smaller hippocampal volume and worse memory function (86); nonetheless, this association is not sustained when cognition is analyzed longitudinally. Other literature suggested that coffee might be slightly beneficial on memory without a dose–response relationship (87). Recent large-scale genetic analysis using mendelian randomization did not find any evidence supporting any beneficial or adverse long-term effect of coffee consumption on global cognition or memory function (88) or AD incidence (89). To our knowledge, there is no literature evaluating the effect on cognition of coffee consumption, specifically for PD patients. Our findings suggest a possible deleterious effect that should be further explored in this population.

After sequence coding familial PD and atypical parkinsonism-associated genes including GBA, SNCA, VPS35, LRRK2, GCH1, PRKN, PINK1, DJ-1, VPS13C, and ATP13A2 and correcting for multiple comparisons, burden and SKAT-O analyses did not show an association of any of the tested genes and PD. We also did not identify any homozygous or compound heterozygous pathogenic variants in any of these genes.

Non-synonymous GBA variants were identified in three individuals including one patient and two unaffected controls. While we cannot rule out that these healthy individuals will develop PD in the future, it is unlikely that GBA variants have a major role in PD among Costa Rican patients especially when compared to other European and Ashkenazi Jewish populations where we find that 8–20% of the patients harbor GBA variants (28).

Finally, one PD patient carried a pathogenic variant, p.K224R, in the GCH1 gene. GCH1 encodes for GTP cyclohydrolase 1, which is a key enzyme for dopamine production in nigrostriatal neurons. Loss-of-function mutations such as p.K224R have been shown to cause Dopa-responsive dystonia (DRD); however, variants in this gene have also been implicated in PD, perhaps through regulation of GCH1 expression (90, 91). It has been suggested that late-onset DRD might present clinically with parkinsonism, or alternatively, pathogenic GCH1 mutations may predispose to both diseases and carriers will develop any or both depending on other genetic or environmental factors (92). Our patient did not present clinical features suggestive of DRD and did not have any family history of PD.

Rare variants in LRRK2 were nominally associated with PD, observed only in affected individuals; six of these rare non-synonymous variants were located between amino acids p.1620 and 1623 in the COR domain of Lrrk2. LRRK2 encodes a multiple domain protein that includes a Roc-COR tandem domain, a tyrosine kinase-like protein kinase domain, and at least four repeat domains located within the N-terminal and C-terminal regions. The Roc-COR domain classifies the Lrrk2 protein as part of the ROCO superfamily of Ras-like G proteins (93). Mutations in LRRK2 are the most common cause of late-onset hereditary PD. Most frequently reported disease-causing mutations are located in the kinase domain (i.e., G2019S), increasing kinase activity, and in the Roc-COR tandem domain (i.e., R1441C/G and Y1699C), impairing its GTPase function. Alterations of both kinase and GTPase activity may mediate neurodegeneration in these forms of PD (94). Of the six patients found to have non-synonymous variants in the COR domain, two had first-degree relatives with dementia, one had a second-degree relative with PD, and one had two sisters with PD diagnosed at a very young age (20 and 30 years old) (see Table 4).

Methodological issues, such as size and composition of the sample (i.e., number of familial and sporadic cases), might explain the variations seen in the frequency of LRRK2 mutations in case series from similar countries. However, there is a clear difference established among geographical regions, where North African Arabs (95), Ashkenazy-Jews (96) and certain Europeans cohorts (97–99) might report a higher prevalence than Latin American and Asian populations for these mutations (15, 100, 101).

The non-synonymous missense mutations described here are all found in the COR domain of Lrrk2. To gain insight into how these mutations may affect the function of Lrrk2, we investigated their locations in the structure of Lrrk2. The high-resolution cryoelectron microscopy (cryoEM) structure of the C-terminal domains of Lrrk2 in different states have recently been reported and shed light on how allosteric interactions between different domains regulate microtubule interactions (102). The structure notably shows interactions between the ROC GTPase domain and the COR-B domain, notably involving the pathogenic mutation sites p.Arg1441 and p.Tyr1699 (Figure 1A). These interdomain interactions enable the kinase activity to be regulated by GTP binding to the ROC domain. The mutations described here, found in the segment a.a. 1619–1623, are all located in a loop of the COR-A domain. This loop, which spans a.a. 1613–1624, is disordered in the cryoEM structure, and thus, no atomic resolution model is available for that segment (Figure 1A). It is therefore not possible to gain detailed insights into the effect of each individual missense mutation.

However, integrative modeling, based on cryoelectron tomography (cryoET) data collected from in situ and in vitro-reconstituted Lrrk2 filaments bound to microtubules, shows how the different domains of Lrrk2 dimerize and associate with microtubules (102, 103). Dimerization is mediated via two sites through reciprocal interactions: one involving WD40–WD40 interactions and another one involving COR–COR interactions. These interactions enable Lrrk2 C-terminal domains to form extended oligomeric filaments that form a helix around the microtubule. Of particular interest here, the COR–COR dimerization interface involves both the COR-A and COR-B domains, with the loop containing a.a. 1613–1624 at the center of this interface (Figure 1B). Mutations in this loop may thus affect dimerization. Given that the kinase activity and conformation affect the ability of Lrrk2 to dimerize through the COR domain via allosteric interactions, it is possible that mutations in the COR-A loop in turn affect the kinase activity. Further experiments would be required to determine how the mutations described here affect the kinase, dimerization, and microtubule-binding activity of Lrrk2.