Introduction

Parkinson’s disease (PD) is a complex multifactorial neurodegenerative disease that affects millions. Dopaminergic neuron (DAN) death in the Substantia nigra causes motor defects in patients. However, 97 % of patients suffer from non-motoric dysfunction decades earlier, including constipation, REM sleep disorder and hyposmia (Ansari & Johnson, 1975; Chase & Markopoulou, 2020; Doty et al., 1988; Parkinson, 2002; Schapira et al., 2017). These non-motor defects in PD are often overlooked because they occur before disease diagnosis and thus, they do not belong to the normal clinical work-up. Furthermore, non-motoric defects are largely non-responsive to dopamine replacement therapy suggesting other cell types are involved (Kalia & Lang, 2015).

Hyposmia, the reduced ability to discern odors, is one of the earliest defects associated with familial and sporadic Parkinsonism (Chase & Markopoulou, 2020; Doty, 2012; Doty et al., 1988; Haehner et al., 2009). These patients show difficulties to classify odors: they often perceive banana odor as motor oil or beer odor as fruit in standardized odor identification tests (UPSIT) (Double et al., 2003). These observations suggest that hyposmia in Parkinsonism is caused by a defect in higher-order perception, consistent with the observations that first-order olfactory receptor cells of the nasal epithelium are not majorly affected in patients (Mueller et al., 2005; Witt et al., 2009). The synaptic site of second-order olfactory neurons, the olfactory bulb, is modulated by basal forebrain cholinergic neurons and local GABAergic and dopaminergic interneurons (Harvey & Heinbockel, 2018). While these DAN do not degenerate in Parkinsonism, it has been suggested there is local impairment of cholinergic transmission in the olfactory bulb of Parkinsonism patients and vertebrate models (Huisman et al., 2004; Mundiñano et al., 2013; S. Zhang et al., 2015), but the molecular defects are elusive.

There are only few models to study hyposmia in Parkinsonism and most are based on the use of toxins or expression of mutant alpha-Synuclein (Taguchi et al., 2020). While this may recapitulate aspects of the disease, they bias the disease mechanisms to known phenomena associated with the specific model. Modelling approaches to Parkinsonism are further hindered due to the large number of mutated genes. These mutations affect various biological functions and therefore it is currently hard to unveil the specific and possibly common pathways and cell types that are at the basis of the early defects in the disease. So far, mutations in more than 25 genes have been identified that cause familial forms of Parkinsonism (Brooker et al., 2024). These genes cause defects in a diverse range of processes including mitochondrial function, endocytosis, lysosomal function, autophagy and protein translation (Kalia & Lang, 2015). Nonetheless, these mutations do result in overlapping motor and non-motor defects, including hyposmia. This raises the question of whether defective mechanisms across the genetic space of Parkinsonism converge on common biological processes and affect similar cell types beyond DAN. While we are in the process of generating a comprehensive collection of fruit flies with PD knock-out and knock-in mutations (see also (Kaempf et al., 2024)), this study focuses on the first five Drosophila Parkinsonism models. Using these models, we investigated cellular and molecular dysfunctions that precede motor defects and find that early cholinergic projection neuron problems contribute to dopaminergic system failure later in life. The mutants included in this study are “classical” Parkinson disease genes, like LRRK2 and PINK1, and Parkinsonism genes that affect vesicle trafficking processes, including SYNJ1, DNAJC6/Aux and RAB39B. Additionally, we also analyzed postmortem brain samples from idiopathic patients with LRRK2 variants and iPSC derived human dopaminergic neurons with a pathogenic LRRK2 mutation.

Results

A new collection of five familial PD knock-in models

The characterization of cells and pathways affected in PD is hindered because patient populations are heterogeneous and because mutations in numerous genes have been linked to familial forms of the disease (Brooker et al., 2024). To overcome these problems, we created a new collection of Drosophila PD knock-in models carrying pathogenic mutations using a genetic standardized and controlled design, enabling direct side-by-side comparative analyses. Additionally, these were extensively backcrossed into one isogenic background (>10 generations) and did not rely on overexpression of proteins (e.g. α-Synuclein). We selected homologues of 5 PD genes and replaced the Drosophila lrrk, rab39, auxilin (aux), synaptojanin (synj) and pink1 genes by wild-type and pathogenic mutant human or Drosophila cDNA sequences (CDS) at the endogenous locus (Fig. 1a; wild-type hPINK1 was a UAS construct; see methods). We confirm that the wild-type and pathogenic mutant human CDS (LRRK2, RAB39B and PINK1) or Drosophila CDS (DNAJC6 (Aux) and SYNJ1 (Synj)) are expressed at endogenous levels (Supplemental Fig. 1a, Supplemental File 1). To determine if our new models recapitulate known PD-relevant phenotypes, we tested them in an array of assays that were previously used to analyze other fly PD models. Wild-type knock-ins are similar to controls (background fly line), indicating the cDNA knock-ins (human or fly) recapitulate normal gene function. Conversely, the PD knock-in models show defects akin to those seen in previously described mutants of these genes. For example, we find disturbances in mitochondrial membrane potential and electrophysiological defects in electroretinogram (ERG) recordings of flies stressed by exposure to constant light for 7 days (Hindle et al., 2013; Morais et al., 2009; Mortiboys et al., 2010; Ng et al., 2012; Vanhauwaert et al., 2017) (Fig. 1b, Supplemental Fig. 1b-d). Note that this specific experimental ERG setup differs from studies involving the progressive aging of PD mutant flies (Jacquemyn et al., 2023). Importantly, we also tested the flies for their performance in a startle-induced negative geotaxis assay (SING). This is a motor behavior that depends on central brain dopaminergic neurons (posterior-anterior-medial dopaminergic neurons (DAN): “PAM”). Previous work showed that PAM-dependent SING was unaffected in young flies that express α-Synuclein, but that SING declines as these flies age (Riemensperger et al., 2013). Similarly, none of our models shows a SING defect at a young age (5 d, young) (Fig.1b). However, as flies get older (25 d, old), the PD knock-ins display impaired performance in this DAN-dependent locomotion assay, while the performance of controls is not affected (Fig. 1b). This is in line with recent work where we show that SING defects in PD mutant models are rescued when the flies are fed L-DOPA but not D-DOPA indicating a very strong correlation between SING defects and defects in dopaminergic synaptic innervation of PAM DAN onto Mushroom body neurons (Kaempf et al., 2024). Taken together, our data suggests that the mutants we used in this study suffer from a progressive locomotion defect that is linked to DAN synapse dysfunction.

Fig. 1.

Cholinergic olfactory projection neurons are impaired across 5 PD mutant fly models

