Pleiotropic effects for Parkin and LRRK2 in leprosy type-1 reactions and Parkinson's disease.

Type-1 reactions (T1R) are pathological inflammatory episodes and main contributors to nerve damage in leprosy. Here, we evaluate the genewise enrichment of rare protein-altering variants in 7 genes where common variants were previously associated with T1R. We selected 474 Vietnamese leprosy patients of which 237 were T1R-affected and 237 were T1R-free matched controls. Genewise enrichment of nonsynonymous variants was tested with both kernel-based (sequence kernel association test [SKAT]) and burden methods. Of the 7 genes tested 2 showed statistical evidence of association with T1R. For the LRRK2 gene an enrichment of nonsynonymous variants was observed in T1R-free controls (P SKAT-O = 1.6 × 10-4). This genewise association was driven almost entirely by the gain-of-function variant R1628P (P = 0.004; odds ratio = 0.29). The second genewise association was found for the Parkin coding gene PRKN (formerly PARK2) where 7 rare variants were enriched in T1R-affected cases (P SKAT-O = 7.4 × 10-5). Mutations in both PRKN and LRRK2 are known causes of Parkinson's disease (PD). Hence, we evaluated to what extent such rare amino acid changes observed in T1R are shared with PD. We observed that amino acids in Parkin targeted by nonsynonymous T1R-risk mutations were also enriched for mutations implicated in PD (P = 1.5 × 10-4). Hence, neuroinflammation in PD and peripheral nerve damage due to inflammation in T1R share overlapping genetic control of pathogenicity.

Leprosy is characterized by chronic infection with Mycobacterium leprae where the progression to disease is strongly dependent upon the host genetic background (1–6). One focus of current leprosy control efforts is the prevention of nerve damage (7). A major contributor to nerve damage are sudden-onset episodes of excessive inflammatory responses termed type-1 reactions (T1R) (8). Without timely intervention T1R can cause irreversible nerve damage due to a pathological cell-mediated inflammatory response directed against host peripheral nerve cells (9). Depending on the epidemiological setting, T1R episodes can afflict 30 to 50% of leprosy cases but are more common among older patients (10). While approximately one-third of T1R are diagnosed at the time of leprosy diagnosis, T1R episodes can occur years after successful completion of therapy presumably in the absence of M. leprae bacilli (10, 11).

Host genetic predisposition is an important aspect of T1R pathogenesis. A prospective study of host responses to M. leprae antigen identified a gene-set signature that predisposed to T1R (12). In addition, candidate gene approaches identified 6 protein-coding genes in association with T1R (13–19). More recently a genome-wide association study identified a long noncoding RNA on chromosome 10 as T1R risk factor in independent populations and also provided suggestive evidence for PRKN (formerly PARK2) as a T1R susceptibility gene (18). Here, we tested if nonsynonymous variants with a likely functional impact at the protein level were enriched in the 7 coding genes previously associated with T1R. We identified protein-altering variants in the PRKN and LRRK2 genes as T1R risk modulating factors. Since PRKN and LRRK2 are established Parkinson’s disease (PD) susceptibility genes, we tested if T1R risk variants had been implicated in PD. Indeed, we found T1R risk modulating amino acid changes significantly enriched among PD susceptibility variants. Loss-of-function mutations within a common subset of amino acids in Parkin were shared risk factors for T1R and PD while a gain-of-function amino acid substitution in LRRK2 had an antagonist pleiotropic effect and was protective for T1R and a risk factor for PD.

Results

We selected from our sample of Vietnamese leprosy patients 237 T1R-affected subjects (cases), which we matched by age at leprosy diagnosis, gender, and clinical subtype of leprosy with 237 T1R-free subjects as controls for deep resequencing of 7 T1R susceptibility genes (). We focused our analysis on the 63 variants that represented nonsynonymous changes, or insertion/deletions, and tested those for association with T1R. None of the studied subjects were homozygous or compound heterozygotes for nonsynonymous rare variants in any given gene. Of the 7 genes, only nonsynonymous changes in PRKN and LRRK2 showed significant evidence for association with T1R (Table 1).

Table 1.

Genewise association of protein-altering variants in T1R-risk loci

Gene	Var	Sing	No MAF cutoff			MAF < 0.05
			PSKAT-O	pVT	Var	PSKAT-O	pVT
TLR1	15	6	0.69	0.62	13	0.68	0.39
TLR2	9	7	0.11	0.47	9	0.11	0.45
PRKN	9	6	7.4 × 10−5	0.01	7	1.1 × 10−3	5.8 × 10−3
TNFSF15	2	0	0.17	0.41	2	0.18	0.44
LRRK2	20	10	1.6 × 10−4	2.2 × 10−4	17	1.5 × 10−4	1.4 × 10−4
NOD2	10	5	0.61	0.97	10	0.62	0.98
Combined	65	34	0.04	0.34	58	0.04	0.37

Var, number of nonsynonymous variants found for the gene; Sing, number of singleton variants found in only 1 individual; MAF, minor allele frequency.

The LRRK2 R1628P Variant Is Associated with Protection from T1R.

