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. In this study, we generated and utilized five well-controlled Drosophila Parkinsonism models to investigate cellular and molecular dysfunction preceding motor defects that occur later in the life of these animals.

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 models that was extensively backcrossed into one isogenic background (>10 generations). We selected homologues of 5 PD genes that represent different affected pathways and clinical manifestations in this disease (Kalia & Lang, 2015) 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). 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 (Fig. 1a, 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 PD-gene knock-ins are similar to controls, 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 reduced lifespan under starvation, disturbances in mitochondrial membrane potential and electrophysiological defects in electroretinogram (ERG) recordings (Hindle et al., 2013; Hoffmann & Harshman, 1999; Morais et al., 2009; Mortiboys et al., 2010; Ng et al., 2012; Vanhauwaert et al., 2017) (Fig. 1b, Supplemental Fig. 1b-d). Importantly, we also tested the flies for their performance in a startle-induced negative geotaxis assay (SING). This is a motor behavior that critically depends on specific 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 alpha-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 age (25 d, old), the mutant knock-ins display impaired performance in this DAN-dependent locomotion assay, while the performance of controls is not affected (Fig. 1b). This suggests that the mutants we used in this study suffer from a progressive locomotion defect that has been linked to DAN dysfunction.

Fig. 1.

Cholinergic olfactory projection neurons are impaired across 5 PD 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). 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 for 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 3). Using the residuals of this model we find eight common predominantly affected cell types across the 5 mutants (Supplemental File 4). 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). Hence, while our mutants are modelling different genetic origins of disease, we find some cell types are commonly affected. 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). Furthermore, hyposmia is one of the earliest defects associated with familial and sporadic Parkinsonism (Chase & Markopoulou, 2020; Doty, 2012) with patients having difficulties to classify odors in standardized UPSIT tests (Double et al., 2003; Mueller et al., 2005; Witt 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 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 neuronal subtypes, including cholinergic neurons (Supplemental Fig. 3c, d) and identified the DEG in these cells. 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’). We then asked how they compare to the DEG found in affected cholinergic projection neurons of our five fly models (Fig. 3a”). 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 5). 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 DEG that have an orthologous gene in the fly are also found among the most deregulated genes across PD fly models.

Fig. 3.

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”). 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 (Gorenberg & Chandra, 2017; Morais et al., 2009; Snead & Eliezer, 2019; Uytterhoeven et al., 2011). In summary, cholinergic neurons of young familial fly PD models and human PD patients both show preferential transcriptional deregulation of genes that regulate synaptic function (Fig. 3c).

Young PD models suffer from synaptic impairment of cholinergic OPN

The transcriptional deregulation of genes encoding synaptic proteins predicts 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 (10 mM Nicotine) results in a robust synaptic Ca²⁺ signal at OPN synapses in controls (Fig. 4c, d (black points), quantified 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 all 5 mutants quantified in Fig. 4g (“OPN>wt CDS -”)). Next, we tested if this defect is cell-autonomous and created PD knock-in mutants expressing wild-type CDS selectively in their OPN (Fig. 4g (“OPN>wt CDS +”)). These animals show synaptic Ca²⁺ signals similar to controls, indicating the PD genes are required in OPN of these young animals to maintain robust synaptic function. Note that while hLRRK2 knock-in flies are very similar to wild-type controls (Fig. 4g (left bar)), the hLRRK2^G2029S knock-in flies are not rescued by OPN-specific expression of hLRRK2 (Fig. 4g (right bar)). 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 pathogenic mutant CDS knock-in flies 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, except in LRRK2^G2019S (Fig. 5c, d), again agreeing with it being a dominant mutant. Note that lrrk^KO flies also show an olfactory performance defect and that expressing hLRRK2 cDNA specifically in the OPN of lrrk^KO rescues this defect (Fig. 5c, d). 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^G193R mutants), does not improve olfactory performance (Fig. 5c, d (5^th row)). Hence, a cell autonomous defect in cholinergic OPN causes olfactory defects across PD mutants.

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 our own unpublished data, showing 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). These findings imply that factors affecting the function of cholinergic neurons might, by absence of insufficient innervation, 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). Fly olfactory projection neurons, like human basal forebrain cholinergic neurons and striatal cholinergic interneurons, are also characterized by large neuronal arbors, numerous synaptic contacts and high activity (C. J. Wilson et al., 1990; R. I. Wilson, 2013; Wu et al., 2014). Such neurons are thought to be particularly vulnerable to synaptic problems (Ashrafi et al., 2014; Vijayan & Verstreken, 2017) and our work now indicates that the defects induced by different familial PD mutations converge onto synapses. This can be through different upstream mechanisms (i.e. the causal genes affect different cellular processes), but what they have in common is that this ultimately induces synaptic failure. This is particularly relevant in such elaborate cholinergic projection neurons where synapses are often at a long distance from their cell body and synapses need to operate relatively independently. As such there can be different mechanisms at play and 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.

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