To define cell types affected by pathogenic PD mutations early in life, we used an unbiased approach and conducted single cell RNA sequencing of entire brains of young (5 d) animals (Supplemental File 2). We pooled 108k newly sequenced cells from these 5 mutants and controls with the 118k wild-type cells from our original atlas leading to 186 cell clusters, of which 81 were annotated to known cell types (Davie et al., 2018; Janssens et al., 2022) (Fig. 2a). PD mutant cells equally mix with the cells of the original fly brain cell atlas and the frequency of recovered cell types is similar between control and mutant brains (including DAN; Fig. 2b-f, Supplemental Fig. 2a, b; Supplemental File 3. We further do not find any cell type that shows higher or lower expression of the knocked-in CDS (Supplemental Fig. 2c). Thus, cellular identity and cellular composition are preserved in young PD fly models.

Fig. 2.

For every cluster, we then compared gene expression between cells from mutants and controls using two independent algorithms (Wilcoxon and DESEQ2, (Wang et al., 2019)) and find both methods correlate well (average Spearman correlation of signed p-value > 0.8 to > 0.95; Supplemental Fig. 2d). However, the number of differentially expressed genes (DEG) is largely dependent on the number of cells present in a cluster. To correct this bias, we fitted a negative binomial model that correlated the number of DEG and cell cluster size in each mutant (Supplemental Fig. 2e, Supplemental File 4). Using the residuals of this model we find eight common predominantly affected cell types across the 5 mutants (Supplemental File 5 and 6). We also find cell types uniquely affected in some mutants, e.g., glial subtypes that were affected in PINK1^P399L. The commonly affected cell types include cholinergic olfactory projection neurons (OPN), Kenyon cells, cholinergic neurons of the visual system (T1 neurons), and other yet-to-be-identified cell types (Fig. 2b-f, black cells; Supplemental File 6). Hence, while our mutants are modelling different genetic origins of disease, we find some cell types are commonly affected. While several cell types are affected (and many are cholinergic), OPN are here of particular interest because they are higher-order projection neurons that control innate olfactory processing and odor classification (Masse et al., 2009).

Synaptic genes and pathways are deregulated in cholinergic brain regions of fly PD models and human PD patients

We next compared the identity of the DEG found in PD fly model cholinergic olfactory projection neurons to the identity of the (homologous) genes differentially expressed in patient brain samples rich in cholinergic neurons. We analyzed nucleus basalis of Meynert (NBM), nucleus accumbens (Nac) and putamen (Put) from five brains of PD patients (all have idiopathic LRRK2 risk mutations) and five brains from unaffected age-matched individuals (without PD-relevant genetic variants). Our analysis identified 37 distinct cell types in these brain regions (Supplemental Fig. 3a), and most of the cell types have a similar frequency when comparing PD vs non-PD condition. An exception is a microglia subtype (MG #20, Supplemental Fig. 3b), which is 900% increased in frequency in PD patient brains compared to non-PD samples. We next analyzed all neuronal subtypes in these samples (Supplemental Fig. 3c, d) and identified the DEG. Many of the DEG are genes encoding synaptic proteins that are also listed in the synapse-specific portal SynGo (Koopmans et al., 2019) (Fig. 3a, a’; Supplemental File 7). We then asked how they compare to the DEG found in affected cholinergic projection neurons of our five fly models (Fig. 3a”; Supplemental File 7). We therefore first listed the deregulated genes in affected cell types for each fly model. There are many DEG unique to each PD mutant, but there is also overlap in DEG between the models. To define this overlap across fly PD mutants, we ranked genes according to their signed differential expression for each of the 5 PD models (up or downregulated) and then summed their rank across the mutants retaining the genes in the top 5% and the bottom 5% as the common highly up- or downregulated genes (Supplemental File 7); i.e. only if the DEGs are mostly commonly up- or commonly down-regulated they end up at the top or bottom of our list. Additionally, when we compare this list of deregulated fly genes to the DEG-orthologues in PD patients, there is remarkable overlap. For example, in the NBM, an area associated with PD (Arendt et al., 1983), >20% of the neuronal DEG that have an orthologous gene in the fly are also found among the most deregulated genes across PD fly models. This is highly significant: of the 2486 significantly differentially expressed human genes, 1149 have a fly orthologue, and of these, 28.46% overlap with the deregulated fly genes (5% top and bottom genes listed in Supplemental File 7). Performing a hypergeometric test confirms that this overlap is significant (p, 9.06e-⁷⁶). Next, we analyzed the list of human-fly overlapping DEG using Gorilla gene ontology enrichment (Eden et al., 2009). Markedly, this shows that the processes most strikingly enriched in both fly models and human NBM are those related to synaptic function (Fig. 3b). Indeed, more than half of the DEG found in human and flies are present in the synapse-specific portal SynGo (Koopmans et al., 2019) (Fig. 3a”; Supplemental File 7). We found enrichment for presynaptic function, including synaptic exocytosis, axon guidance, synapse organization, ion transport and protein homeostasis (Fig. 3b, c). Additionally, we also found enrichment of genes involved in RNA regulation and mitochondrial function that are also important for the functioning of synaptic terminals (Supplemental File 7) (Gorenberg & Chandra, 2017; Morais et al., 2009; Snead & Eliezer, 2019; Uytterhoeven et al., 2011; Verstreken et al., 2005). Of note, while the human samples all have LRRK2 variant mutations, comparing the vulnerable gene signatures from each of the PD fly models to the DEGs from the human samples does not show any greater similarity between the LRRK2 mutants compared to the other PD mutants (Supplemental File 7). In summary, cholinergic olfactory neurons of young familial fly PD models and human PD patients both show preferential transcriptional deregulation of genes that regulate synaptic function.

Fig. 3.

Young PD models suffer from synaptic impairment of cholinergic OPN

The transcriptional deregulation of genes encoding synaptic proteins suggests synaptic alterations in the affected cell types. We therefore expressed the Ca²⁺-sensor GCaMP3 in the OPN, one of the majorly affected cell types in our PD flies, and monitored synaptic terminal function through a small window in the head capsule of living flies upon stimulation (Fig 4a, b). Neuronal stimulation results in a robust synaptic Ca²⁺ signal at OPN synapses in controls (Fig. 4c, d (black points), quantified in Fig. 4g (gray bar)) and in animals with the wild-type PD gene knocked-in (e.g. hLRRK2 in Fig. 4c, d (dark red points), quantified in Fig. 4g (left bar)). In contrast, the five PD mutant knock-ins show significantly weaker Ca²⁺-responses (e.g. hLRRK2^G2019S in Fig. 4c, d (pink points), and quantified in Fig. 4g (middle bar; “OPN>wt CDS -”)). Next, we tested if this defect is cell-autonomous and created PD mutant knock-ins expressing the wild-type CDS in their OPN using GH146-Gal4 (Fig. 4g (right bar; “OPN>wt CDS +”)). Even though this Gal4 driver also drives expression in one pair of inhibitory APL neurons (GABAergic neurons that we did not recover in our single cell sequencing), it mostly expresses in OPN enabling us to assess PD gene function in OPN (Li et al., 2017; Liu & Davis, 2009). PD mutants expressing wild-type PD genes using GH146-Gal4 show synaptic Ca²⁺ signals similar to controls, indicating the PD genes are required in OPN of these young animals to maintain robust synaptic function (Fig. 4g (right bar; “OPN>wt CDS +”). Note that while hLRRK2 wild-type knock-in flies are very similar to wild-type controls (Fig. 4g (right bar)), the hLRRK2^G2029S knock-in flies are not rescued by OPN-specific expression of hLRRK2 (Fig. 4g (right bar); “OPN>wt CDS +”). This is in agreement with the G2019S mutation being dominant (West et al., 2005; Zimprich et al., 2004).

Fig. 4.

Previous work indicated that functional defects of OPN often cause morphological changes at the level of the synaptic area in the calyx (Kremer et al., 2010). We thus also assessed the synaptic area of OPN in the calyx (based on the GFP fluorescence area of the GCaMP3 sensor). Interestingly, the PD mutants have a significantly larger area compared to controls, despite very similar GFP expression levels (Fig. 4e, f, h). Again, this defect is cell-autonomous, because it is rescued by OPN-selective expression of the wild-type CDS, except -again- for LRRK2^G2019S agreeing with it being a dominant mutant (Fig. 4h). Thus, OPN of young PD mutants display cell-autonomous presynaptic defects: even though the synaptic area is increased, synaptic Ca²⁺ signals are diminished.

Hyposmia is prevalent in young PD models

While old (25 d) PD mutants show PAM-dependent locomotion dysfunction in the SING assay (Fig. 1b (SING old)), young (5 d) PD mutants do not show SING defects (Fig. 1b (SING young)). These data indicate PAM neurons are functional in young PD mutants. Yet, these young mutants suffer from OPN-synaptic dysfunction (Fig. 4). We thus wondered if our young PD models have olfactory performance defects. We subjected flies to a choice between two odors (test 1: motor oil and banana or test 2: beer and wine). Control and wild-type knock-in flies show a robust odor preference in these two odor tests (Fig. 5a, b). Conversely, the knock-in flies carrying the pathogenic mutations behave indecisive (Fig. 5a, b). Hence, young PD mutants show olfactory behavior defects.

Fig. 5.

We then asked if this defect is cell-autonomous and expressed the relevant wild-type CDS selectively in the OPN of the pathogenic knock-in mutant flies. This manipulation restores olfactory performance back to control levels (Fig. 5c, d; “OPN>wt CDS +”, green label) when compared to PD mutants lacking OPN expression (Fig. 5c, d; “-”). However, the LRRK2^G2019S mutant remains unaffected, consistent with its classification as a dominant mutant. Note that lrrk^KOflies also show an olfactory performance defect (Fig. 5c, d; “-”) and that expressing hLRRK2 cDNA specifically in the OPN of lrrk^KO rescues this defect (Fig. 5c, d; “OPN>hLRRK2 +”, green label). This indicates hLRRK2 is a functional protein that, in this context, can compensate for the function of the fly orthologue. Finally, we demonstrate that the rescue of the olfactory performance in PD mutants is specific for OPN, because expressing wild-type CDS in another affected cell type, “T1-interneurons of the visual system”, of knock-in PD animals (we chose to test this in Rab39B^G192R mutants), does not improve olfactory performance (Fig. 5c, d; “T1>hRab39B +”, red label). Hence, a cell autonomous defect in cholinergic OPN causes olfactory defects across PD mutants.

Rescuing OPN dysfunction in young PD mutants prevents DAN defects at older age

Our data indicates that young PD mutants (aux, synj and LRRK2) have olfactory defects while displaying normal DAN function (SING). Conversely, aged aux^R927G, synj^R258Q and hLRRK2^G2019S knock-in animals exhibit both olfactory preference defects and progressive DAN (SING) dysfunction (Fig 6a; line graphs; “DAN-dependent movement behavior”). We also assessed the synaptic afferent innervation area of PAM DAN onto mushroom bodies by quantifying the area of anti-TH labeling. In line with the SING behavioral defect, all three 25-day-old knock-in mutants show a significant reduction in DAN afferent (Fig. 6a, left bar; “OPN>wt CDS -”; “DAN synapse morphology”). We then wondered if rescuing OPN at young age affects DAN at later age. We therefore expressed the wild-type PD gene in (a.o.) OPN of aux^R927Gand synj^R258Q mutants using the GH146-Gal4 (that does not drive expression in DAN) (Fig. 6a, right bar; “OPN>wt CDS +”; “DAN synapse morphology”). Given that hLRRK2^G2019S is dominant, we did not express the wild-type hLRRK2 in OPN, but endoA^S75D that we showed previously, genetically interacts with lrrk (Matta et al., 2012) (Fig. 6a, right bar; “OPN>endoA^S75D +”; “DAN synapse morphology”). We find these genetic manipulations rescue the DAN-associated SING and synaptic connectivity defects in old PD mutants (Fig. 6a). Altogether, these results suggest that functional OPN have a cell non-autonomous effect to prevent DAN dysfunction that emerges at old age in the PD mutants. To further substantiate this observation, we also conducted a set of independent experiments where we fed hLRRK2^G20190S mutants nicotine. Nicotine activates the nicotinic acetylcholine receptors that, under normal circumstances, are also activated by the release of acetylcholine from cholinergic neurons, including OPN (Chambers et al., 2013; Changeux, 2010; Zhou et al., 2001). While feeding nicotine does not rescue the olfactory preference defect of hLRRK2^G2019S mutants, it also does not rescue the OPN synapse morphology defect or the OPN-associated defects in synaptic Ca-imaging (Fig. 6b). Interestingly, nicotine does rescue the DAN-associated defects, including SING, synapse loss and defects in Ca-imaging at DAN synapses (Fig. 6c). These results are consistent with the role of OPN function, specifically acetylcholine release, in preventing DAN damage in PD mutants.

Fig. 6.

In a final experiment we assessed if this effect of nicotine on DAN health is conserved across species. We therefore generated human induced neurons derived from iPSC in which we engineered a LRRK2^G2019S mutation and differentiated the cells into DAN (Supplemental Fig. 4) (Kriks et al., 2011; Nolbrant et al., 2017). Our protocol results in >50% DAN in both the LRRK2^G2019S mutant line and the isogenic control (Fig. 6d, Supplemental Fig. 4f). To assess neuronal activity, we expressed GCaMP6f and find that LRRK2^G2019S mutant DAN are significantly less active and the amplitudes of their Ca²⁺-spikes are smaller than those in isogenic controls (Fig. 6e). Next, we incubated the differentiated DAN with nicotine and found this rescues the LRRK2^G2019S-induced neuronal activity defects. This effect is specific to nicotine as no rescue was observed when cells were co-incubated with mecamylamine, a non-competitive antagonist of nicotinic acetylcholine receptors, trumping the effects of nicotine (Fig. 6e). Hence, the positive effect of nicotine on PD mutants is conserved between flies and iPSC-derived dopaminergic neurons.

Discussion

Parkinsonism is currently diagnosed based on motor defects that only occur following significant DAN degeneration. However, patients suffer from non-motoric problems decades earlier, including sleep problems and hyposmia (the inability to properly discern odors) (Yamakado & Takahashi, 2024). Most of these non-motoric problems are not responsive to dopamine replacement therapy, suggesting other transmitter systems are involved. Using unbiased methodology, we find cholinergic system failure as one of the earliest defects across 5 different Parkinsonism fly models and we show that synaptic failure of specific cholinergic projection neurons (OPN) causes non-motor problems, including hyposmia. Also in patients, cholinergic defects in several brain regions have been observed, including the olfactory bulb and the occipital cortex (Kuhl et al., 1996; Shimada et al., 2009). These defects correlate with several non-motor symptoms, including impaired olfactory and visual processing (Bohnen et al., 2010). We now place cholinergic failure firmly ahead of dopaminergic system failure in flies, and we suggest that dysfunction of specific cholinergic projection neurons may be a prodromal signal of disease. Our work opens the possibility to explore the early detection of specific cholinergic interneuron dysfunction, e.g. using PET tracing, as an entry point to move diagnosis and intervention ahead of dopaminergic degeneration.

The connection between early cholinergic neuron dysfunction and subsequent failure of the DAN system in disease progression remains uncertain. However, dopaminergic neurons have a high density of nicotinic acetylcholine receptors, and evidence suggests that stimulating these receptors enhances DAN survival (Brimblecombe et al., 2018). This is supported by epidemiological data indicating that nicotine exposure from chronic smoking protects against PD (Janssens et al., 2022; Wang et al., 2019), and various studies, including the work we present here, show that direct nicotine supplementation mitigates DAN dysfunction in mouse and fly PD models (Chambers et al., 2013; Fares et al., 2023; Olsen et al., 2023). We also demonstrate this with a genetics experiment where we express wild-type PD genes in OPN and observe rescue of DAN dysfunction in older flies. These findings imply that factors affecting the function of OPN or cholinergic neurons might, by absence of insufficient innervation or activity, lead to DAN problems and degeneration, warranting further exploration of the underlying molecular mechanisms.

Our findings reveal a common trigger relevant for both familial and idiopathic cases. Observations of Lewy body pathology in postmortem samples suggest the disease may nucleate at specific sites and then propagate. This includes spreading of the disease between the gut and the brain as well as early aggregation of alpha-Synuclein, a synaptic protein, in the olfactory bulb (Braak, Tredici, et al., 2003). The presence of protein aggregates in these areas is consistent with protein homeostasis and synaptic defects that ultimately contribute to the aggregation of alpha-Synuclein (Chu et al., 2009; Friedman et al., 2012; Sarkar et al., 2007; Webb et al., 2003). In idiopathic disease, there is some evidence this is because the gut and olfactory bulb synapses are more easily exposed to environmental toxins that cause protein misfolding and aggregation (Braak, Rüb, et al., 2003). Furthermore, we also find genetic interactions with the protein turnover machinery that warrant further investigation (this study and (Matta et al., 2012)). However, different mechanisms may be at play since our single-cell sequencing shows the deregulation of several synaptic pathways in the PD knock-in mutants. While several genes mutated in familial forms of PD have proposed synaptic functions (endocytosis, synaptic autophagy, etc.), others may be affecting the synapse indirectly (e.g. mitochondria, transport, etc.). Future work will now need to elucidate which synaptic processes are affected and how this has an impact on cholinergic projection neuron function and cholinergic transmission.

Strengths and caveats to our approach

Strengths

In this work, we generated new Drosophila PD models using a controlled and standardized design to introduce pathogenic mutant PD variants at the endogenous locus. With this approach, the gene regulatory environment is as much as possible preserved ensuring that our models are physiological relevant, circumventing overexpression artifacts or ectopic effects. By studying these mutations in their endogenous context, we can better understand their specific contributions to transcriptional and cellular phenotypes. Given that the mutants are generated in a controlled ‘isogenic’ fashion, it is now possible to conduct careful side-by-side comparisons of different PD mutant genotypes.

Analyzing these newly generated PD models at young age using single-cell RNA sequencing enables us to identify affected cell populations and uncover early molecular events that might precede overt neurodegenerative disease phenotypes. Furthermore, the comparison between our fly PD models and human patient brain samples provides evidence that some of the molecular mechanisms of disease are conserved across species.

Finally, our in vivo investigations enabled us to uncover cell-non-autonomous effects of cholinergic projection neurons onto the fly dopaminergic system, a feat that would have been difficult to discover in in vitro (mono-) human neuron cultures.

Caveats

One of the caveats in our study involves the use of the GH146-Gal4 line. GH146-dependent Gal expression includes OPNs (that are cholinergic) and one pair of inhibitory APL neurons (that are GABAergic). There are only 2 APL per fly brain and our single cell sequencing experiment does not have the resolution to allow us to test if these neurons had a significant number of DEG. However, we are able to rescue DAN dysfunction by mimicking cholinergic output (using nicotine). These data do not exclude that APL-neuron problems contribute to the defects we observe in our PD mutants, but they do suggest that cholinergic output is critical to maintain normal DAN function.

Another caveat involves the temporal mismatch of comparing single-cell RNA sequencing data obtained from young PD models versus postmortem end-stage human disease tissue. Indeed, there may be changes in early stages that may not completely align with the late-stage disease profiles typically observed in human postmortem samples. Unfortunately, obtaining human brain samples from PD patients at an early stage are rare, these samples would have been more relevant to be compared to our young PD models. Nonetheless, we do also find that cholinergic neurons from postmortem brain samples have similar affected DEG pointing to similar affected molecular pathways when comparing to our young PD fly data.

Methods

Resource availability

Lead contact

For further information and requests for resources and reagents contact Patrik Verstreken (patrik.verstreken@kuleuven.be).

Material availability

Plasmids and fly lines generated in this study are available on request.

Data and code availability

Codes and data sets generated in this study are stored online https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE228843 (full fly dataset) and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE235332 (full human dataset) with restricted access until publication. Access will be given to reviewers via token.

Experimental model details

The “control” we refer to is a semi-isogenized w¹¹¹⁸backcrossed to CantonS for 10 generations (cantonized w¹¹¹⁸). We sequenced the entire genome of the resulting fly line and defined there are no deleterious mutations in homologues of known causative Parkinson’s disease genes or in homologues of human genes located close to GWAS loci. All our knock-ins (see below), UAS (see below) and Gal4 lines were backcrossed for >10 generations to this control. All generated fly lines were verified by multiple PCRs and PCR-product sequencing.

Flies were reared on standard yeast food medium at 25 °C and kept on 12/12 h light/dark cycles. Flies were kept in mixed populations of analogous density and were flipped onto fresh food every three days. Brains for RNA experiments, including single cell sequencing, were dissected between 11 and 12 h Zeitgeber time, behavioral and physiology experiments were performed between 18 and 20 h Zeitgeber time. Males and females were tested together, quantified separately, and results pooled if there was no statistical difference. Male hPINK1 pathogenic mutant knock-in flies are sterile, requiring balancing of the stock and consequently they produce only male hemizygotes to be quantified in experiments.

Generation of Drosophila PD models

A complete list of primers, gRNAs, oligos, gBlocks and fly strains is listed in Supplemental File 1. To clone inserts into plasmids NEBuilder® HiFi DNA Assembly (NEB) was used.

Unique gRNAs for LRRK2, PINK1, Rab39B, Aux and Synj were identified by http://crispr.mit.edu/ and introduced into pCFD4: U6:1-gRNA U6:3-gRNA (Port et al., 2014) by following an established protocol: http://www.crisprflydesign.org/wp-content/uploads/2014/06/Cloning-with-pCFD4.pdf.

Homology arms of 1 kB surrounding the first common exon of the different transcripts from the gene of interest were cloned into pWhite-STAR (Choi et al., 2009). For Aux, homology arms surrounded the entire Aux gene. HiFi DNA Assembly (NEB) was performed with four fragments: mini-white IMCE-cassette, left homology arm, right homology arm and the plasmid backbone. The homology arms were either generated by gBlocks (IDT), which were modified in order to make them compatible with the requirements of IDT, or amplified by Q5 PCR from genomic fly DNA of the target genotype (cantonized w¹¹¹⁸). The mini-white IMCE-cassette and vector backbone were generated by restriction digest on pWhite-STAR (Choi et al., 2009) with AvrII and XhoI. In order to avoid cutting of the donor plasmid, the PAM sequence of the target was mutated, if located in the homology arms. Mostly NGG to NAA modifications were used, however different modifications were used when necessary (presence of PAM in coding region). The donor plasmids pDoC-miniwhite-gene (Donor cassette) contains an Integrase mediated exchange (IMCE) cassette surrounded by two homology arms that allows for HDR (homology directed repair). HDR occurs when two double strand breaks are introduced in the DNA surrounding the first (or second) common exon of all transcripts. This HDR replaced the first common exon by the IMCE-cassette. Site of HDR was directed by the choice of the homology arms in a way such that the two attP sites for IMCE were situated inside evolutionary non-conserved regions.

Rescue plasmids were used to exchange the mini-white IMCE-cassette with the coding DNA sequence (CDS) of target genes. All rescue plasmids contain: 1) a stretch of genomic sequence that will be referred to as ‘utr’: this sequence was deleted together with the removed exon of the target gene and now is being placed back to repair the locus; 2) CDS of the target gene and 3) a poly(A)-tail. These 3 elements are surrounded by two attB sites. The rescue plasmids are named pReC-gene or pReC2-gene and are target gene-specific. The backbone: pReC was generated by linearizing pUC19 with SapI and EcoRI followed by insertion of two gBlocks: attB-MCS-L and MCS-attB-R. We also created pReC2, an optimized version that allows for easier cloning of the utr and CDS as it already contains the SV40 poly(A), by linearizing pUC19 with SapI and EcoRI followed by insertion of two gBlocks: attB-2xMCS-L and 2xMCS-attB-R.