Nonsynonymous LRRK2 variants were significantly enriched among T1R-free controls (Fig. 1 and Table 1). Common variants did not contribute to this effect since restriction of the analysis to variants with minor allele frequency (MAF) < 0.05 did not cause a decay in the strength of LRRK2 association with T1R (Table 1). In contrast, the 2 low-frequency variants P755L and R1628P were highly enriched (>3-fold excess) in T1R-free subjects. When we performed single-variant analyses we observed a strong protective effect of the LRRK2 R1628P mutation for T1R (P = 0.004; odds ratio [95% CI] = 0.29 [0.13 to 0.67] for 1628P carriers). None of the other variants, including P755L, displayed significant univariate association with T1R (). We conditioned the genewise association of LRRK2 with T1R on the R1628P mutation and observed a large drop in the strength of association (PSKAT-O = 1.6 × 10−4 to PSKAT-O = 0.01). By also adjusting on the P755L variant, LRRK2 association with T1R protection became nominally nonsignificant (PSKAT-O = 0.06). However, given the remaining trend in favor of association it is possible that additional rare LRRK2 variants contributed to protection from T1R. Of note, the MAF for the R1628P was similar between leprosy cases (T1R-affected + T1R-free) and 472 Vietnamese healthy subjects that were genotyped independently (3.4% vs. 3.3%, respectively; P = 0.83), supporting a T1R specific effect.

Fig. 1.

Overview of the deep resequencing for the LRRK2 and PRKN genes. At the top, the population mean depth of coverage is plotted according to the exons encoding the LRRK2 gene (A) and PRKN (B). Blue circles indicate the depth of coverage per base pair with the light blue shade representing the SD of the mean. A red horizontal line marks the average depth of coverage for the 2 genes in the studied population. The protein functional domains are shown in the center and linked to their respective coding exons. The location of known coding mutation of LRRK2 and Parkin reported by either dbSNP or gnomAD, mutation associated with Parkinson’s disease in PDmutDB and LOVD and mutations associated with T1R are shown as yellow with blue bars denoting variants that overlap between PD and T1R. At the bottom, the allele count or the frequency of LRRK2 and PRKN mutations is shown for the T1R-affected and T1R-free samples.

To validate the detected LRRK2–T1R associations we decided to test for a possible impact of the 1628P and 755L mutations on LRRK2 functional activity. Hence, we genome-edited RAW cells to express hetero- and homozygously the 1628P and 755L mutant forms of LRRK2 (). We also generated a LRRK2 knockout line as negative control. Mutant and wild-type RAW cells were infected with Mycobacterium bovis (bacillus Calmette–Guérin) or M. leprae. Infection with bacillus Calmette–Guérin induced LRRK2 expression in all cells, but no significant difference of expression levels was observed between the LRRK2 wild-type and mutant proteins (Fig. 2 ). Following infection with bacillus Calmette–Guérin 2 independent 1628P/P mutant cell lines (C1 and C2) mediated a significantly stronger respiratory burst while production of reactive oxygen species (ROS) by the 755L line was indistinguishable from the response of wild-type cells (Fig. 2 ). Compared with wild-type cells, the 1628P mutant also triggered a stronger respiratory burst in response to M. leprae (Fig. 2). Interestingly, the response of heterozygous LRRK2 1628P/R cells was indistinguishable from the one mounted by homozygous 1628P/P cells (Fig. 2). Next, we evaluated if the LRRK2 mutants impacted the proinflammatory response to bacillus Calmette–Guérin infection by measuring levels of tumor necrosis factor (TNF). We observed that TNF secretion was strongly up-regulated in 1628P/P cells after 6 h of bacillus Calmette–Guérin infection (). However, at the 24-h timepoint TNF levels were no longer statistically different between wild-type and 1628 P/P mutant cells ().

Fig. 2.

The LRRK2 R1628P rare variant displays an increased respiratory burst in response to challenge with both live bacillus Calmette–Guérin and M. leprae. A and B show the relative abundance of LRRK2 for wild-type (WT) RAW cells and 3 CRISPR/Cas-edited mutant cell lines in both resting state and under bacillus Calmette–Guérin stimulated conditions in a representative Western blot (A) and as triplicates with SDs (B). C–F show the kinetics of total reactive oxygen species (ROS) produced by RAW cells upon challenge with bacillus Calmette–Guérin (C and D) or M. leprae (E and F). LRRK2 knockout (KO) cells and the nonsynonymous mutations (1628P/P, 1628R/P, and 755L/L) were compared with WT using two-way ANOVA. *P < 0.05; **P < 0.01 and ***P < 0.001. The experiments were repeated at least 3 times with similar results. Two independent clones for the 1628P/P were tested with similar results. Compared with WT, LRRK2 KO cells consistently produced less ROS upon bacterial challenge with bacillus Calmette–Guérin or M. leprae (C–F). In contrast, there was no significant difference in the kinetics of ROS production following exposure to bacillus Calmette–Guérin between WT and LRRK2 755L/L-expressing cells (D). Compared with cells expressing LRRK2 WT, cells expressing the LRRK2 1628P/P variant consistently showed increased ROS production following exposure to both bacillus Calmette–Guérin and M. leprae (C and E). Compared with LRRK2 WT, cells heterozygotic for LRRK2 1628R/P also displayed higher ROS production upon M. leprae stimulation which was not significantly different from the response of homozygous LRRK2 1628P/P carriers.