Declarations

Ethical approval

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

Data and materials availability

Plasmids and fly lines generated in this study are available on request. Data sets generated in this study are stored in public repository GEO https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE235332 with restricted access until publication. Access will be given to reviewers via tokens.

Competing interests

Authors declare that they have no competing interests.

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.

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 Prashberger, 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).

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

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=GSE235332.

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, oligos, gBlocks and fly strains used will be provided. To clone inserts into plasmids NEBuilder® HiFi DNA Assembly (NEB) was used.

Unique gRNAs (Supplemental File 1) 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 (CantonS 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.

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 of all Parkinson’s disease models, Pink1 we used the expression of human PINK1 with UAS-hPINK1 and nsyb-Gal4. This is because, unlike the two pathogenic variants, we were unsuccessful in targeting the wild-type human PINK1 CDS to the fly pink1 locus.

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^R258Q were 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 the 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 and transcribed using Maxwell RSC instrument and Maxwell RSC RNA kit (Promega) and quantified by subsequent qPCR using a LightCycler 480, 480 SYBR Green master mix (Roche).

Survival assay

Flies were collected into populations of 10 males and 10 females the day after eclosion. Five days after eclosion, these populations were transferred into vials containing 1 % agarose for wet starvation. Dead flies were scored and removed every 8 h, the remaining flies were flipped on fresh agarose every 24 h. The assay was repeated with at least two different batches of a given genotype. Median survival and survival curves were determined using GraphPad Prism.

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 MgCl2 • 6H2O, 10 NaHCO3, 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.

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.

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, (Jenett et al., 2012)). 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 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.

Seurat

Dimensionality reduction and clustering were performed using the Seurat package in R (Satija et al., 2015). The processed and filtered matrices of all samples were combined together with an annotated dataset of 57k cells of the adult brain containing two wild-type genotypes (Davie et al., 2018), leading to a total number of ∼170k cells. In the first steps of the Seurat pipeline, each genotype was processed separately. The data were log-normalized with a scale factor of 104 and scaled, regressing out the number of UMIs and the percentage of mitochondrial genes as latent variables using a negative binomial model. For further downstream analysis, the 2000 most variable genes were selected using FindVariableGenes and method=’vst’ with default parameters. Then 3000 integration features were determined across all genotypes with SelectIntegrationFeatures and used in the rest of the analysis. At this point we looped over a number of components to test the dimensionality of the dataset, ranging between 30 and 100 components. These components were used as input for the integration functions, PCA, tSNE and clustering. Integration was performed by linking cells between datasets (so-called anchors, function FindIntegrationAnchors), followed by integrating the datasets one by one (IntegrateData). The scaled and integrated matrix was used as input for PCA to reduce the dimensions. The principal components were then used to calculate the non-linear dimensionality reduction tSNE (RunTSNE), and to perform clustering with Louvain Shared-Nearest Neighbor (SNN) graph algorithm (FindNeighbors, FindClusters). Different resolutions ranging between 1 and 8 were tested. In general, results are only slightly affected by the number of components used (Seurat tutorial), so we decided to compare the tSNE and clusters qualitatively, settling on 70 as the final number.

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; Jenett et al., 2012) 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 Seurat 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). 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 to correct the principal components for batch effects. 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_li the log of the number of cells of cluster i in the mutant and X_2i the number of cells of cluster i in the wild type.

Models were fit using the statsmodels package in Python and goodness of fit was measured using the Log-Likelihood (Supplemental File 3). 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.

Quantification and statistical analysis

GraphPad Prism was used for visualization and statistical analysis of results from all but the sequencing experiments. Males and females 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