Note that while for LRRK2, RAB39B and PINK1 knock-in models we used human CDS, for aux and synj we resorted to Drosophila CDS, since knock-in of the human homologues at the locus of these genes did not rescue its loss-of-function (lethality).

LRRK2

pReC was linearized with SapI. gBlock LRRK2:utr-P2A-SapI-SV40, containing the utr a P2A site (since we knocked-out exon 2), SapI restriction site (for subsequent cloning of the CDS) and an SV40 poly(A)) was introduced. The resulting plasmid was linearized with SapI. Next the CDS of LRRK2 was amplified from 2XMyc-LRRK2-WT with F_hLRRK2_cDNA and R_hLRRK2_cDNA and introduced into the linearized plasmid. The resulting plasmid is called pReC-hLRRK2. pReC-hLRRK2G2019S and pReC-hLRRK2Y1699C were created by using Q5® Site-Directed Mutagenesis (NEB) with primers F_mut_LRRK2_G2019S, R_mut_LRRK2_G2019S and Fmut_LRRK2_Y1699C, R_mut_LRRK2_Y1699C respectively on pReC-hLRRK2.

PINK1

pReC was linearized with SapI. gBlock: PINK1:utr-SapI containing the utr and SapI restriction site (for subsequent cloning of the CDS) is cloned into pReC. The resulting plasmid was linearized with AflIII. Next, the gBlock PINK1:SV40poly(A) containing the human Sv40 poly(A) was introduced into the linearized plasmid. This plasmid is linearized with SapI. Next the CDS of human PINK1 was amplified from pLenti6-DEST PINK1-V5 WT with F_hPink1_cDNA and R_hPINK1_cDNA and introduced into the linearized plasmid. The resulting plasmid is called pReC-hPINK1. pReC-hPINK1R145Q and pReC-hPINK1L347P were created by using Q5® Site-Directed Mutagenesis (NEB) with primers F_mut_PINK1_L347P, R_mut_PINK1_L347P and F_mut_PINK1_P399L, R_mut_PINK1_P399L respectively on pReC-hPINK1. While we created wild-type CDS knock-ins for rescue experiments and knock-ins carrying the pathogenic mutant variants of all Parkinson’s disease models (including the Pink1 mutant knock-ins), we were unsuccessful in targeting the wild-type human PINK1 CDS to the fly pink1 locus, unlike the two pathogenic variants. Therefore, the Pink1 control we used was generated by neuronally expressing human PINK1 with UAS-hPINK1 and nsyb-Gal4 (Yang et al 2003) in the Pink1^KO background (full genotype: w¹¹¹⁸ Pink1^K.O.;UAS-hPINK1/+; nSyb-Gal4/+).