Since LRRK2 had been linked to protection from apoptosis, we assessed the impact of hetero- and homozygotic clones of 1628P as well as 755L LRRK2 homozygous cells on bacillus Calmette–Guérin-triggered apoptosis. We infected cells expressing LRRK2 wild type, 1628P/P (C1 and C2), 1628R/P, 755L/L LRRK2 mutant forms, and LRRK2 knockout cells with bacillus Calmette–Guérin for 24 h. We then determined in bacillus Calmette–Guérin-infected and uninfected controls the extent of apoptosis. Following infection with bacillus Calmette–Guérin for 24 h there was an induction of apoptosis in cells of all lines (Fig. 3). Compared with LRRK2 knockout cells, all other cell lines displayed a reduction in the extent of apoptosis triggered by bacillus Calmette–Guérin. Cells expressing wild-type and 755L forms of LRRK2 were indistinguishable by their bacillus Calmette–Guérin-triggered apoptotic activity. However, both clones expressing the 1628P mutant form of LRRK2 showed significant abrogation of apoptosis relative to wild-type cells (Fig. 3). Interestingly, we observed that the 1628P allele effect on apoptosis inhibition was significantly additive (trend test: P = 0.003 including LRRK2 KO, and P = 0.04 without LRRK2 KO). Hence, our results revealed a consistent impact of LRRK2 on bacillus Calmette–Guérin-induced apoptosis and identified LRRK2 1628P as gain-of-function variant and strong apoptosis-reducing factor.

Fig. 3.

The LRRK2 R1628P mutation abrogates apoptosis induced by live bacillus Calmette–Guérin. Fluorescence-activated cell sorting (FACS) was used to evaluate the impact of LRRK2 nonsynonymous variants on apoptotic cell death upon bacillus Calmette–Guérin challenge. (A) Cellular staining as revealed by FACS analyses for RAW cells expressing LRRK2 WT and genome-edited proteins after 24 h of coincubation with live bacillus Calmette–Guérin and unchallenged baselines. (B and C) Proportion of cells undergoing apoptosis. The bar plots display the mean estimate for 2 experiments done in triplicate with their respective SEs. There were no significant differences in the proportions of live or apoptotic cells under unstimulated condition. Infection with bacillus Calmette–Guérin triggered a significant (P < 0.001) increase of apoptosis in wild-type, 755L LRRK2, and LRRK2-depleted (KO) cells with a corresponding decline of live cells. In contrast, 2 independent clones expressing 1628P/P LRRK2 protein did not respond with a significant increase of cell death or apoptosis. When stimulated with bacillus Calmette–Guérin, cells depleted for LRRK2 (KO) expressed significantly higher apoptosis compared with mutant clones or wild-type cells (P < 0.001). The levels of total apoptosis showed a significantly dose-dependent trend with the 1628R/P line situated between wild-type and homozygous 1628P/P lines (P = 0.003 including LRRK2 KO and P = 0.04 without LRRK2 KO). Statistical tests were performed using one-way ANOVA with Tukey correction and a linear mixed model for the trend test.

We had previously identified the common M2397T amino acid change as a T1R risk factor (19). Interestingly, we observed in the present sample that the 1628P protective allele was always on the background of the T1R-protective 2397T allele (D′ = 1). To evaluate the independent effect of the 2397M variant, we removed the effect of the 1628–2397 P–T haplotype (). In this case, the risk effect of the 2397M allele was highly nonsignificant (P = 0.40), indicating that the LD with the 1628P allele might explain, at least to some extent, the previously detected effect of the 2397M allele.

Rare Nonsynonymous Variants in PRKN Are a T1R Risk Factor.

For the PRKN gene we observed 9 missense variants of which 2, rs1801474 (S167N) and rs1801582 (V380L), were common and 7 were rare, including the uncatalogued Parkin mutation K299R (Fig. 1). PRKN rare variants were observed only in T1R-affected cases, resulting in strong evidence for genewise association of PRKN with T1R susceptibility (Fig. 1). To evaluate if the PRKN association with T1R was only driven by rare variants, we reanalyzed the data removing the 2 common PRKN polymorphisms and observed an approximately 1-log drop in strength of association (PSKAT-O = 7.4 × 10−5 to PSKAT-O = 1.1 × 10−3). Next, we tested if the 2 common PRKN missense variants were associated with T1R in a single-variant model. We observed a trend for association of T1R and the Parkin V380L variant (P = 0.09; odds ratio [95% CI] = 1.44 [0.94 to 2.20] for 380L carriers), but no association of T1R and S167N (P = 0.22; ). To confirm that the V380L variant contributed to the genewise association we conditioned the genewise model with all nonsynonymous variants of PRKN on the V380L mutation. We observed the same drop in the P value as in the genewise analysis with only rare variants (PSKAT-O = 7.4 × 10−5 to PSKAT-O = 5.7 × 10−3). This observation indicated that although the rare variants were the main driver of the PRKN association with T1R the common V380L variant also contributed to T1R risk.

The PRKN rare variants were dispersed across the primary structure of Parkin, with 3 mutations in the central RING1 domain, 1 in the Ubl domain, 1 in the IBR domain, and 2 in the REP linker between the IBR and RING2 (Fig. 1). Alterations in 3 Parkin residues (D243N, R275Q, and R366W) were predicted to be damaging in all 5 databases curated (Dataset S1). Two other substitutions (R42C and K299R) were predicted to impair Parkin function in more than 1 database while the remaining 2 rare variants located toward the terminal end of Parkin (R396G and K408R) were predicted to be tolerable (Dataset S1). Structural analyses of the 5 amino acid changes shared by PD and T1R cases predict impairment in either phosphorylation by PINK1 (R42C), ubiquitination activity (D243N), or stability (R275Q, R366W, and R396G) (). To assess the functional impact of Parkin mutations, we tested these 5 amino acid changes for protein stability, ubiquitination, phosphorylation, and mitophagy (Fig. 4 and ). Using a thermal shift assay, all mutations except D243N showed a significant lower protein melting point compared with wild-type Parkin (Fig. 4). The Parkin R275Q and R366W mutations had the lowest melting temperature while mutations R42C and R396G had a less-pronounced impact on protein stability (Fig. 4). The variant R42C decreased Parkin phosphorylation by PINK1 while R275Q and R366W increased phosphorylation (Fig. 4 ). The R42C and D243N variants displayed significantly reduced ubiquitination activity against UbcH7 in the presence of PINK1 (Fig. 4 ). In addition, the polyubiquitin smear of the R275Q variant was clearly different from wild type, with shorter chains being formed (Fig. 4). The D243N and R275Q mutations showed significantly reduced mitophagy in U2OS cells transfected with green fluorescent protein (GFP)-Parkin mutants (Fig. 4 ). The R275Q mutation, and to a lesser extent R396G, showed lower protein levels in U2OS cells (Fig. 4). Of note, different Parkin mutants displayed variable impact in the different assays used to probe Parkin function.

Fig. 4.

Characterization of Parkin mutations shared by PD and T1R cases. (A) Melting temperatures (T) of recombinant human Parkin mutants, measured using differential scanning fluorescence. (B) Mass spectrometry analysis of Parkin phosphorylation by PINK1. Human Parkin (20 μM) is phosphorylated by PINK1 (1 μM) for 15 min in the presence of ATP (5 mM) and ubiquitin (20 μM). The arrows indicate a 79.9-amu (atomic mass unit) shift resulting from phosphorylation. (C) The fraction (percent) of Parkin and ubiquitin phosphorylation was calculated from the MS peak intensities. (D) SDS/PAGE analysis of Parkin ubiquitination reactions following PINK1 phosphorylation. In step 1, human Parkin (2 μM) is phosphorylated by PINK1 (0.1 μM) for 15 min in the presence of ATP (5 mM) and ubiquitin (50 μM). In step 2, E1 (50 nM) and UbcH7 (2 μM) are added to allow ubiquitination for another 15 min at 37 °C. Parkin auto-ubiquitinates to form smears and also ubiquitinates UbcH7. (E) Ubiquitination activity was quantified using loss of the unmodified UbcH7 band, compared with the no ATP control. (F) Mitophagy was examined using a FACS-based analysis of mitochondrially targeted Keima (mt-Keima). Representative FACS data of mt-Keima expressing WT or PRKN mutants untreated (control) or treated with CCCP (20 μM) for 4 h from at least 3 independent experiments. (G) Bars showing mitophagy percentage for either untreated (control) or after 4 h of CCCP treatment quantified from mt-Keima signal in U2OS cells expressing GFP-Parkin variants normalized to wild-type (WT) Parkin. (H) Bars showing GFP intensity of untreated cells expressing GFP-Parkin variants normalized to WT Parkin. The quantification of average percentage for the mitophagy was a combination of at least 3 independent experiments. P values for all of the plots were determined by one-way ANOVA with Dunnett’s post hoc tests. *P < 0.05; **P < 0.01; ***P < 0.001. n.s., not significant.

Pleiotropic Effect of LRRK2 and PRKN Amino Acid Changes on T1R and PD.

Since mutations in LRRK2 and PRKN are main causes of PD we assessed if the amino acid substitutions identified in our study were curated in 2 PD mutation databases (20, 21). The main evidence for association between LRRK2 and T1R was provided by variants P755L and R1628P and, while not causally linked with PD, both variants are risk factors for PD (Table 2). Interestingly, antagonistic pleiotropy of the 1628P variant was also observed for PD and Alzheimer’s disease (22). For the PRKN gene, 5 of the 7 amino acid mutations identified as risk factors for T1R had previously been implicated in PD (Fig. 2 and Table 2). Moreover, 1 of the 7 variants had not been detected in PD or healthy subjects. This shared co-occurrence of amino acid mutations between PD and T1R was strongly nonrandom (P = 8.7 × 10−4; ). When we disregarded the uncatalogued variant sharing became more significant (P = 1.5 × 10−4).

Table 2.

Reported association for variants observed in our study

Gene	Mutations	Phenotype reported	gnomAD	1000 Genomes
			All pop	EAS	KHW
	R42C (rs577232474)	PD (54, 55); EOPD (R42P) (56–60); glioblastoma (42)	0.0002	0	0
	S167N (rs1801474)	PD (61, 62)	0.07	0.39	0.29
	D243N (rs146173584)	PD (63, 64); lung cancer (42)	0.0003	0	0
	R275Q (rs769230602)	EOPD (R275W) (37, 38, 65); glioblastoma (42)	0.00002	0	0