Synj

pReC was linearized with SapI. gBlock Synj_utr-SapI containing the utr and SapI restriction site (for subsequent cloning of the CDS) is cloned into pReC. The resulting plasmid was linearized with AflIII. gBlock Synj:SV40poly(A) containing the human Sv40 poly(A) was introduced into the linearized plasmid. This plasmid is linearized with SapI. Next, the CDS of Drosophila synj and synj^R258Qwere amplified from respectively UAS-Synj+ and UAS-SynjRQ (Vanhauwaert et al., 2017) with FW Synj and RC Synj and are introduced into the linearized plasmid by utilizing HiFi assembly. The resulting plasmids are called pReC-dSynj and pReC-dSynjR258Q.

Rab39B

pReC2 was linearized with XhoI and XbaI. A gblock: Rab39_utr was inserted into the pReC2. Next, this plasmid was linearized with SapI and a gBlock containing the human RAB39B CDS was cloned. The resulting plasmid is called pReC2-hRab39. pReC2-hRab39T168K and pReC2-hRab39G192R were created by using Q5® Site-Directed Mutagenesis (NEB) with primers F_mut_hRab39_T168K, R_mut_hRab39_T168K and F_mut_hRab39_G192R, R_mut_hRab39_G192R respectively on pReC2-hRab39.