PRKN	K299R (novel)	PD (I298S) (38)	0	0	0
	R366W (rs56092260)	PD (66)	0.0003	0.003	0.005
	V380L (rs1801582)	EOPD (67–69)	0.16	0.08	0.08
	R396G (rs539917500)	EOPD (70)	0.00007	0.001	0.005
	K408R (rs562362828)	—	0.000004	0.00006	0.005
LRRK2	N551K (rs7308720)	PD (25), Crohn’s disease (25)	0.08	0.10	0.07
	K616R (rs281865042)	PD (71)	0.000004	0.00006	0
	P755L (rs34410987)	PD (72, 73)	0.0007	0.009	0.01
	P1446L (rs74681492)	PD (74, 75)	0.0002	0.003	0.005
	R1628P (rs33949390)	PD (28, 76–78); essential tremor (79); Alzheimer’s disease (22, 80)	0.002	0.02	0.05
	S1647T (rs11564148)	PD (28)	0.29	0.34	0.36
	N2081D (rs33995883)	PD (25, 81); Crohn’s disease (25)	0.016	0.0005	0
	G2385R (rs34778348)	PD (28, 81)	0.002	0.02	0.01

All pop, allelic frequencies estimated in 141,456 individuals from all cataloged populations in gnomAD; EAS, allelic frequencies for East Asian populations in the 1000 Genomes database; KHW, allelic frequencies for the Kinh in Ho Chi Minh City, Vietnamese population in the 1000 Genomes database.

Discussion

Despite the availability of an efficacious treatment for leprosy, the incidence of leprosy has shown only a moderate global decline over the past 15 years (23). This unexpected epidemiological situation has given rise to 2 major areas of research. The first theme concerns the identification of the unknown mode of transmission of the M. leprae bacillus. The second theme is focused on the early detection of signs of leprosy and the prevention of nerve damage in leprosy patients (7). A main contributor to nerve damage are T1R and the question remains why only a minority of leprosy patients develop T1R and how to identify those who are at risk. While the clinical manifestation of this predisposition requires the trigger by M. leprae it is noteworthy that T1R can occur after microbiological cure of leprosy.

Acute episodes of T1R are characterized by a cytokine storm which is indicative of dysregulated inflammatory responses. Prospective transcriptomic analyses of the immune response to M. leprae antigen and a genome-wide association study support the central role of host genetics in the breakdown of coordinated gene regulation in the T1R phenotype (12, 18). Additional candidate gene T1R association studies identified T1R susceptibility genes that were investigated in the present study. However, these previous genetic studies focused on the role of common genetic variants which are more likely to impact gene expression levels and not the functional capacity of encoded proteins. Here, we show that for 2 of 7 T1R susceptibility genes, uncommon amino acid mutations that impact protein function are significant contributors to T1R risk that act independently of common regulatory variants. Our findings that a gain-of-function mutation in the LRRK2 protein is protective for T1R while loss-of-function variants in Parkin are T1R risk factors are consistent with the concept that a genetically controlled miscommunication between pro- and antiinflammatory host responses underlies the T1R phenotype.

Of the 7 T1R candidate genes studied (), we identified rare and low-frequency nonsynonymous variants only in the E3 ligase Parkin and the polyfunctional kinase LRRK2 as significantly associated with M. leprae-triggered pathological inflammation. While the PRKN mutations were pleiotropic for T1R and PD, the 1628P LRRK2 mutation was protective for T1R but had previously been reported as risk-variant for PD, and hence displayed antagonistic pleiotropy for the 2 disease entities. The 1628P variant has been shown to be a gain-of-kinase-function mutation (24). Interestingly, the LRRK2 G2019S and N2081D gain-of-kinase-function mutations that are risk factors for PD displayed defects in autophagic cargo clearance and lysosome acidification relative to LRRK2 wild type (25). It is therefore tempting to speculate that the 1628P mutation has similar impact on LRRK2 function. This may imply that T1R is more strongly impacted by gain-of-function pathways (e.g., apoptosis and ROS production) while PD is more strongly impacted by loss of function pathways (e.g., autophagy and lysosome acidification).

One of the unexpected findings was the observation that the 1628P mutation blocked bacillus Calmette–Guérin-triggered apoptosis. However, in sepsis, another disease of dysregulated inflammation, increased apoptosis of key immune cells underlies the pathological inflammatory response and it is possible that similar effects may occur in T1R (26, 27). For example, the increased the respiratory burst of RAW cells in response to bacillus Calmette–Guérin is consistent with a previous report that the LRRK2 1628P allele promoted increased kinase activity (28). ROS generated during a respiratory burst down-regulate inflammation (29), and systemic ROS does induce apoptosis which is generally considered an antiinflammatory cellular response since apoptotic cells release multiple antiinflammatory mediators (30). However, if apoptotic cells are not cleared efficiently, this may lead to increased inflammation (31). This dual role of apoptosis may underlie the antagonistic pleiotropy of 1628P in PD and T1R. It is possible that the reduction of antiinflammatory molecules resulting from abrogated apoptosis mediated by the 1628 mutant protein is disease-promoting for PD patients while the resulting lower yield of apoptotic debris in leprosy patients may protect against T1R. A previous prospective study of leprosy patients at risk for T1R identified higher pro- and antiinflammatory gene expression in those patients destined to suffer T1R (12). This suggested a breakdown in the regulation of the inflammatory response in T1R patients, opening the possibility that the 1628P mutation, at least in the periphery, reestablishes the equilibrium between pro- and antiinflammatory pathways.

An interesting finding was the complete linkage disequilibrium between the LRRK2 1628P and 2397T alleles. In an Ashkenazi Jewish sample, the PD and Crohn’s disease protection allele 551K also occurred on the background of the 2397T allele as the 1628 P mutation in Vietnamese (25). However, the N551K had no effect on T1R in our study, suggesting 2 independent effects of amino acid changes in LRRK2 (N551K and R1628P) that are captured by another common functional variant, i.e., M2397T.

Parkin is an evolutionary conserved modulator of innate immunity to infectious agents (32, 33) and genetic variants of PRKN were associated with increased susceptibility to both leprosy and typhoid fever in ethnically diverse populations (2, 34, 35). Moreover, mutations in Parkin are the most common known cause of early-onset PD (EOPD) (36). However, it was not resolved if the role of Parkin in infectious and neurodegenerative diseases reflected the promiscuous nature of E3 ligase functions or was based on shared mechanisms of dysfunctional host immunity. Here, we showed that rare amino acid changes associated with T1R in Parkin are predicted loss-of-function mutations that display a striking overlap with established risk factors for PD. Indeed, 5 out of the 7 Parkin rare variants associated with T1R susceptibility were substitutions at the same amino acid residue as observed in PD, while a sixth uncatalogued T1R risk variant (K299R) was adjacent to the annotated I298L and I298S PD mutations (Table 2) (37, 38). As the shared Parkin mutations are mainly connected to EOPD our findings strongly suggested that EOPD shares a mechanistic overlap with T1R, an inflammatory disorder triggered by M. leprae.

To assess the functional relevance of the PRKN mutations we employed a series of commonly used assays of Parkin function. Given the spectrum of Parkin substrates, it is possible that the PRKN mutations also impact aspects of Parkin function that were not captured by the employed assays. However, all Parkin variants evaluated showed significant deviation from wild type in at least one of the assays. The R42C mutant displayed only a minor reduction in protein stability but reduced Parkin autoubiquitination. The reduced activity is most likely the result of impaired PINK1-dependent phosphorylation, which is required for Parkin activation (39). Furthermore, we had previously found that the charge-reversal mutation R42E in the Ubl impairs phosphorylation, since Arg42 is required for binding PINK1 (40). The R275Q and R366W mutations displayed a more pronounced impact on thermal stability and were linked to increased phosphorylation, probably because these mutations induce the dissociation of the Ubl domain that makes it more available for phosphorylation. The D243N variant reduced E2 ubiquitination and displayed impaired mitophagy. Asp243 is located in close proximity to 2 arginine residues in the E2 enzyme UbcH7 (39), and we had previously observed a strong reduction in activity for the D243A mutant (41). The R275Q mutation had the strongest impact on mitophagy. The side chain of Arg275 interacts with a helix in RING1 that acts as an allosteric switch between the inactive (Ubl-bound) and active (pUb-bound) forms (). The R275Q mutation therefore strongly impacts protein stability, which affects its levels in cells and perhaps also its binding to phospho-ubiquitin and ability to activate. Finally, the R396G mutation, located in the Repressor Element of Parkin (REP) helix, has a mild impact on thermal stability and cellular level, as predicted from the structure. Among the mutations identified by our study, amino acid changes at Parkin positions 42 and 275 are the most common PRKN point mutations found in PD cases. Moreover, R42C and R275Q somatic mutations have been associated with higher cell proliferation in cancer (42). Taken together, these results supported the overall conclusion of our study that T1R and PD share a set of molecular lesions in the PRKN and LRRK2 genes.

While mutations within the same amino acids in Parkin are implicated in neurodegeneration in the brain and nerve cell destruction in the periphery, the question remains if related infectious triggers are necessary for the manifestation of both disease phenotypes. Parkin, together with PINK1, is critically important for mitochondrial quality control by modulating mitophagy, a form of autophagy which delivers damaged mitochondria to lysosomes for degradation. However, there is a large body of experimental evidence which extends the role of PRKN and PINK1 in mitochondrial quality control beyond the control of mitophagy (43). Given that mitochondria are key players of adaptive and innate immunity, and the suggested role of Prkn and Pink in autoimmune-mediated neuronal destruction in a mouse model, it is possible that unrecognized infection provides the missing link between impaired mitochondrial quality control and dysregulated neuroinflammation in EOPD (44). The genetic link between EOPD and T1R leads to one additional consideration. Since a majority of T1R cases present within 3 years of diagnosis and treatment of leprosy it seems prudent to consider the role of temporally remote and largely resolved infections as possible triggers for some cases of EOPD. Careful screening for the immunological footprint of such remote infections in EOPD cases may provide additional support for a role of infections in EOPD.