Aux

The entire Drosophila auxilin gene region was cloned from BAC CH322-22D05 into pRec using BbsI, with an HA-tag at the N-terminal of the gene. AuxR927G was produced by Q5® Site-Directed Mutagenesis (NEB). See also (Jacquemyn et al., 2023).

UAS-aux

pUAST.attB was linearized with XhoI and XbaI; gBlock: Aux was cloned into this plasmid. See also (Jacquemyn et al., 2023).

Constructs were injected in house or by BestGene (BestGene Inc, CA, US).

Quantitative RT-PCR

For each sample, 20 heads of 5 day (5 d) old males were collected, RNA was extracted using Maxwell RSC instrument and Maxwell RSC RNA kit (Promega) and transcribed using the SuperScriptIII synthesis system (ThermoFisher Scientific). Subsequently, cDNAs were quantified by qPCR using a LightCycler 480, 480 SYBR Green master mix (Roche) using qPCR primers enlisted in Supplemental File 1. mRNA levels were determined using the Δ-Δ-CT method, where Ct values were first normalized to the housekeeping gene Rp49, and next expressed as a percent of endogenous Drosophila gene expression.

Mitochondrial membrane potential (TMRE) assay

Measurements of mitochondrial membrane potential were performed with the potentiometric dye tetramethyl rhodamine ethyl ester (TMRE). Male third-instar larvae were dissected in HL3 (in mM: 110 NaCl, 5 KCl, 10 MgCl₂ • 6H2O, 10 NaHCO₃, 30 sucrose, 5 trehalose, and 10 HEPES, pH 7.2). Larval fillets were incubated for 15 min in the presence of 50 nM TMRE (Abcam). Subsequently, the TMRE solution was discarded and fillets were rinsed three times with normal HL3 without TMRE. Mitochondrial labeling of TMRE was imaged on a Nikon spinning disk confocal microscope with a 40x water dipping objective 0.8 NA, excitation wavelength of 561 nm and an emission filter 595/50 nm. We imaged at 0.5 µm Z-steps entire NMJ branches. Larvae were stained a mutant and control larva were processed and imaged side-by-side. TMRE labeling intensity was quantified in individual synaptic boutons using Fiji (Schindelin et al., 2012). Because of variability between sessions, every mutant was normalized to the controls stained and measured in the same session.

Electroretinograms (ERG)

ERGs were recorded as previously described (Slabbaert et al., 2016). Flies were immobilized on glass microscope slides, by use of liquid Pritt glue. For recordings, glass electrodes (borosilicate, 1.5 mm outer diameter) filled with 3 M NaCl were placed in the thorax as a reference and on the fly eye for recordings. Responses to repetitive 1 s light stimuli were recorded using Axosope 10.7 and analyzed using Clampfit 10.7 software (Molecular Devices) and Igor Pro 6.37. To assess the ERG response after exposure to constant light, 3 day old flies were placed under continuous illumination (1300 lux) for 7 days at 25 °C prior to ERG data acquisition as described (Soukup et al., 2016; Vanhauwaert et al., 2017).

Startle-induced negative geotaxis