Methods

Population Sample and Deep Resequencing.

Written informed consent was obtained from all subjects participating in the study. Minors included in the study assented to participation while a parent or guardian gave informed consent on their behalf. The study was approved by the Research Ethics Board at the Research Institute at McGill University Health Centre in Montreal (REC98-041) and the regulatory authorities of Ho Chi Minh City in Vietnam (So3813/UB-VX and 4933/UBND-VX). We selected from our records of Vietnamese leprosy patients 240 T1R-affected subjects as our cases. Since borderline leprosy subtype and older age at leprosy diagnosis are strong T1R risk factors (10), we selected 240 T1R-free subjects as controls to match the proportion of leprosy subtypes and the age distribution of T1R-affected cases (). For both T1R cases and leprosy controls, 3 samples failed quality control. The paired design was used to avoid confounding the genetic control of T1R with genetic control of leprosy subtype. The follow-up for leprosy-affected individuals without T1R was greater than 5 y. An Illumina TrueSeq Custom Amplicon v1.5 was used to create libraries containing the exonic sequences of 7 targeted genes. Individual indexed libraries were pooled in batches of 96 samples randomly mixing cases and controls for paired-end sequencing using Illumina MiSEq. 600 cycles kit v3. Next, the paired-end sequences were aligned to the reference build hg19 of the human genome using a banded Smith–Waterman algorithm from Illumina. In all steps, standard parameters were used except for the exclusion of duplicated read as our approach was based on amplicons. Deep resequencing for the 7 genes resulted in a population mean sequence coverage of 262× ranging from 93× for the sample with lowest coverage to 646× for the highest coverage (). Of the 42,265 exonic bases sequenced, 18,354 were in protein coding region. Focusing on coding sequences 92.6% of the bases were covered with more than 30× ().

Alignment and Variant Calling.

Haplotype caller was used to create individual GVCF files that were combined to perform joint genotyping using the GATK pipeline (45). Variants were hard-filtered using the suggested GATK thresholds for amplicon sequencing employing the following parameters for single-nucleotide polymorphisms (SNPs) (QD < 2, MQ < 40, SOR > 3, MQRankSum < −12.5 and ReadPosRankSum < −8.0) and for INDELs (QD < 2, SOR > 10, ReadPosRankSum < −20.0 and InbreedingCoeff < −0.8). Additional thresholds were used to ensure the quality of the variant calling by 1) setting as missing genotypes with GQ < 30 and/or with depth of coverage per individual <10× and 2) excluding variants with a population mean depth of coverage < 30× and call rates < 90%. The function of the variants that passed quality control were annotated with ANNOVAR (46). Transcribed variants (untranslated regions, splice sites, and synonymous and nonsynonymous variants) and the subset focusing only on coding nonsynonymous variants (missense, stop gain/loss, and frameshifts) were selected for further evaluation.

Statistical Analyses.

The enrichment of protein altering variants per gene was evaluated using 2 different statistical tests: the variable threshold (VT) (47) method and the optimal unified burden and SNP-set Kernel Association Test (SKAT-O) (48) both implemented in the software EPACTS (49). Briefly, VT is a burden test which collapses rare variants in a gene into a single burden variable and tests for the cumulative effect of rare variants in the gene by regressing the burden variable on the phenotype. Burden tests are the most powerful rare variant association tests when most variants are causal and the effects are in the same direction. However, in the presence of a mix of protective and deleterious variants in a gene, the nonburden SKAT is more powerful. To maximize power under both scenarios, SKAT-O adaptively combines a weighted burden and the SKAT statistics. SKAT-O with default parameters was used to weight variants according to their MAF adding the optional correction for studies with fewer than 2,000 samples. These approaches can be applied to both rare and low-frequency variants. Bonferroni was used for multiple testing correction by considering 2 tests (SKAT-O and VT) in 7 genes and 2 conditions (transcribed variants and protein altering) resulting in a cutoff of P < 0.002 to be significant at the genewise level. To evaluate if the effect per gene was driven by a particular variant, a univariable analysis was performed for nonsynonymous polymorphisms with MAF > 0.01 under additive model in PLINK (50). For the genes associated with T1R in the genewise analysis that also presented a variant with P < 0.1 in the univariate model a conditional genewise association analysis was performed with EPACTS.

A hypergeometric test was used to determine if T1R and PD cases shared substitutions at the same amino acid residues for the LRRK2 and Parkin proteins more often than expected by chance. As not all amino acids of a protein have the same probability of being mutated in the general population, 2 databases, dbSNP (51) and gnomAD (52), were accessed to catalog both LRRK2 and Parkin mutated residuals. In the enrichment analysis, only protein mutated residuals with frequency lower than 1% in the curated databases were considered as baseline since common mutations are likely to be present in both T1R and PD samples merely due to their frequency. For instance, in the curated databases Parkin presented 287 out of its 465 amino acids residuals mutated at least once in the general population. The number of Parkin residues associated with PD in the PDmutDB (20) and LOVD (21) databases was 72. The hypergeometric test calculated the statistical significance of randomly selecting 5 PD-associated residues in the T1R samples when 6 out of the 287 known mutated positions were observed in our population (). The same test was performed including novel variants for PRKN and applied to the LRRK2 gene.

ROS Detection.

Different cell lines were seeded in 96-well plates at a concentration of 3 × 104 cells per well and stimulated with interferon-γ (100 ng/mL) for 24 h. Cells were then infected with bacillus Calmette–Guérin Russia or M. leprae at a multiplicity of infection (MOI) of 10:1. At indicated time points following infection, ROS production was detected using ROS-ID total ROS detection kit (Enzo Life Science) according to the manufacturer’s instructions. Briefly, cells were carefully washed with 200 μL per well of 1× wash buffer. Following wash buffer removal, 100 μL per well of ROS detection mix (4 μL of 5 mM oxidative stress detection reagent/10 mL of 1× wash buffer) was added before incubation of plates in a humidified incubator (37 °C, 5% CO2) for 30 min. Readings were acquired at wavelengths 488/520 nm on a plate reader. Each experiment was performed at least 3 times, each in triplicate.

Apoptosis Assay.

One the day before infection, RAW264.7 cells expressing LRRK2 wild-type or mutant protein or depleted for LRRK2 gene expression by knockout were separately cultured at a concentration of 6 × 105 per well in a 6-well plate for 16 to 18 h. The cells were then infected with or without bacillus Calmette–Guérin (MOI 10:1) for 24 h followed by apoptosis analysis using Annexin V-APC and Zombie Aqua fixable viability dye as detailed in the Biolegend manual. Briefly, cells were washed twice with cold PBS and incubated with 1 mL of 5 mM EDTA/PBS on ice for 30 min followed by pipetting to detach and dissociate the cells. Cells were washed twice with PBS and then stained with Zombie Aqua in PBS for 30 min at 2 to 8 °C, avoiding light. After incubation cells were washed twice with PBS/0.2% BSA and once with 1× binding buffer followed by incubating with Annexin V in 1× binding buffer for 15 min at room temperature in dark conditions. Finally, cells were washed once with 1× binding buffer and resuspended in 200 μL of 1× binding buffer for flow cytometry with BD FACSCanto II (BD Biosciences). The data were analyzed on FlowJo v10.4.2 (FlowJo, LLC) with viability and Annexin V single stains as fluorescence minus one (FMO) controls.

Thermal Shift Assay.

Parkin at 10 μM (0.5 mg/mL) was mixed with 2× Sypro Orange (5,000× stock, S6650; ThermoFisher). The mix was distributed in a 96-well black-coated plate and analyzed with a QuantStudio 7 Flex analyzer (Life Technologies). The temperature was ramped from 15 °C to 70 °C with continuous fluorescence readout at excitation/emission 488/570 nm. The melting temperature was derived from the peak of the first derivative from the temperature–fluorescence profile. Measurements were performed in quadruplicates and analyzed using one-way ANOVA with the Dunnet’s test to calculate significance of difference to the wild-type group.

Intact Mass Spectrometry.

Protein samples were diluted at 0.5 mg/mL in 1% formic acid and 2 μL (1 μg) was injected on a Waters C4 BEH 1.0/10-mm column and washed 5 min with 4% acetonitrile, followed by a 10-min 4 to 90% gradient of acetonitrile in 0.1% formic acid, with a flow rate of 40 μL/min. The eluate was analyzed on a Bruker Impact II Q-TOF mass spectrometer equipped with an Apollo II ion funnel electrospray ionization source. Data were acquired in positive-ion profile mode, with a capillary voltage of 4,500 V and dry nitrogen heated at 200 °C. Spectra were analyzed using the software DataAnalysis (Bruker). The total ion chromatogram was used to determine where the protein eluted, and spectra were summed over the entire elution peak. The multiply charged ion species were deconvoluted at 5,000 resolution using the maximum entropy method.

Mitophagy Assay.

Mitophagy was examined using a fluorescence-activated cell sorting (FACS)-based analysis of mitochondrially targeted Keima (mt-Keima) as previously described (36, 53). Briefly, U2OS cells stably expressing an ecdysone-inducible mt-Keima were induced with 10 μM ponasterone A, transfected with GFP-Parkin for 24 h, and treated with 20 μM carbonyl cyanide m-chlorophenyl hydrazone (CCCP) for 4 h. For FACS analysis, cells were trypsinized, washed, and resuspended in PBS before their analysis on an LSR Fortessa (BD Bioscience) equipped with 405- and 561-nm lasers and 610/20 filters (Department of Microbiology and Immunology Flow Cytometry Facility, McGill University). Measurement of lysosomal mt-Keima was made using a dual excitation ratiometric pH measurement where pH 7 was detected through the excitation at 405 nm and pH 4 at 561 nm. For each sample, 50,000 events were collected, and single, GFP-Parkin–positive cells were subsequently gated for mt-Keima and the percentage of cells with an increase in the 405 nm:561 nm ratio of mt-Keima was quantified. Untreated, single GFP-Parkin–positive cells were gated for quantification of the geometric mean of the GFP signal as a measure of steady-state Parkin protein levels. Data were analyzed using FlowJo v10.1 (Tree Star). For statistical analysis, the data represent the average percentage of mitophagy from at least 3 independent experiments, and P values were determined by one-way ANOVA with Dunnett’s post hoc tests. *P < 0.05.