Geotaxis locomotion was quantified in mixed populations of ∼20 flies (Benzer, 1967) with modifications described by (Inagaki et al., 2010). Flies were flipped into the apparatus without anesthesia and allowed to adjust for some minutes to experimental environment (24 °C, 50 % humidity). They were tapped down into the first lower tube and allowed to climb into the upper tube for 30 s, at which time the flies that reached the upper tube were moved to tube 2. The procedure was repeated five times and the SING scored: (#flies in tube 1 + (#flies in tube 2 * 2) + (#flies in tube 3 * 3) + (#flies in tube 4 *4) + (#flies in tube 5 * 5)) / (#total flies *5).

Modeling to define the number of cells to be sequenced

We randomly sampled between 0 and 50,000 cells, in steps of 100, from the annotated brain cell dataset generated in (Davie et al., 2018). For each sample we calculated the number of cell types that were retrieved for different thresholds of detection. This process was repeated 50 times for convergence and Michaelis-Menten kinetics were fit through the data points.

Brain dissection for single cell sequencing

30 fly brains (15 females and 15 males, except Pink1 mutants, where only males were used) were dissected on ice and collected in 100 µl ice cold DPBS solution, centrifuged at 800 x g for 5 min, and placed into 50 μL of dispase (3 mg/mL, Sigma-Aldrich, D4818-2mg) and 75 µl collagenase I (100 mg/mL, Invitrogen, 17100-017). Brains were dissociated in a Thermoshaker (Grant Bio PCMT) for 2 h at 25 °C and 500 rpm. The enzymatic reaction was reinforced by pipetting every 15 min. Subsequently, cells were washed with 1000 µl ice cold DPBS solution and resuspended in 400 µl DPBS 0.04 % BSA. Cell suspensions were passed through a 10 μM pluriStrainer (ImTec Diagnostics, 435001050), cell viability and concentration were assessed by the LUNA-FL™ Dual Fluorescence Cell Counter, and cells were immediately used for the 10x run. Samples from the different genotypes were separately prepared and sequenced, but they were all processed in parallel to avoid batch effects. The control line (cantonized w¹¹¹⁸) and two extra wild types strains (DGRP-551 and w¹¹¹⁸) were included.

10x Genomics high-throughput sequencing

Before sequencing, the fragment size of every library was analyzed on a Bioanalyzer high sensitivity chip. The libraries were diluted to 2 nM and quantified by qPCR using primers against p5-p7 sequence. 10x libraries were sequenced twice: shallow on NextSeq500 (Illumina) to determine library quality and more in depth using NovaSeq6000 (Illumina). The targeted saturation for each sample was around 60 %, and additional sequencing runs were performed if needed. Each sequencing run used the following sequencing parameters: 28 bp read 1 – 8 bp index 1 (i7) – 91 bp read 2.

10x Data Preprocessing

The 10x fly brain samples were each processed (alignment, barcode assignment and UMI counting) with Cell Ranger (version 3.0.2) count pipeline. The Cell Ranger reference index was built upon the 3^rd 2017 FlyBase release (D. melanogaster r6.16) (Gramates et al., 2017). To count the human mutated genes, sequences from the construct were added as artificial chromosomes to the index. For every sample all sequencing runs were aggregated using CellRanger and the – recovered-cells parameter was specified as 5000.

Data filtering and clustering

CellRanger generated a gene expression matrix, which was used as input for VSN pipelines (Zenodo, v0.27.0) (Flerin et al., 2021) for the first filtering and data clean up. These pipelines are written in Nextflow DSL2 making them highly reproducible, and the config file used for the fly samples can be found in Supplemental File 2. The VSN pipeline removed predicted doublets using scrublet (Wolock et al., 2019), followed by a cleaning of ambient RNA reads using DecontX (S. Yang et al., 2020), similar to the procedures used in the Fly Cell Atlas (Li et al., 2022).

Afterwards the data was loaded into scanpy (1.8.2, (Wolf et al., 2018)). Cells were filtered using min_genes=200, max_genes=7000, min_counts=500 and max_counts=30000, and finally with percentage of mitochondrial counts below 15. The data were log-normalized with a scale factor of 10⁴ and scaled (maximum values = 10), regressing out the number of UMIs and the percentage of mitochondrial genes as latent variables using a linear regression model. Highly variable genes were selected using sc.pp.highly_variable_genes with default parameters.

We used sc.tl.pca to calculate 250 principal components on the highly variable genes, from which the first 100 components were selected for follow-up analysis based on an elbow plot. Since the data contains many samples (n=65) covering a variety of different ages and genotypes, we noticed batch effects and low mixing with default parameters. Therefore, we used Harmony (Korsunsky et al., 2019) to correct the principal components. For this we used sce.pp.harmony_integrate with max_iter_harmony=20, using the genotype + age combination as batch variable. The batch effect-corrected components were then used as input to calculate neighbors using sc.pp.neighbors with n_neighbors=30, followed by Leiden clustering for resolutions between 2 and 20. Subsequently, tSNE and UMAP visualizations were calculated.

Label transfer

In the previous step different resolutions were used in clustering, leading to an increasing number of clusters detected. To start annotating these and to find a final clustering resolution, we compared the detected clusters with the annotations from the co-clustered published dataset (Davie et al., 2018) and annotations from the Fly Cell Atlas (Li et al., 2022). In addition, we also used the neural net classifier (Davie et al., 2018) using default parameters as described. Finally, we compared the clusters with sorted cell types from (Davis et al., 2020) using lasso regression as described below. These data were used to annotate clusters in resolution 20, after which clusters with the same annotation were merged. For unknown cells, we chose to use Leiden resolution 5 as basis annotation to not split clusters too much and keep enough cells to perform meaningful differential expression analysis.

In total we found 188 clusters, of which 83 were annotated (median size = 1,011 cells, total = 140,778 cells) and 105 unknown (median size = 810 cells, total = 145,010 cells).

Lasso regression

To expand the published atlas, we analyzed the bulk transcriptomes of multiple FAC-sorted neuronal cell types (Davie et al., 2018). Regression models have been previously used to compare bulk transcriptome profiles to aggregated cluster profiles, including non-negative least squares (Davie et al., 2018). Here we use lasso regression to predict single cell cluster transcriptome profiles using the bulk profiles as variables. Therefore, each single cell cluster is decomposed as a weighted sum of bulk profiles. The regularization of the lasso model, allows for a higher signal to noise ratio and clearer assignments. Clear matches between a cluster and a bulk profile are detected as high coefficients (>0.1) in the sum. The regression model was fit using gene signatures derived from Scanpy as features. Before fitting, all datasets were scaled to counts per million CPM using the common genes, independent of previously applied normalization techniques. The calculations were performed using the SciPy (Virtanen et al., 2020) and Scikit-learn packages in Python.

Human brain samples

Usage of post-mortem brain samples was ethically approved (EC reference NH019 2019-02-01). All tissue samples were obtained from the Parkinson’s UK Brain Bank at the Imperial College London. The Parkinson’s disease cases used in this study all carried a mutation in the LRRK2 gene (P1542S, M1646T or R1325Q). One case carried a mutation in PARKIN in addition and one case carried a mutation in EIF4G in addition. Control samples are from cases without pre- or post mortem Parkinson’s disease-associated pathology and diagnosis.

Human nuclei isolation

Nuclei of human samples were prepared as described in (Hodge et al., 2019; Slyper et al., 2020).

Calling SNPs

In the human data, 10 different genotypes (5 PD cases, 5 controls) were mixed across 3 regions. To perform demultiplexing, the SNPs are required per sample. Since no genotype information was available, we performed bulk RNA-seq on all samples and called SNPs on the transcriptome. For this we used BCFFtools (v1.11, (Li et al., 2022)). We performed a pileup using the common loci of the 1000 Genomes project (1000 Genomes Project Consortium et al., 2015; Gramates et al., 2017) on the Hg38 iGenomes genome (-d 8000, -Ou), calling SNPs with the bcftools call function (-mv, -Ob). Using the reheader function we then updated the bcf file to make it compatible with demuxlet. We detected 146,848 SNPs, ranging between 25k and 51k per genotype (see table below). Next, the bcf file was sorted in the same order as the bam files and bcftools +fill-tags was used to add frequencies.

Demuxlet

Bam files were filtered using SAMtools (1.11, (Li et al., 2022)) to only contain reads overlapping with the vcf file. Then popscle dsc-pileup was used to create a pileup file for input in popscle demuxlet (Kang et al., 2018). All cells labeled as doublets were removed (∼4-7 %), and cells assigned to singlets were kept (25-43 %). Additionally, a large proportion of cells (∼49-69 %) were labelled as ambiguous. These cells were assigned to their most likely singlet identity and kept in the remaining analysis. In a next step, we only kept possible genotype-experiment combinations, ending up with 97,489 cells.

Data filtering and clustering

CellRanger generated a gene expression matrix of 151,368 cells. After removing doublets and wrongly assigned cells, 97,489 cells remained. The data was loaded into scanpy (Wolf et al., 2018) and further filtered using min_genes=500, max_genes=7000, min_counts=1000 and max_counts=30000, and finally with percentage of mitochondrial counts below 5. We lowered the allowed percentage of mitochondrial counts since the human data was obtained from snRNA-seq. The data were log-normalized with a scale factor of 10⁴ and scaled (maximum values = 10), regressing out the number of UMIs and the percentage of mitochondrial genes as latent variables using a linear regression model. Highly variable genes were selected using sc.pp.highly_variable_genes with default parameters.

We used sc.tl.pca to calculate 200 principal components on the highly variable genes, from which the first 75 components were selected for follow-up analysis based on an elbow plot. We again used Harmony (Korsunsky et al., 2019) to correct the principal components for batch effects, using genotype as the batch variable. The batch effect-corrected components were then used as input to calculate neighbors using sc.pp.neighbors with n_neighbors=30, followed by Leiden clustering for resolutions between 1 and 10. Subsequently, tSNE and UMAP visualizations were calculated.

Using marker genes for main neuronal and glial populations we chose to combine Leiden resolution 1 for glia and Leiden resolution 3 for neurons.

Gene set enrichment was performed using the AUCell function from the pyscenic package (0.10.4) in Python. The genes detected in affected cell types of Drosophila models were converted to human homologues using the FlyBase database, selecting homologues with a DIOPT score>=5. Significant enrichment was determined using a two-tailed t-test.

Differential expression

For every cluster, cells from each mutant genotype were compared against cells from the wild-type control. In the fly, OPN subtypes were merged before calculating differentially expressed genes (DEG). To control for biases in methods, we used Wilcoxon and DESEQ2 (Love et al., 2014). The Wilcoxon test was performed with the scanpy wrapper sc.tl.rank_genes_groups. For DESEQ2 we used rpy2 to call the R package DESEQ2 in Python. DESEQ2 was originally designed to be used on bulk RNA-seq data. Therefore, we grouped cells into pseudobulks based on their experimental run by summing their gene expression profiles. Only pseudobulks consisting of at least 10 cells were kept. We did not select any additional variables to regress.

To detect any method biases, we calculated the spearman correlation of the signed p-value (log-foldchange multiplied with log(p-val)), showing that both methods differ in the exact number of DEG, but agree quite well in their ranking of up- and downregulated genes. We found larger differences in the number of DEG called in the human data (0-10 in DESEQ2, 200-1800 in Wilcoxon) but the rankings remained highly correlated (r∼0.8). The Wilcoxon test is more sensitive to genes expressed in all cells in the cluster, while DESEQ2 can be oversensitive to genes expressed in only a minor fraction of cells.

Additional runs were made in which the number of cells per cluster was down sampled to correct against cluster-size based biases.

Modeling of DEG

Since the number of DEG has been shown to be highly correlated with the number of cells in single-cell RNA-seq, we used a modeling approach to find outlier clusters with more than expected DEG. Since the data is over dispersed count data, we fit a negative binomial model to the data, with number of differential genes (nDEG) as dependent variable and number of cells in the mutant and number of cells in the control as independent variables:

With y_i being the number of differential genes detected in cluster i, X_1ithe log of the number of cells of cluster i in the mutant and X_2i the number of cells of cluster i in the wildtype.

Models were fit using the statsmodels package in Python and goodness of fit was measured using the Log-Likelihood (Supplemental File 4). The number of UMIs per cluster was not included as independent variable, since it has been implicated with cell size and thus could remove biological effects. Residuals were calculated as the difference between the measured number of differential genes and the predicted number, wherein a positive residual means that the cluster has more differential genes than expected.

Olfactory behavior assay

Olfactory behavior was quantified in mixed populations of ∼50 5±1 day old flies. Flies were flipped into the experimental setup without anesthesia and allowed to adjust for some minutes to the experimental environment (24 °C, 50 % humidity). They were placed into a custom-built T-Maze (modified from (Tully & Quinn, 1985)), where they were allowed to freely walk between two tubes for one minute. Throughout the test, an aquarium pump delivers two separated, constant airstreams (1 l/min) which are directed through an odor vial before entering the tube. Odor vials contained: 2 ml motor oil, 2 g fresh banana, 2 ml undiluted beer (Jupiler), or 2 ml wine (LesTruffiers, 2015) diluted 1:1000. Afterwards the olfactory preference index was scored: (#flies in tube with odor A - #flies in tube with odor B)/#total flies. Tests were pseudo-randomized by switching the position of odor tubes every trial. Control flies were tested each experimental session.

In vivo Ca²⁺-imaging and synaptic area

5±1 day old (unless noted otherwise) male flies were placed into a custom-build chamber and the brains were imaged at 5 Hz through a window in the head capsule using a wide field fluorescence microscope (Nikon) equipped with a Hamamatsu camera and a 20×/NA=0.9 water-immersion objective. 100 mM nicotine was injected into the Ringer’s solution covering the brain (final concentration ∼10 mM). Images were analyzed using Fiji (Schindelin et al., 2012). Fluorescence intensity was quantified in a region of interest (ROI, d=10 µm) that was placed in the center of the synaptic area (calyx for OPNs, mushroom body lobes for dopaminergic PAM neurons) and normalized by subtracting background fluorescence. Basal fluorescence F₀ in the synaptic region was determined by averaging five frames before stimulus onset and the change of fluorescence over time was calculated by F-F₀/F₀. Flies of related genotypes and controls were always measured within the same sessions, in alternating order.

The area of the calyx was determined by measuring the (2D) surface of the entire visible fluorescent area.

Nicotine treatment in flies

Nicotine (N3876, Sigma) was mixed with fly food to reach a final concentration of 0.2µl/ml. A fly weighs ∼0.1 mg (Wu et al., 2016) and consumes ∼1 µl food a day (Deshpande et al., 2014). This corresponds to 0.002 µg nicotine uptake per mg body weight daily and corresponds to the amount a smoker absorbs with 10 cigarettes a day. Fly vials containing nicotine and their controls were kept in constant darkness, as nicotine is light-sensitive.

Generation of LRRK2^G2019S hIPSCs

hiPSC from the KOLF2-1J (Pantazis et al., 2022) were used as control. LRRK2 p.G2019S was engineered in the KOLF2-1J iPSC line, by means of CRISPR/Cas9 gene editing following (Skarnes et al., 2019). Briefly, KOLF2-1J iPSCs were nucleofected using the Lonza 4D Nucleofector with Cas9 HiFi protein (IDT), a synthetic guide RNA (Synthego) and a single strand oligodeoxynucleotide (ssODN; IDT) as a donor template. Four days post-transfection, cells were dissociated to single cells and seeded at low density in 10 cm plates. Single colonies were isolated in 96-well plate and sampled for molecular analysis by PCR followed by Sanger sequencing. More than 90% of the colonies harbored the LRRK2 p.G2019S mutation in homozygosity. Selected clones were karyotyped by means of comparative genomic hybridization (CGH) and stained for pluripotency markers.

To build the transfer vector pLenti-hSynI-mRuby2-T2A-GCaMP6f (Addgene plasmid # 197595), the hSynI-mRuby2-T2A-GCaMP6f-Wpre fragment was PCR amplified from pAAV-hSyn1-mRuby2-GSG-P2A-GCaMP6f-WPRE-pA, which was a gift from Tobias Bonhoeffer & Mark Huebener & Tobias Rose (Addgene plasmid #50943 (Rose et al., 2016)) and was introduced into a regular 3^rd generation lentiviral backbone.

iPSC differentiation to ventral midbrain dopaminergic neurons

iPSCs are first differentiated to ventral midbrain neural progenitor cells before they are terminally matured to ventral midbrain dopaminergic neurons using established protocols (Kriks et al., 2011; Nolbrant et al., 2017). In summary, on day -1, 400.000 hiPSCs/cm2 were seeded on Matrigel-coated wells of 6-well plates in StemFlex medium supplemented with 10 µM Rho kinase inhibitor (RI). On day 0, medium was switched to a 1:1 mix of DMEM/F-12 and Neurobasal supplemented with 0.5x N2, 0.5x B27 without vitamin A, GlutaMAX, Penstrep, non-essential amino acids, β-mercaptoethanol (10 μM) (all Life Technologies), insulin (5 μg/ml, Sigma), LDN193189 (500 nM, Sigma), SB431542 (10 μM, Tocris), SHH-C24II (200 ng/ml, Miltenyi Biotec) and Purmorphamine (1 μM, Sigma). CHIR99021 (1.5 μM, Stemcell Technologies) was added to the medium from day 3 to day 13. SB431542, SHH-C24II and Purmorphamine were withdrawn from the medium at day 9. FGF8b (100 ng/mL, Miltenyi Biotec) was added to the medium from day 9 until day 16. Cells were split 1:1 at day 11 and medium was switched to Neurobasal supplemented with 0.5x N2, 0.5x B27 without vitamin A, Penstrep, GlutaMAX. At day 16, cells were split 1:2 and the medium was switched to terminal differentiation medium consisting of Neurobasal-A supplemented with 1X B27 supplement without vitamin A, GlutaMAX, Penstrep containing BDNF (brain-derived neurotrophic factor, 10 ng/ml; Miltenyi Biotec), ascorbic acid (0.2 mM, Sigma), GDNF (glial cell line-derived neurotrophic factor, 10 ng/ml; Miltenyi Biotec), dibutyryl cAMP (0.2 mM; Sigma), SR11237 (100 nM, Tocris), and DAPT (10 μM; Tocris). On day 20, ventral midbrain neural progenitors were cryopreserved and quality controlled. To obtain mature neurons, neural progenitor cells were terminally differentiated on coverslips previously coated with poly-D-lysine (50 µg/mL, Life Technologies) and mouse laminin (1 µg/mL, Sigma) in terminal differentiation medium for an additional 40 days.

Characterization of iPSCs, progenitor cells and mature neurons

For the characterization of iPSCs, ventral midbrain floor plate progenitor cells and mature dopaminergic neurons were fixed for 15 min in 4 % formaldehyde at the pluripotency stage, day 20 of differentiation or day 60 of differentiation, respectively. Cells were blocked for 1 h at room temperature with 3 % normal goat serum + 0.3 % Triton X-100 (Sigma) in DPBS supplemented with Ca²⁺ and Mg²⁺ (Life Technologies). Primary antibodies were incubated overnight at 4 °C in blocking solution. Secondary antibodies were incubated for 1 h at room temperature in blocking solution. Coverslips were mounted in Mowiol (Sigma) and imaged on an upright Nikon A1R confocal microscope equipped with a DIC N2 20X lens (NA 0.75). Z-stacks were acquired with pinhole of 1 Airy unit, a Galvano scanner with line averaging of 2, image size of 1024 x 1024 pixels and step intervals of 2 µm.

The following antibodies were used: mouse IgG1 anti-SOX2 [1:200 (Santa Cruz)], rabbit anti-OCT4 [1:50 (Abcam)], mouse IgG1 anti-NANOG [1:50 (Santa Cruz)], mouse IgM anti-TRA-1-81 [1:100 (Abcam)], rabbit anti-LMX1A/B [1:1000 (Millipore)], mouse IgG2a anti-FOXA2 [1:250 (Santa Cruz)], mouse IgG2a anti-OTX2 [1:100 (Santa Cruz)], rabbit anti-TH [1:500 (Sigma)], mouse IgG1 anti-MAP2 [1:1000 (Sigma)], Alexa Fluor-488/ Alexa Fluor-555/ Alexa Fluor-647 conjugated secondary antibodies [1:500 (Invitrogen)].

In vitro ca²⁺ imaging and drug treatment

Day 40 dopaminergic neurons were treated for the following 19 days with 1 µM nicotine or with 1 µM nicotine plus 1 µM of the nicotinic receptor antagonist Mecamylamine hydrochloride (Sigma). At day 50, dopaminergic neurons were transduced with a lentiviral vector expressing mRuby2-T2A-GCaMP6f (Addgene plasmid # 197595). Transduction efficacy was assessed by directly visualizing mRuby2 red fluorescence. Drug treatment was interrupted 1 day before analysis. For calcium imaging, cell culture medium was switched to BrainPhys Imaging medium (Stem Cell Technologies) supplemented with B27 plus (Life Technologies), PenStrep (Life Technologies), GDNF (10 ng/mL, Miltenyi Biotec), BDNF (10 ng/mL, Miltenyi Biotec), dibutyryl cAMP (0.2 mM, Sigma) and SR11237 (100 nM, Tocris).

Cells were imaged at 37 °C at a wide field fluorescence microscope (Nikon) equipped with a Hamamatsu camera and a 20×/NA=0.9 water-immersion objective. For each plate, at least 2 areas were imaged for six minutes at a framerate of 5 Hz. Images were analyzed using Fiji (Schindelin et al., 2012). All visible cell bodies in the field of view defined a region of interest (ROI), background fluorescence was subtracted to correct for bleaching. Basal fluorescence F₀ in each ROI was determined by the minimal average of five frames and the change of fluorescence over time was calculated by F-F₀/F₀. Values above 2xStdDev of the baseline frames were considered as activity. Activity is expressed as active frames/total measured frames. The maximal amplitude was counted as amplitude for each neuron.

Quantification and statistical analysis

GraphPad Prism was used for visualization and statistical analysis of results from all but the sequencing experiments. Male and female flies were quantified separately, and results pooled later because they did not show significant differences. Data sets were checked for normal distribution and then analyzed with ANOVA and subsequent multiple pairwise comparisons using Dunnett’s test when comparing to a general control, or Bonferroni-correction of p-values when comparing different control-test pairs.

Key resources table

Acknowledgements

We thank Ann Geens, Eren Can Eksi, Vinoy Vijayan, Liesbeth Deaulmerie, Willem van den Bergh, Marianna Decet and Sara Aibar, Bart De Strooper, Joris de Wit, Roman Praschberger, the members of the Verstreken and Aerts labs and the VIB BioImaging core for help and comments. We thank Atefeh Pooryasin and the Bloomington Drosophila Stock Center (NIH P40OD018537) for fly lines. Human tissue was kindly provided by the Parkinson’s UK Brain Bank at UCL, funded by Parkinson’s UK, a charity registered in England and Wales (258197) and in Scotland (SC037554).

Additional information

Author contributions

Conceptualization: UP, SV, PV

Methodology: UP, JJ, NS, SK, SM, GH, SP, JS, SV

Investigation: UP, JJ, AB Visualization: UP, JJ, PV

Funding acquisition: PV, SA, UP, SV, AB

Project administration: PV

Supervision: SA, PV

Writing – original draft: UP, PV Writing – review & editing: all

Ethical approval

Usage of post-mortem brain samples was ethically approved (EC reference NH019 2019-02-01).

Funding

Research support was provided by ERC CoGs (PV and SA) and ERC AdG (PV), the Research Foundation Flanders (FWO), a Methusalem grant of the Flemish government, IMI2, Opening the Future and VIB. UP was supported by a fellowship from DFG; AB was supported by an EMBO long term fellowship and SV was supported by a fellowship from FWO. PV is an alumnus of the FENS Kavli Network of Excellence.

Abbreviations

• DAN: Dopaminergic neurons

• DEG: Differential expressed genes

• ERG: Electroretinograms

• KC: Kenyon cells

• Nac: Nucleus accumbens

• NBM: Nucleus basalis Meynert

• OPN: Olfactory projection neurons

• PD: Parkinson’s disease

• Put: Putamen

• SING: Startle induced negative geotaxis

• iPSC: Induced pluripotent stem cells

Additional files

Supplemental Text

Supplemental Figure 1

Supplemental Figure 2

Supplemental Figure 3

Supplemental Figure 4

Supplemental File 1

Supplemental File 2

Supplemental File 3

Supplemental File 4

Supplemental File 5

Supplemental File 6

Supplemental File 7

Supplemental File 8