Behavioral screening defines the molecular Parkinsonism-related subgroups in Drosophila

Domain: biology
Published: 2026-03-10
Editorial outcome: major_revision

Abstract
Parkinson’s disease (PD) and related familial Parkinsonism are defined by motor dysfunction, but the specific upstream molecular causes of these clinical symptoms can vary widely. We hypothesize that these causes converge onto a limited number of core cellular pathways. To investigate this, we created a collection of 24 genetically well-controlled Drosophila models of familial forms of PD and related mono-genic forms of Parkinsonism. Using unbiased behavioral screening and machine learning we identify clusters of mutants that converge on (1) mitochondrial function; (2) retromer/vesicle trafficking and proteostasis/autophagy. Genes within each cluster have a similar genetic interaction profile and compounds that target specific molecular pathways ameliorate dopaminergic neuron dysfunction in a cluster-specific manner. Together, our data indicate that familial PD and related forms of Parkinsonism may fall into two broad functional groups, and may inform further work toward targeted biomarker discovery and therapeutic development.

Introduction
Parkinson’s disease and related genetic causes of Parkinsonism (for simplicity referred to as “Parkinsonism” hereafter) are characterized by motor symptoms such as bradykinesia, rigidity and tremor, caused by the progressive dysfunction and loss of dopaminergic neurons in the substantia nigra. However, patients often report non-motor issues decades before the onset of these motor symptoms, including hyposmia, and various sleep disturbances. In some cases, patients may also experience dementia and epileptic seizures1–5 (Supplementary Table 1). While therapies to restore dopamine tonus can relieve motor symptoms6, they do not address non-motor problems and they do not stop disease progression5,7,8. These treatments can also lead to side effects including hallucinations, psychosis and agitation9. This underscores the need for developing disease-modifying strategies that are grounded in the biological mechanisms of Parkinsonism.

While there are common motor and non-motor symptoms in Parkinsonism, the clinical presentation may vary10–13. There is work using these variations in an attempt to define disease subtypes14–16, but these strategies have met limitations. Factors contributing to this challenge include the inherent genetic heterogeneity among patients and the limited availability of prodromal and longitudinal quantitative data. Furthermore, these clinical classifications pose experimental difficulties for therapeutic development as they are not associated with specific genetic mutations and thus molecular dysfunction17.

An alternative approach to classify Parkinsonism involves exploring its genetic complexity, providing a pathway-based perspective for stratification. Over 20 genes have been linked to familial forms of Parkinsonism, and recent genome-wide association studies (GWAS) have identified a growing number of loci that modify the risk of Parkinsonism18–21. These genetic findings offer valuable insight into the molecular mechanisms that potentially drive the disease22–24. However, the diversity of the potentially affected pathways poses the question of how, or whether, these Parkinsonism-related genes converge on common molecular networks to result in the same clinical disease. These molecular dysfunctions could manifest before dopaminergic-dependent motor symptoms given that patients experience complex behavioral and functional disturbances in the prodromal phase25. Thus, stratifying Parkinsonism based on molecular dysfunction could have significant implications for clinical trials and personalized treatment approaches already at an early stage.

In this study, we established a collection of 24 genetic Drosophila models for Parkinsonism and leveraged these to identify distinct disease subgroups. Using high-throughput behavioral monitoring, machine learning and pattern recognition algorithms, we were able to classify individual mutants in an unbiased manner. Our analysis revealed two major groups, whose existence was confirmed through extensive genetic interactions observed within each group. These groups share common molecular pathways related to mitochondrial function, vesicle trafficking/retromer and autophagy/proteostasis, and we demonstrate that targeted small molecule treatments can be tailored to these specific subgroups. This behavioral screening approach thus enabled us to categorize Parkinsonism into two major non-overlapping molecular pathways, providing a framework for understanding the disease and developing targeted treatments.

Results
Parkinsonism fly collection recapitulates dopamine responsive, progressive motor impairments.

Using CRISPR/Cas9-mediated gene editing we created a comprehensive collection of Drosophila loss-of-function Parkinsonism models. We identified the 24 fly homologues of genes known to be causative or conveying risk for Parkinsonism in humans, based on extensive literature research and clinical insights. We replaced the first common exon of possible transcripts by a attP-flanked w+ reporter cassette in each of the genes (Supplementary Table 1, 2, Supplementary Fig.1a-b, using the donor plasmid pWhiteSTAR (“WLS”), methods). These knockouts (KO) were backcrossed to isogenic control flies for >10 generations, allowing for detailed comparisons between these models independent of genetic variation - for park we were unable to generate such a knock-out and thus we resorted to the established park1/Δ21 mutant 26,27. Most of the Parkinsonism genes we included are thought to cause Parkinsonism or increase risk of the disease by a loss-of-function mechanism with some exceptions; nonetheless, we included the analysis of KO flies only as a way of comparing the genetic factors to one another without confounding our analysis at this stage with the effects of specific pathogenic mutations. These KO flies are loss of function mutants as confirmed by RT-PCR (Supplementary Fig.1c) and are collectively referred to as ‘Parkinsonism mutants’.

To assess whether these Parkinsonism mutants display locomotion phenotypes, we used a well-established climbing assay (SING - startle-induced negative geotaxis). SING defects in Parkinsonism fly models have previously been shown to involve dopaminergic impairments28–30. The performance of young Parkinsonism KO flies is similar to controls, except for gbaKO flies that already exhibit a reduced SING score at young age (Fig.1a). When aged, we detect a significant decline in motor performance in 80% of the Parkinsonism mutants (Fig.1b), consistent with age-dependent dopaminergic impairments. To test this, we labeled dopaminergic neurons with anti-TH, which allowed us to assess synaptic connectivity onto their downstream target, the mushroom body (MB). In line with the SING results, we find that young Parkinsonism KO flies have normal innervation onto MB compared to controls (measured by the area of positive TH labeling, normalized to MB area as a measure of individual mutant brain size, Fig.1d, Supplementary Fig.2a) except for gbaKO flies. This TH-labeled innervation area progressively declines (Supplementary Fig.2b) and ultimately 90% of the Parkinsonism mutants with motor defects also show a significant dopaminergic neuron innervation defect (Fig.1c,d and Supplementary Fig.2c-f). We find that impairments in dopaminergic innervation onto MB significantly correlate with the defects in SING motor performance across our Parkinsonism mutants; R2: 0.53 (Fig.1e).

We next tested if Parkinsonism mutant flies with a progressive motor defect (p≤ 0.001) are L-Dopa responsive6,31. The motoric defects of most Parkinsonism mutants are significantly alleviated when flies are fed L-Dopa in the 10 days prior to the SING assay. We do not observe such rescue when flies are fed the stereoisomer D-Dopa (Fig.1f) or H2O (solvent control, for selected mutants, Supplementary Fig.2g)27,32. Thus, the majority of our Parkinsonism mutants suffer from a progressive L-Dopa-responsive motor defect.

In addition to typical motor symptoms, Parkinsonism patients sometimes experience additional progressive neurological defects, including epileptic seizures4 (Supplementary Table 1). Although seizure-like activity in flies does not reproduce the full pathophysiology of human epileptic seizures, it serves as an accessible model of hyperexcitability.33–36 We assessed seizure-like behavior using an established protocol based on mechanical hyper-stimulation of sensory inputs and assessing seizure-like activity (methods;36). While controls are fine, eleven aged Parkinsonism mutants show such behavior (Supplementary Fig.3a-c, Supplementary videos 1-6). While SING defects in these flies are rescued by L-Dopa, seizure-like activity is not (Supplementary Fig.3e). We then compared the familial Parkinsonism cases reported to suffer from epileptic seizures to the fly mutants exhibiting seizure-like behavior: 7 of the 24 familial forms of Parkinsonism were described to suffer from seizures37–41, and 6 of these 7 homologous fly mutant models also display seizure-like activity (Supplementary Fig.3d). There are also 5 other fly models with seizure-like behavior where clinical assessment of patients has not (yet) reported epileptic behavior (Supplementary Fig.3d).

Taken together, our collection of Parkinsonism mutants recapitulates hallmarks of the disease, including typical age-dependent and dopamine-dependent neurological defects and defective dopaminergic neuron innervation. Several of the models also show epileptic-like behavior.

High throughput behavior monitoring reveals Parkinsonism subgroups with distinct sleep and activity patterns.

We used the Parkinsonism mutants to stratify the disease based on comparing complex behavioral phenotypes across our genetic models. We placed individual flies in in-house-built ethoscopes that enable automated video tracking of movements over multiple diurnal cycles (Fig.2a)42–44. We measured day-time activity, but also multiple other parameters such as features relevant to sleep. Sleep is an evolutionary conserved behavior45,46 and various aspects of the sleep cycle are affected in the prodrome of Parkinsonism such as insomnia, excessive daytime sleepiness, and other alterations47,48. We first calibrated our system using the pigment-dispersing factor (pdf) null mutant, which is a well-studied circadian deficient fly line. As expected, we observe defects in evening and morning anticipation (Fig.2bi and 2biii), reduced latency and increased sleep at night compared to control (Fig.2bii) thus validating our experimental system49.

We then screened the entire Parkinsonism fly collection and tested ≥69 individual flies per genotype at a young age (2-5 days post eclosion) before the onset of dopaminergic motoric defects (Fig.1a), and analyzed their activity profiles comprising basic motor performance, circadian rhythms (morning anticipation) and the sleep architecture (sleep fraction in light and dark periods; number and duration of sleep bouts; onset of longest bout) (Fig.2a and Supplementary Fig.4a-d), aspects of which are also frequently found affected in patients2,50. We find specific behavioral defects and none are common across all the mutants, indicating this type of analysis is capable of distinguishing specific features between Parkinsonism models (Supplementary Fig.4). Furthermore, we also show that these defects are not the result of basic motor performance impairment (e.g. velocity if awake) as this feature does not generally correlate with the behavioral parameters we quantified (Supplementary Fig.5).

To analyze these results also in an unbiased and holistic manner, we employed unsupervised pattern recognition analysis, enabling us to decompose the ethoscope data of individual flies. Using non-negative matrix factorization (NMF)51,52 we identified distinct behavioral components (Fig.2c, d). The elbow method and Calinski-Harabasz Index showed that 5 components explained most of the biological meaningful variance (Supplementary Fig.4e,f, methods)51,53: ‘Component 0’ is defined by a sleep bout in the morning; ‘Component 1’ delineates increased sleep as the night progresses, especially in the early morning; ‘Component 2’ constitutes sleep late during the day; ‘Component 3’ is defined by a normal sleep pattern and ‘Component 4’ constitutes a strong sleep bout right at the start of the night. Each individual fly is represented by multiple behavioral components. For example, control flies are represented mostly by component 3 (Supplementary Fig.4g-k). In contrast omiKO-WS is mostly represented by component 1, and less by components 0 and 3, while Pink1KO-WS and park1/Δ21 flies are characterized chiefly by components 0, 2 and 4, but also to a lesser extent by the other components (Supplementary Fig.4g-k). The mean principal component value of these features that represent each fly per genotype provides an unbiased quantification of the behavioral phenotype for each model. Also in these analyses, the components do not generally correlate with basic motor performance defects and with the behavioral parameters (Supplementary Fig.5, 6, 7), meaning they can be treated as independent variables.

We then combined the unsupervised NMF-defined phenotypes with the explicit behavioral features and used hierarchical clustering to identify subgroups characterized by similarities in behavioral features (Fig.2e-g). This delineates 2 major groups, each comprising 3 and 2 subgroups, respectively. Providing credence to our approach, Pink1KO-WS/y and park1/Δ21, whose protein products are known to interact to regulate mitophagy54,55, cluster closely in group B1 (orange label). Furthermore, they also cluster together with another mitochondria-related gene Chchd256, suggesting behavioral clustering reveals common underlying biology. Furthermore, the adjacent group (B2, light orange) also includes genes involved in mitochondrial function such as nutcracker (FBXO7) and Coq257,58. Among other elements, their behavior is characterized by increased sleep in the morning (component 0, Fig.2e). The sleep traces of park1/Δ21 flies and Pink1KO-WS reveal a shortened latency to fall asleep (Fig.2gi) and, as also previously reported, impaired anticipation of morning “lights-on” timing59 (Fig.2gii). On the other hand, gbaKO, iPLA2-VIAKO-WS, DJ-1a/bKO-WS, Vps35KO-WS/+ and Rab39KO-WS flies cluster in ‘group A1’ (dark blue label) and show a reduced sleep fraction (Fig.2fi). Similarly, the closest neighboring group (A2, medium blue label) also shows reduced sleep, especially during the day. Finaly, group A3 (light blue label) comprises Parkinsonism mutants such as VAC14KO-WS/+, VhaM8.9KO-WS/+, LrrkKO-WS, Rme-8KO-WS/+ and eIF4GKO-WS/+ showing e.g. reduced sleep late during the day (component 2) and increased numbers of sleep bouts per night.

We verified the existence of these groups using an independent methodology: principal component analysis (PCA). This reveals similar groups of Parkinsonism-gene KO animals (Supplementary Fig.8, groups outlined as in Fig.2e) as those found using hierarchical clustering, except for gbaKO that in the PCA associates with the genes in behavior group B (orange). This suggests the behavior of gba mutants shares characteristics of both groups A and B (Supplementary Fig.8a). PC1 and PC2 together explain 58% of the total variance (Supplementary Fig.8b) and are strongly influenced by individual behavior features (Supplementary Fig.8c). PC1 is influenced by component 3 and 4 of our NMF analysis, the number of sleep bouts and the velocity if awake , while PC2 is mostly defined by sleep fraction, bout length and component 2 (Supplementary Fig.8c). PC1 delineates the difference between the subgroups in groups A and B (light versus dark shades), while PC2 is most specific in delineating the difference between the mutants of groups A and B . Only the mutants in group A3 (light blue) are not discriminated well by our PCA. Hence, both types of unbiased analyses of the behavior parameters define 2 broad groups across 24 Parkinsonism mutants.

To assess the stability of the hierarchical clusters we removed one behavioral feature or NMF component at a time from the original analysis (Supplementary Fig.9a, b) and compared the newly generated groups (group 1-5) to the originally determined clusters (behavior A1 - B2) using Jaccard similarity (Supplementary Fig.9c-n), which computes the intersection-over-union51. This shows that in more than 50% of the cases when one feature is removed we find a Jaccard similarity higher than 0.35 (indicating substantial overlap60–62) (Supplementary Fig.9o). This indicates the robustness of the subgroups we defined.

To test this further we included the pdfKO mutant in our clustering. Circadian defects in park and Pink1 mutant flies have been associated with impaired biogenesis of dense core vesicles containing the neuropeptide PDF59. Hence, we recorded the behavioral features of pdfKO animals, and also defined the components by NMF. pdfKO flies show behavioral similarities with park1/Δ21, Pink1KO-WS/y and Chchd2KO-WS/y mutants (green, Supplementary Fig.10a-b) and, consistent with our prediction, pdfKO co-clusters with these mutants in behavior group B1 (orange, Supplementary Fig.10b).

Parkinsonism gene expression patterns do not correlate with subgroups.

The convergence of Parkinsonism mutants into subgroups based on similar behavioral phenotypes could find its origin in shared molecular mechanisms or expression patterns across cell types of the genes in one group. To test the latter, we utilized single-cell sequencing datasets from our young control flies30 and assessed the expression level of each of the fly Parkinsonism genes across all identified cell types in the brain (Supplementary Fig.11a and b) and the entire body (Supplementary Fig.11c). We then performed hierarchical clustering analysis on these gene expression profiles. For both the brain cell type dataset and the entire body cell dataset, this does not recapitulate the groups we identified based on behavior. For example, the expression profiles of Pink1, park and Chchd2 in the brain do not co-cluster and are expressed broadly (Supplementary Fig.11a and b). Similarly, gba1a/b, DJ-1a/b, Rab39, Vps35 and iPLA2-VIA that comprise group A1 do not exhibit similar or specific brain or body expression patterns or expression levels across cell populations (Supplementary Fig.11a-c). While we cannot exclude that genes could act within multicellular functional circuits comprising distinct cellular populations, we conclude that the clustering of genes based on behavior is not easily explained by shared expression profiles in the brain or body.

Genes in behavioral subgroups show similar genetic interaction profiles.

We assessed the validity and robustness of the behavior-defined subgroups by performing genetic interaction (GI) analyses across all Parkinsonism mutants in vivo at a young age, similar to the age used in the behavior screening. In this independent approach, our strategy was to combine all mutants one-on-one with each other and to use a sensitive and quantitative electrophysiological readout of neuronal function: electroretinogram recordings (ERGs) that were found to be disrupted in many neurological disease models in flies. ERGs are the extracellular voltage response of the retina to a transient light pulse63,64. We analyzed the amplitude of the ERG-depolarization of single heterozygous mutants and the combinations of double trans-heterozygous Parkinsonism mutants (Fig.3a,b). For Pink1KO-WS, Chchd2KO-WS and Rab39KO-WS we used hemizygous flies, as their genes are located on the X chromosome. We then used a previously established methodology to define the GI as the difference between the prediction from a non-interacting model using a “multiplicative function” of 2 single heterozygous mutant flies and the observed value of the double heterozygous gene pair65,66. This allowed us to define the GI of gene pairs as “aggravating”, “alleviating” or “not interacting”. For example, within the 276 combinations we tested (Supplementary Data 1), we find the phenotype of double heterozygous mutants of iPLA2-VIA and DJ-1a/b that are co-clustering in behavior group A (Fig.2e, blue), to be less severe (alleviating) as compared to either of the heterozygous mutants alone (Fig.3c). Thus suggests that they act in the same molecular pathway. Similarly, we find an alleviating GI of Rme-8 and VhaM8.9 (Supplementary Fig.12a). Combing mutants of Vps35 and omi or Pink1 and nutcracker, which co-cluster in behavior group A and B, respectively (Fig.2e), aggravates the phenotype seen in the single heterozygous mutants (Fig.3e and Supplementary Fig.12b, Supplementary Data 1), highlighting a strong GI between these genes. The latter is consistent with nutcracker’s function in mitophagy and its molecular interaction with the protein product of Pink157. In contrast, we also find a gene pairs from behavior group A, such as Gdh and Vps13, that show the same ERG amplitude as predicted from single heterozygous mutants suggesting no molecular interaction (Fig.3f).

When we perform hierarchical clustering on the GI profiles of the Parkinsonism mutants, we find 2 subgroups: GI 1 and 2 (Fig.3f). Interestingly, ‘GI 1’ comprises 6 out of 8 genes that also co-cluster in behavior group B (Fig.2e, orange), including Pink1, park, Chchd2, nutcracker, anne and auxilin. Furthermore, the majority of mutants that cluster in behavior group A (Fig.2e, blue) are found in ‘GI 2’, including gba, Synj, DJ-1a/b, iPLA2-VIA and Tango14 as well as Rme-8, VhaM8.9, Lrrk, VAC14 and eIF4G. We quantified the similarity of the behavior clustering approach and the GI approach using the Jaccard similarity index. This confirms that behavior group A significantly overlaps with GI 2 and behavior group B with GI 1 (Fig.3g). Thus, both the unbiased behavioral analyses and the GI profiles define similar Parkinsonism subgroups, which ultimately define shared molecular handles into different disease subtypes.

Cluster-specific rescue of dopaminergic neuron defects.

Next, we tested if the common age-dependent defects in dopaminergic neurons across the Parkinsonism mutants (Fig.1 and Supplementary Fig.2) have a cluster-specific pathomechanistic origin. We targeted pathways based on the proposed functions of the genes in a subgroup. We supplemented mutants of behavior group B that includes Pink1KO-WS, park1/Δ21, anneKO-WS/+, auxilin KO-WS/+ and nutcrackerKO-WS (Supplementary Fig.13a), with Coenzyme Q10, an electron carrier of the mitochondrial respiratory chain previously shown to alleviate pink1 phenotypes67. Indeed, the reduced TH positive area of these mutants is significantly rescued; likewise, the SING score of Pink1KO-WS and nutcrackerKO-WS is also significantly improved by Q10 (we did not test the other mutants - Fig.4a,b, Supplementary Fig.13b). In contrast, Q10 does not rescue the reduced TH positive area of mutants in behavior group A1 and A2 (dark blue and medium blue) and the SING score of iPLA2-VIAKO-WS and DJ-1a/bKO-WS (Fig.4a-b, Supplementary Fig.13b), in line with the idea of group-specific rescue mechanisms.

We next targeted endo-lysosomal trafficking in behavior group A68,69. This group comprises the trafficking proteins Rab39, Synj1 and Vps13 and the retromer regulators Vps35 and iPLA2-VIA (Supplementary Fig.13a). Defects caused by mutations in iPLA2-VIA were shown previously to be rescued by stabilizing the retromer complex with the chemical chaperone R5570 and we show that also dopaminergic synaptic innervation defects and SING are rescued (Fig.4c,d, Supplementary Fig.13c). Similarly, the dopaminergic defects in other members of group A, such as DJ-1a/bKO-WS, gbaKO, Rab39KO-WS, omiKO-WS, Vps13KO-WS/+, Synj KO-WS/+, and Tango14 KO-WS/+, are also rescued. In contrast, R55 does not rescue the dopaminergic neuron and SING defects of group B mutants (Fig.4c,d, Supplementary Fig.13c). Interestingly, animals of the subgroup A3 (light blue) including VAC14KO-WS/+, VhaM8.9KO-WS/+, Gdh KO-WS/+ and eIF4GKO-WS/+ show a rescue of the reduced TH positive area when supplemented with Q10 and some with R55 (for VhaM8.9KO-WS/+). These genes encode proteins (loosely) involved in “proteostasis/autophagy” and further work is required to identify tools to target this subgroup more specifically. The significant subgroup dependent rescue by Q10 and R55 (Fisher’s exact test, p < 0.001) supports that dopaminergic defects originate from cluster-specific pathomechanisms. Hence, our work delineates at least 2 Parkinsonism subgroups of which one is defined primarily by “mitochondrial-defects” (group B) and another primarily by “endo-lysosomal trafficking-defects” and “proteostasis/autophagy defects” (group A).

Methods
Resource availability

Lead Contact

Further information and requests for resources and reagents are available upon request and should be directed to the lead contact, Patrik Verstreken (patrik.verstreken@kuleuven.be).

Experimental model and subject details

Fly stocks and rearing conditions
The work in this study was overall compliant with institutional ethical standards. Using CRISPR/Cas9-mediated gene editing, we created a comprehensive collection of Drosophila loss of function Parkinsonism models. We selected the 24 fly homologues of genes known to be causative or conveying risk for Parkinsonism in humans and inserted an attP-flanked w+ reporter cassette in each of the genes substituting the first common exon of possible transcripts.
The "control" strain denotes a semi-isogenized w1118 strain that underwent backcrossing to Canton-S for 10 successive generations, resulting in a strain termed Canton-S-w1118. Additionally a Canton-S-w1118 strain was augmented with a transgenic mini-white integration on the X-chromosome by microinjection of pFL44S{w+}104 into the su(Hw)attP8 landing site (BestGene). This modification results in a red-eyed “control” strain named Canton-S-w1118-w+. The entire genome of Canton-S-w1118 was sequenced, revealing the absence of deleterious mutations in homologs of known causative Parkinson's disease genes, as well as in homologs of human genes proximal to GWAS loci. In order to maintain genetic consistency within the experimental framework, all Parkinson's disease KO lines (see below, except for park1/Δ21) were subjected to over ten rounds of backcrossing with the established Canton-S-w1118 control strain. A thorough verification protocol was implemented for all generated fly lines, which included conducting multiple PCR assays followed by sequencing of the PCR products. Detailed information regarding the primers utilized for this verification process can be found in the supplementary file (Supplementary Document 1). For park, the null-mutant fly lines park1 and parkΔ21 from the Bloomington Drosophila Stock Center were used 26,27.
Flies were maintained on standard corn meal and molasses food. For experiments, flies were raised in parallel on the same batch of food, which was exchanged every three to four days, in a temperature and light-controlled incubator at 25°C and 12-hour light-dark cycle and kept in mixed populations of similar density. Control and experimental genotypes were aged to 6±1, 26±2 and 42±2 days (d) as indicated for immunohistochemistry and behavior assays except for some Parkinsonism mutants: nutcrackerKO-LWS were tested at 15 d, Pink1KO-WS at 22±2 d, iPLA2-VIAKO-WS and gbaKO mutants at 25±1 d due to their shorter survival.
For Q10 and R55 treatments male flies of control and experimental genotypes were collected at 1-3 days post eclosion and aged on treatment food, which was renewed every two to three days. On the day of the experiment flies were flipped on fresh treatment food at ZT zero to one. All treatment food was prepared every one to two weeks and kept at 4°C in darkness. Q10 (Merck) was prepared as a 50 mM stock solution in ethanol (99.8%), aliquoted in single-use portions and stored at -20°C. R55 (Calbiochem/Millipore) was dissolved in sterile distilled water as a 12 mM stock solution, aliquoted in single-use portions and stored at -20°C. Q10 and R55 with their respective solvent controls were added to fresh fly food heated up to around 40°C to a final concentration of 1 mM for Q10 and 12 µM for R55.
Both L-Dopa (Merck) and D-Dopa (Merck) were prepared as a 15.2 mM stock in sterile distilled water on ice and stored at -80°C in single-use aliquots. Both L-Dopa and D-Dopa as well as an equal volume of the solvent were added to fresh fly food as described above to yield a 1 mM final concentration. Aged male flies were selected and transferred to L-/D-Dopa food 10 days prior to the behavior assay.105
The fly stocks and genotypes were used and abbreviated in this study are summarized in Supplementary Table 3.
w[1118]; gmRPL11/CyO Tb[1] was a gift from Nicolas Tapon and the pdf null mutant, park1 and parkΔ21 were obtained from BDSC.

Generation of Parkinsonism knock-out Drosophila models
Parkinsonism KO D. melanogaster strains were established employing the CRISPR/Cas9 methodology. This involved the utilization of both a donor plasmid (pWhite-STAR) and a tandem gRNA expressing plasmid. The donor plasmids were equipped with an Integrase mediated cassette exchange (IMCE) cassette, which upon integration into the genome, expressed mini-white. Two homology arms facilitating homology-directed repair (HDR) flanked this cassette. HDR is initiated subsequent to the introduction of two double-strand breaks within the DNA encompassing the initial common exon across all possible transcripts. This process culminates in the substitution of the first common exon with the IMCE-cassette. The positioning of homology-directed repair (HDR) was strategically guided by the selection of homology arms, ensuring that the IMCE cassette resided between two non-evolutionarily conserved regions. The donor plasmid also possesses I-SceI endonuclease sites. In the rare event that the donor plasmid integrates fully into the fly genome by the recombination of only one homology arm, we can to correct the integration by expressing ectopically I-SceI in the fly. This enzyme induces a double-strand DNA break at the I-SceI site, initiating a new DNA repair process where the appropriate homology regions can undergo recombination to ultimately yield the correct fly.

Cloning of Tandem gRNA constructs
Unique guide RNAs (gRNAs) targeting all Parkinsonism genes used in this study were identified through http://crispr.mit.edu/ and subsequently introduced into pCFD4: U6:1-gRNA U6:3-gRNA constructs as per established methodologies outlined in Port et al. (2014)106 and detailed in http://www.crisprflydesign.org/wp-content/uploads/2014/06/Cloning-with-pCFD4.pdf.

Cloning of donor plasmids
Homology arms, spanning 1 kilobase (kB) encompassing the initial shared exon of all gene transcripts, were cloned into either pWhite-STAR107 or pWhite-STAR2 (see Cloning of pWhite-STAR2) (Supplementary Fig.1a-b).
Utilizing HiFi DNA Assembly (NEB), four distinct fragments were combined: the mini-white IMCE-cassette, left homology arm, right homology arm, and the plasmid backbone. The homology arms were either synthesized via gBlocks (IDT), which underwent modifications to adhere to IDT specifications, or were amplified through Q5 (NEB) PCR utilizing genomic fly DNA from the target genotype (Canton-S-w1118). The mini-white IMCE-cassette and vector backbone were derived through restriction digestion of pWhite-STAR using AvrII and XhoI enzymes107. To obviate cleavage of the donor plasmid, mutations were introduced to the PAM sequence within the homology arms, predominantly altering NGG to NAA, although alternative modifications were implemented when necessitated, particularly in instances where the PAM sequence resided within coding regions.
To enhance the functionality of pWhite-STAR for backbone integration detection and to prevent the SV40 poly(A)-tail to become functionally silent108, modifications were made to the original pWhite-STAR construct107. Specifically, one FRT site and one I-SceI site were replaced with a 3xP3-eGFP fluorescent eye marker, enabling the selection against backbone integration. In pWhite-STAR2 the homology arms were cloned such that, upon integration into the genome, the mini-white gene is transcribed in the opposite direction to the Parkinsonism gene. In detail, pWhite-STAR was subjected to IXhoI-HindIII digestion, resulting in the excision of one FRT and one I-SceI site. Subsequently HiFi DNA Assembly (NEB) was conducted using two distinct fragments: a 3XP3-eGAFP marker109 using a gBlock (IDT) and an SV40 poly(A)-tail amplified via PCR from the pUAST.attB plasmid110 to create the pWhite-STAR2.
For two fly lines, gbaKO and omiKO-WS the strategy had to be adjusted:
The GBA orthologs in Drosophila consist of two genes, gba1a and gba1b, located adjacent to each other in the fly genome, with a third gene, CG31413, located between them. To simultaneously knock out both Gba1a and Gba1b, we engineered homology arms for the donor plasmid such that, after recombination, a stop codon is introduced in Gba1a, and CG31413 along with the first exon of Gba1b is replaced with our IMCE cassette. In a subsequent step, we reintroduced CG31413 using PhiC31 integrase-mediated cassette exchange with the plasmid pReC-CG31413 to generate gbaKO (gba1a,gba1bKO). The pReC (Rescue Cassette) was generated by linearizing pUC19 with SapI-EcoRI, followed by the insertion of two gBlocks synthesized by IDT. Subsequently, pReC was linearized with XhoI-XbaI and the entire genomic CG31413 region, amplified via PCR from Canton-S-w1118 was cloned into the linearized plasmid trough HiFi DNA Assembly (NEB).
The omiKO-WS stock is homozygous lethal, a phenomenon previously reported by Tain et al. (2009)111: The imprecise P-element excision of G4907 (Genexel Inc.) between omi and mRpL11 disrupted not only the genomic region of omi but also affected the promoter region and transcription start site of mRpL11. Similarly our strategy also affected mRpL11 by generating a null allele for omi. To address this, we utilized the gmRpL11 rescue construct described by Tain et al. (2009)111 in our omiKO-WS flies, which successfully rescued the homozygous lethality of our line.
A comprehensive array of primers, oligos, gBlocks, and fly strains utilized is delineated in Supplementary Document 1.
All plasmids were sequence verified, and plasmid injections were performed by BestGene (BestGene Inc., Chino Hills, CA, US).

Quantitative RT-PCR
For each sample, 10 flies: 5 males and 5 females were collected except for pink1KO-WS, for which we only used hemizygous males, at 1-5 days after eclosion. RNA was extracted using Maxwell RSC instrument and Maxwell RSC RNA kit (Promega) and transcribed using the SuperScript VILO cDNA Synthesis Kit (ThermoFisher Scientific). Subsequently, cDNAs were quantified by qPCR using a LightCycler 480, 480 SYBR Green master mix (Roche) using qPCR primers enlisted in Supplementary Table 4. qPCR primers were designed using the PrimerQuest Tool (Integrated DNA Technologies). mRNA levels were determined using the Δ-Δ-CT method, where Ct values were first normalized to the housekeeping gene Rp49. Gene expression levels were then plotted as a relative gene expression compared to the endogenous Drosophila gene expression in control samples.30

Immunohistochemistry and confocal imaging
Immunohistochemistry of 6±1, 26±2 and 42±2 day old adult brains in at least 2 independent experiments were dissected in ice-cold PBS and fixed for 30 min in freshly prepared 3.7 % paraformaldehyde (in 1x PBS, 0.2% Triton X-100 (PBX)) at room temperature (RT). After 3x washes for 15min in PBX at RT on a shaker, they were incubated in blocking solution (PBX, 10 % normal goat serum (NGS)) for 1 h at RT. Subsequently, the brains were incubated in primary antibodies (rabbit -TH (Millipore Cat# AB152, RRID: AB_390204), 1:200, mouse -DLG (DSHB Cat# 4F3 anti-discs large, RRID: AB_528203), 1:100) in blocking solution at 4°C for 1.5-2 days and washed 3x 15 min in PBX at RT. Secondary antibodies (goat -rabbit IgG Alexa Fluor™ 488 (Life Technologies Cat# A-11034, RRID: AB_2576217), goat -mouse IgG Alexa Fluor™ 555 (Life Technologies Cat# A-21424, RRID: AB_141780), at 1:500) in PBX with 10 % NGS were applied overnight at 4°C and brains were finally washed 3x 15 min in PBT at RT on a shaker and mounted, anterior facing up in RapiClear 1.47 (SUNJin Lab).112 Z-stacks of entire brains were imaged at a Nikon A1R confocal microscope with a 20x (NA0.95) water immersion lens. Acquisition was performed using a galvano scanner, a zoom factor of 1, scan speed 0.5 and line averaging of 2. All images were acquired with a pinhole of 2.3 Airy unit and a resolution of 1024 × 1024. Z-stacks (step intervals of 3 μm) were used in data acquisition and the same image settings were maintained across the genotypes and imaging sessions. Fiji 113 was used for image analysis. Dopaminergic neuron innervation of the mushroom body was quantified as follows: Anti-DLG was used to determine the five z-stacks which include the synaptic region of the mushroom body (MB) lobes. The outline of the MLB (region of interest, ROI) was defined in the sum projection of the five z-planes. Within the selected z-stacks, the area of the anti-TH fluorescence was thresholded (default threshold- using the same threshold for all the 5 z-planes) excluding background signal comparable to the control. The thresholded area within the ROI was quantified in every z-plane, summed and normalized to the area of the MB ROI for every brain individually (TH+ area/MB area). For every experiment, the individual TH+ area/MB area values were normalized to the mean of the control. The maximum age of every genotype is plotted in “aged flies”. For representative images the maximum projection of five z-planes as well as the thresholded middle z-plane is shown.

Startle induced negative geotaxis (SING)
Negative geotaxis was assessed at Zeitgeber ZT three to five in mixed fly populations or in experiments with treatments (L-/D-Dopa, Q10, R55 and their corresponding controls) in groups of male flies as described previously114,115. Flies at indicated ages (Fig. legends) were transferred into the first tube of the apparatus without prior anesthesia and allowed to adjust for some minutes to experimental environment (24±1ºC, 40-60% humidity). They were tapped down five times and allowed to climb up for 30 sec, at which time the flies that reached the upper tube were moved to the next tube. The procedure was repeated four times, the number N of male flies with the corresponding genotype were counted per tube to calculate the SING score, which was normalized to the mean SING score of the control116.
(1) 𝑆𝐼𝑁𝐺 𝑠𝑐𝑜𝑟𝑒 = ((𝑁1*0)+(𝑁2*1)+(𝑁3*2)+(𝑁4*3)+(𝑁5*4))/4(𝑁1+𝑁2+𝑁3+𝑁4+𝑁5)
𝑁𝑘 = number of flies in the kth tube

Seizures – like behavior assay
Seizure-like behavior induction was tested in groups of 5 ± 1 male flies at 6 ± 1 d and 26 ± 2 d post eclosion except for nutcrackerKO-WS and pink1KO-WS, which were tested at 15 d and at 22±2 d, respectively. Flies were transferred into transparent vials, rested for at least 30 min after anesthesia exposure and mechanically stimulated by vortexing the vial for 10 s at maximum intensity36. Numbers of resistant flies and flies with seizure-like behavior (involuntary movements or paralysis) were quantified within 10 s post stimulation.117 Graphs represent the % of flies without seizures. To visualize the convergence of familial Parkinsonism cases reported to have seizures with the Parkinsonism mutants showing seizure-like behavior, we used the R (v3.6.3) and R package VennDiagram (v1.7.3).

Electroretinograms (ERGs)
ERGs were recorded from young flies 4±1 d after eclosion as previously described63,64. To limit pigmentation variability due to locus dependent w+ activity in pigment cells of mutant flies, we crossed all control heterozygous flies with CSw1118w+ to achieve w+/w+ (i.e. wild type red pigment). Moreover, all gene pairs will also have a w+/w+ (one w+ from each mutant). Flies were immobilized on glass microscope slides, by use of double-sided tape. 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. The flies were exposed to darkness for 3 seconds, followed by 1 second of LED light illumination. This was repeated 5 times for each fly. Light-evoked signals were amplified by a DC amplifier and the amplified signal was processed by a data acquisition device (Clampex) and Axosope 10.7, connected to a PC running Clampfit 10.7 software (Molecular Devices). ERG traces were analyzed in IGOR Pro 6.37 (WaveMetrics) using a custom-made macro.118

Fly behavior monitoring
The screening of sleep and activity behavior of young male flies was performed using the fly behavior video recording platform ethoscope119. Briefly, at 2-5 days post eclosion flies were introduced into glass tubes (65 mm length, 5 mm external and 3 mm internal diameter) with standard corn meal and molasses food (sealed with wax) on one end and closed with cotton wool on the other. In incubators at 25°C with 12:12 light-dark conditions, 20 flies per ethoscope were recorded for at least 5 consecutive days at 2 frames per second with infrared light and infrared light filters (transmission > 730 nm, LEE Filters)120. We screened n≥69 animals for every genotype in at least 3 independent experiments. The position of each animal was saved at each time point in SQLite files and subsequently analyzed with R (v3.6.3) and rethomics with adjusted R packages behavr, scopr and sleepr121 (v0.3.99). Using the automatic behavior annotator at a resolution of 10 s, the behavior of individual flies was determined. Thus, in 10 s intervals flies were scored as moving or asleep. Sleep was defined following the 5-min rule, in which immobility bouts longer than 5 min were counted as sleep bouts, including the first 5 min. A fly was scored as immobile if the distance between 2 consecutive frames in the 10 s interval was less than 0.25 mm (corresponding to a velocity correction factor 𝑎 = 0.0042). Flies that died during the monitoring as well as the first day of recording due to habituation were excluded from further analysis. The following parameters were quantified to assess the activity, the circadian rhythm and sleep architecture of individual flies: Within 5 days the mean sleep amount per light condition (sleep fraction/day) and per dark condition (sleep fraction/night) were calculated. For the weighted mean latency to longest bout, the time between ZT 12 (= lights off) and the longest sleep bout was determined per night. The weighted mean bout length (mean bout length/night) as well as the weighted mean number of sleep bouts (nr bout/night) were scored during the dark condition. Morning anticipation were calculated as described previously59. The mean velocity if awake (velocity if awake) was assessed based on the velocity determined within the 10 s interval annotation ‘moving’. The cumulative distance within the 10 s interval was calculated and summed per 24 h to determine the travelled distance (total distance) for each animal per day. Representative sleep traces of genotypes were visualized using the ggetho package in R (v0.3.6).

Pattern recognition using non-negative matrix factorization (NMF)
The pattern recognition for behavioral analysis was done in three steps. First, the (x,y)-coordinates of each individual fly are transformed into a one-dimensional vector that captures the fly’s temporal behavior. Specifically, we counted how often the fly performs a certain type of behavior during each period of the day, aggregated over the period of the experiment. Second, the vectors of all flies are concatenated to form a matrix and non-negative matrix factorization (NMF)51,122 is applied to automatically cluster common behavioral patterns. This method effectively identified shared features in the behavior of (a subset of) all flies. Ultimately, the behavior of individual flies can be described by a weighted combination of these shared features (e.g., waking up early, more activity during daylight periods). The weights are used to compare the behaviors across different flies, enabling the identification of groups of flies exhibiting similar behaviors. The implementation of the NMF method is provided by Scikit-Learn 122 version 1.0.1 for Python 3.9.6.
Constructing vectors of temporal behavior:
For each fly 𝑓, a vector is created that summarizes during which periods of the day the fly performs a certain type of basic behavior. Therefore, we first categorized the behavior of a fly in each individual frame based on the fly’s corrected velocity, as discussed in Geissmann et al. (2017) 119. For each frame 𝑖 sampled at timestamp t, the corrected velocity 𝑣 is computed as:
(2) 𝑣 = 𝑑𝑖 / ((𝑡𝑖−𝑡𝑖−1)∗𝐹𝑃𝑆∗𝑎) = 𝑑𝑖 / 𝑎,
where 𝑑 is the log-transformed Euclidean distance between successive frames. The velocity correction factor 𝑎 = 0.0042 was empirically defined.
Second, we created bins for each minute of the day (24h x 60min = 1440 bins) and aggregated the behavior of each fly in each bin over the duration of the experiment (5 days and nights; excluding the first day of recording due to habituation) obtaining a vector 𝑥 ∈ ℝ1440. Similar to the rethomics based analysis we used the aggregation approach based on “Time Spent Sleeping”, in which the percentage of frames where the basic behavior was “immobile” for 99% of the frames in the preceding 5-minute period is calculated. To promote smoothness in the values of successive one-minute bins, a Gaussian blur123 with 𝜎 = 5 is applied to the vector 𝑥. This process replaced the value of each bin in 𝑥 with a weighted average of itself and its neighboring bins, resulting in the blurred vector 𝑥′ ∈ ℝ1440.
Identifying behavioral patterns:
For each fly and behavior, we now have a single one-dimensional vector 𝑥′ ∈ ℝ1440 that summarizes the behavior of the fly in function of the time of the day. We then construct the matrix 𝑀 = [𝑥 𝑥 ...𝑥 ] that contains as columns the vectors of all 𝑙 flies in our dataset for behavior type 𝑏. Next, we compress matrix 𝑀 by applying NMF, which is a form of principal component analysis where the resulting components contain numbers ≥ 0. This results in two matrices 𝑊 and 𝐻 such that:
(3) 𝑀 ≈ 𝑊𝐻,
where 𝑀 ∈ ℝ1440×𝑙, 𝑊 ∈ ℝ1440×𝑘, and 𝐻 ∈ ℝ𝑘×𝑙. Here, 𝑘 is a user-defined parameter that refers to the number of principal components for behavior type 𝑏.
The columns of 𝑊 are the principal components that represent basic temporal patterns of behavior. The rows of 𝐻 are compressed versions of the vectors in 𝑀. In other words, if the vector 𝑥 was the 𝑖-th column in matrix 𝑀, then the 𝑖-th row of 𝐻 is its compressed vector. By multiplying the compressed vector with the principal component matrix 𝑊, a vector similar to the original can be reconstructed from only 𝑘 features.
To determine the optimal value for 𝑘, we used the elbow method53. This standard approach plots the reconstruction error (explained variance score) between 𝑊 ×𝐻 and the original matrix 𝑀 over different values of k and looks for an “elbow” (i.e., inflection point), which revealed 5 meaningful and distinct patterns in this case (Supplementary Fig.4e). Additionally, we explored the use of NMF for clustering Parkinsonism mutants based on their dominant behavioral pattern. For this, we assigned each mutant to the cluster corresponding to the maximum component in its encoding vector (i.e., cluster(i) = argmaxW ), then computed the Calinski–Harabasz Index124 across different k for the obtained clustering. The score was maximized at k = 2, indicating that two clusters optimally separate the mutants, but k = 5 offered a reasonable balance between more fine-grained clusters and cluster quality (Supplementary Fig. 4f).
Analyzing behavior as a combination of patterns:
Each value in a fly’s compressed vector is interpretable in the sense that its numeric value quantifies how often the fly executes behaviors of type 𝑏 at periods of the day in the temporal group of a specific principal component. Thus, we quantified how strong each pattern is reflected in the behavior of each fly by aggregated the distribution of principal component values for every fly within each genotype group (Supplementary Fig.4g-k).

Correlation, cluster and principal component analysis of behavior screen
Correlation analysis for the SING score and TH+ area/MB area as well as velocity if awake vs sleep features and NMF components were computed based on the mean values of every genotype in GraphPad Prism 10.1.2 (San Diego, USA). Linear regression plots show the 95% confidence intervals (grey area) and R2 and p-values.
For the cluster analysis, the morning anticipation was excluded, as it would be highly correlative with the NMF components. The means of the activity and sleep parameters as well as NMF components were calculated for each genotype, scaled and hierarchically clustered based on correlation distance and mcquitty-linkage with the R pheatmap package (v1.0.12). Similarly, the dataset was subjected to dimension reduction principal component analysis using the prcomp function and R ggbiplot package (v0.6.2). The scree plot describing the proportion of variance explained by each principal component was determined by dividing the squared standard deviations of the principal components by the total variance and the cor2-values were obtained with R corrplot package (v0.92).

Parkinsonism gene expression clustering across cell-types of the fly brain and body
To assess whether Parkinsonism gene expression differences across the fly brain might explain the different sleep and locomotion phenotypes, we re-mined our previous single-cell RNAseq atlas125. We selected only control cells at 5 days of age, since we were interested in potential baseline differences at a young time point, where also our behavioral analysis was carried out. To perform hierarchical clustering, we first normalized the counts and took the natural logarithm +1, followed by taking the average for each cell-type. After scaling we then used the scipy linkage and dendrogram functions with the wpgma method and correlation distance. To analyse PD expression across cell-types of the fly body, the Fly Cell Atlas “body” dataset was downloaded from https://flycellatlas.org/#data126. The log1p normalized counts of each cell-type were first averaged and then clustered as described above.

Genetic interaction analysis
In a trans-heterozygous screen the amplitude of the ERG depolarisation of every heterozygous Parkinsonism mutant (except for pink1KO-WS, chchd2KO-WS and Rab39KO-WS we resorted to hemizygous flies, as their genes are located on the X chromosome) as well as of all 276 gene pairs were recorded and normalized to the control (CSw1118w+) of the individual experiment. To compute the expected depolarization amplitude for every gene pair in a non-interacting model, we multiplied the normalized mean depolarization amplitude of the single mutants 65,66. Next, the genetic interaction (GI) strength was calculated by subtracting the modelled expected depolarization from the individual observed depolarization amplitude of double heterozygous gene pairs. If the expected and the observed amplitudes are the same (GI strength = 0) the genes are not genetically interacting. Other gene pairs exhibit an “aggravating” or “alleviating” genetic interaction when expected and observed values diverge from another (Fig. 3a). Following prior genetic interaction studies we consider a genetic interaction only when |GI strength|≥ 0.3. Gene pairs with scores within -0.3 to +0.3 were treated as no/weak interaction.127,128 The GI profiles of Parkinsonism mutants were hierarchically clustered based on Euclidean distance and mcquitty-linkage with the R pheatmap package (v1.0.12).

Bayesian analysis
Because the normalized mean depolarization amplitude of the single mutants ignores the variance observed in biological replicates, we opted to additionally run a Bayesian analysis approach aiming to quantify the uncertainty of the results by modelling the consistency of the biological replicates. We defined the following random variables to build our Bayesian model
− α observed depolarization of the first single mutant, relative to the mean of the controls of the same day.
− β observed depolarization of the second single mutant, relative to the mean of the controls of the same day.
− ℩ value which models deviations from the expected value under the no-interaction assumption as an additive amount (so deviations from 0 indicate an interaction).
We built our model so that:
(4) α ~ Beta(2, 5) * 1.05
(5) β ~ Beta(2, 5) * 1.05
(6) ℩ ~ Normal(0, 1)
We decided to model α and β following a Beta distribution with parameters 2 and 5 because it has the following characteristics:
− A posterior probability distribution between 0 and 1 (because they represent fractions of the control). The multiplication by 1.05 makes it possible to be slightly above 1, to model cases where minor increases may be possible
− Asymmetry, with a shift to values higher than 0.5, to capture fractions closer to 1
We decided to model ℩ following a Normal distribution centered at 0, because we wanted to impose very little prior knowledge about whether the gene pair interacts, and with a standard deviation of 1 to model noise in the readings.
This way, the contribution to the observed depolarization by each genotype, normalized with the WT control is as follows
• Single mutant 1 = α
• Single mutant 2 = β
• Double mutant = α * β + ℩
Techniques like Markov Chain Monte Carlo (MCMC) sampling allow us to estimate the posterior probability distribution for each of these variables (i.e. an estimate how likely every possible value is that these variables can take), given the observed data. This way, if 0 is not within the area with the most probability for ℩, the model is telling us the observed data is consistent with some sort of interaction between the genes. Our definition of the area with most probability is the 95% High Density Interval (HDI), i.e. the smallest interval of possible values of ℩, which contains 95% of the area under the posterior distribution.
The implementation of the model and the MCMC sampling is provided by the PyMC package version 5.10.4 for Python 3.10.14.

Jaccard similarity and random distribution analysis
To describe the similarity of behavior-defined groups and genetic interaction groups as well as the stability of behavior-defined groups in the absence of one behavioral feature at a time, we used the Jaccard similarity function, which computes the intersection over union (IoU) for any two groups51. The similarity for each combination was calculated and visualized in GraphPad Prism 10.1.2 (San Diego, USA).
The stability of clusters was determined and expressed in [%] by counting the Jaccard similarity index for every subgroup (A1-B2) in categories of 0-0.34, 0.35-0.69 and 0.7-1.
To estimate the significance of the observed similarity, we generated 10000 simulated sets of random groups of genes with the same size distribution as the genetic interaction groups, and computed how frequently a behavior-defined group had a higher similarity with one of these groups than with a real genetic interaction group. The fraction of simulations (out of 10000) where the similarity is higher in the random group constitutes our p-value (the expected probability of a more extreme similarity under the null-hypothesis). The code for this analysis was implemented in R 4.2.3.

Quantification and statistical analysis
GraphPad Prism 10.1.2 (San Diego, USA) was used to determine statistical significance. Datasets were tested for normal distribution using the D’Agostino-Person Omnibus and Shapiro-Wilk normality tests. For non-normally distributed datasets ANOVA Kruskal-Wallis test followed by a Benjamini-Hochberg post hoc test (q = 0.05) was used for multiple datasets. When multiple parameters were compared (genotypes and treatments) a two-way ANOVA was used, followed by a post hoc Tukey test for multiple comparison correction. Significance levels are defined as ∗∗∗p < 0.0001, ∗∗p < 0.01, ∗p < 0.05 and ns, not significant. ‘n’ in the legends indicates the number of animals used and analyzed. In the SING and seizure assay ‘n’ indicates the number of times the experiment was performed with groups of animals. The association between genotype subgroups (A1–B2) and small molecule treatment (Q10 vs R55) was tested in R (version 3.6.3). A contingency table was generated, and significance was assessed with Fisher’s exact test, appropriate for small sample sizes. Data are plotted as mean ± SD or median and IQR. Specifics on the statistical test used for analysis are reported in the figure legends.

Data Availability
The data generated in this study are provided in the Source Data file with this paper. Further data are available upon request.

Code Availability
The code for data analysis can be found at https://github.com/verstrekenlab/drosophila-Parkinsonism-subgroups.129

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Kaempf et al. have performed a targeted screen using Drosophila to systematically study several dozen genes implicated in PD and other genetic/neurodegenerative disorders that cause parkinsonism. They characterize locomotor impairment, recovery from vortexing, dopaminergic neuron loss, sleep, and electroretinograms. They classify genes into several groups based on the pattern of phenotype observed. They also interrogate pairwise genetic interactions to define potential pathways, and test drugs predicted to modulate mitochondrial energetics and vesicular trafficking. Overall, this work provides a valuable compendium of reagents and data for further cross-species genetic dissection of familial forms of parkinsonism, and begins to define genetic subgroups sharing common molecular mechanisms.\n\nPlease consider the following suggestions for improving the manuscript:\n-I believe that the term “parkinsonism” needs additional qualification for use throughout this manuscript. “Parkinsonism” refers to a rather non-specific human motor syndrome—a pattern that usually localizes to basal ganglia dysfunction. This is similar to other human neurologic motor syndromes, such as “spastic paraparesis” or “ataxia”, which usually result from upper motor neuron or cerebellar dysfunction, respectively. Thus, parkinsonism has anatomic localizing value but does not imply a specific molecular etiology. Drugs / toxins, vascular lesions (stroke, bleeds), and even brain tumors can also cause parkinsonism. This manuscript is focused on genetic/neurodegenerative causes of parkinsonism rather than “parkinsonism” more generally. This is where the Drosophila model really is most helpful due to genomic conservation to examine consequences of manipulation of gene homologs for dopaminergic neuronal health. By contrast, there are no basal ganglia in the invertebrate nervous system—making it a poor model to examine higher-level circuit level contributions to parkinsonism. I would therefore recommend that this be better explained in the introduction. Perhaps clarfiy that the most immediate interest is in “Parkinson’s disease and related genetic/familial causes of parkinsonism”. The discussion of non-motor manifestations (Lines 56-57) is also confusing, and perhaps a distraction—consider deleting. These features are not part of “parkinsonism”, nor do they localize to dopaminergic neuronal / basal ganglia dysfunction, but may be seen accompanying parkinsonism in PD. The term is less relevant in thinking about other non-PD forms of parkinsonism. In the discussion, “the specific causes of most idiopathic parkinsonism cases remain unknown” (line 397)—Are the authors intending here to mean idiopathic PD? This would make the most sense. Similarly, I recommend revising the text later in this paragraph (lines 405-409): “…capture most of the molecular complexity underlying all parkinsonism cases”. For clarity, it would best to specify “PD and related genetic causes of parkinsonism”, rather than other causes of parkinsonism, such as cerebrovascular pathology.\n-Very limited information is provided on how genes were selected? Were there specific search criteria or was this based on certain review articles? The diseases under consideration are rather heterogeneous, including both pediatric and adult conditions. Some are considered primary dystonias, in which parkinsonism is a secondary or a more variable feature. Given the importance of gene selection to the design of the study, I would consider including a descriptive table in the main text, perhaps a condensed version of Supp Table 1. I think it is also important to mention what (if any) human genes were excluded due to poor conservation in Drosophila. It is useful to classify the different diseases being considered under the umbrella of genetic causes of parkinsonism, such as Parkinson’s disease (GBA, LRRK2), autosomal recessive juvenile parkinsonism, dystonia-parkinsonism, NBIA, etc. I’m not certain that the distinction of typical vs. atypical in Supp. Table 1 is particularly useful—I would consider most of the genes/disorders being included here as “atypical” causes of parkinsonism, with PD being the most common cause of neurodegenerative parkinsonism, and therefore “typical” (e.g., GBA, LRRK2). Some of the genes included, such as COQ2, NUS1, eiF4G1 and HTRA2/Omi are controversial with limited replication or non-replication. It may be reasonable to include them as part of the screen, but it would be best to note the caveats. If these are indeed not true genetic causes of parkinsonism, it could skew / bias the clustering.\n-Related to the above, I believe it is incorrect to talk about a fly recapitulating “hallmarks of parkinsonism”, since flies do not have basal ganglia, and except for movement speed, other essential features that define parkinsonism (rigidity, tremor, postural reflexes) are of course not assessed. Perhaps note instead that flies recapitulate “dopamine responsive, progressive motor impairment”.\n-Based on Fig1a, it appears that different genes are being tested in heterozygosity versus homozygosity (or hemizygosity), with most being tested as heterozygotes? What is the rationale / justification? Was this based on inheritance model? This does not seem to be the case. For example, ATP13A2 causes recessive parkinsonism, but the fly mutant, anne, is tested in heterozygosity? GBA on the hand usually is associated with PD in heterozygosity, but it appears that homozygous KO is tested. What about embryonic or adult lethality?—this information does not appear to be provided anywhere.\n-lines 160-171: seizure is not a common accompaniment of PD and neurodegenerative parkinsonism. Moreover, “seizure-like behavior” following vortexing of flies is not necessarily representative of the pathophysiology underlying an epileptic seizure—caveats should be noted. Why is this considered a “seizure” assay, rather than recovery from a kind of trauma—head trauma is also a cause of parkinsonism…\n-lines 272-283: It is possible that genes participate in common pathways, yet may not be expressed in the same cells but rather in the same multicellular functional circuits comprised of distinct cell types?? Perhaps consider this as a discussion point.\n-Fig 3b: why was ERG rather than locomotor testing used for the gene-gene interaction analysis? Some additional justification/rationale is needed here? A preponderance of genes have ERG deficits in heterozygosity—is this result expected? It would be interesting to test whether any of the interactions are consistent using the SING assay.\n-lines 317-319: I would recommend softening the claim that behavior and gene interaction analysis using ERG “define the same 3 parkinsonism subgroups”. While modest overlap was found, clearly the classifications also diverge in important ways.\n-340-358: Were additional genes tested with drug feedings for rescue within each group (A, C) besides the examples shown?\n-As the authors are aware, the main pathways / gene clusters, such as mitochondrial quality control and endocytic trafficking are well established, and this provides an opportunity to show how the method applied re-demonstrates known relationships among genes. Thus, the classification of genetic causes of parkinsonism into groups would be most compelling if it immediately allowed novel predictions that could be validated from human data. While it may not be reasonable to include this in the current manuscript, maybe the authors speculate on what type of prediction might be made for consideration in future analyses? Could this guide gene discovery or therapy? One potential future approach to consider in the discussion is to apply the method to candidates from PD GWAS, where dozens of potential genes might be assigned to the clusters defined from Mendelian forms of PD / parkinsonism.\nSome addition minor questions / suggestions:\n-Lines 111-113: the description gives the impression that the SING phenotype may be specific for DA neuronal defects; however, many other CNS (or PNS) lesions can cause SING phenotypes.\n-Fig 1: rather than “young” and “aged”, it would be best to note the exact age ranges in the fig.\n-lines 155-157: perhaps note that only “selected mutants” were tested with H2O controls (Supp Fig 1g)\n-lines 184-185: it may be helpful to many readers to include additional discussion of the sleep disruptions that occur in PD/parkinsonism and what is know about anatomy / pathophysiology to put the screening based on this phenotype into a bit better context.\n-line 352: “We challenged this result”--unclear what this phrase means?\n-lines 359-362: these final sentences seem better suited to discussion. More context and citations are suggested for the claim that group B is related to autophagy/proteostasis.", "reviewer_2": "The manuscript by Kaempf, et al. is extremely impressive and exceptionally significant. The authors generate comparable knockout mutations in 24 genes in Drosophila that are orthologous to genes strongly associated with Parkinson Disease in humans. They then performed unbiased assays investigating a variety of behavioral phenotypes of these mutant flies and using them to cluster the mutations into groups, as well as systematically performing a genetic interaction test among the mutations based on neuronal function. Quite remarkably, the gene clusters suggested by the multidimensional behavioral analysis show very good concordance with the pattern of genetic interactions, suggesting that the two orthogonal approaches are clustering PD risk genes consistently into the same three discrete classes. Finally, a pharmacological test provides an initial suggestion that the classes defined by phenotype and genetic interaction may also serve to stratify individual mutations into classes that are differentially sensitive to particular drug treatments.\n\nA large number and variety of genetic loci have been linked to PD, as is also true of other neurodegenerative diseases. Like those other diseases, while the genetics have highlighted a variety of molecular pathways as potentially being implicated in PD it has been exceedingly difficult to see an underlying order in that diversity of genetic causes and predisposing factors in PD. The current work potentially provides a unique insight into the genetic structure of the disease by giving unbiased, systems-based grounds for pulling together different genetic loci into a small, defined collection of experimentally-motivated groupings. In addition to the value of the specific results, it provides a template that can and should be used with other datasets to test and expand the analysis here.\n\nI have one more significant comment and one or two minor comments\n1. I was disappointed that results of the pharmacological test (Fig 4) were reported for only 4 of the 24 mutations. The concordance of treatment response with gene groupings is certainly suggestive, but the number of genes tested is sufficiently limited to make one worry that this is just happenstance. If possible, it would certainly be preferable to see the results of these pharmacological tests across the whole collection of mutants for at least one of the two assays of Fig 4, if not for both. One doesn’t expect to see perfect agreement, just as the concordance between the behavioral and ERG assays is not perfect and is not expected to be, but as it stands it is not clear how strongly to interpret this result.\n2. Regarding the ERG assay, it sounds like the flies were w- at the genomic locus, with w+ introduced in the knock-in procedure (and therefore with different amounts of w+ activity in pigment cells in different lines, depending on the locus). Is that correct? At the age the assay is done, is the w+/- status of the fly expected to affect the ERG results? It would be worthwhile to address this explicitly in the text for the sake of readers who are not familiar with the assay.\n3. This is a very minor point, but my download of the supplements had a blank page for Suppl Table 2, the values of the genetic interaction screen.", "reviewer_3": "The authors generated a new collection of fruit fly mutants for 24 PD-related genes and used it to group these genes into three subgroups with similar functions and behavioral activities. Since Drosophila has been used very efficiently in previous studies to understand the pathological functions of several PD-related genes, this study can be a good starting point for understanding the molecular mechanisms of all currently known PD-related genes whose molecular mechanisms have not yet been clearly elucidated.\n\nMajor Points:\n1. It is commonly accepted that disrupting the first exon of a gene typically leads to the downregulation of its expression. However, as the techniques the authors employed were novel and these flies were firstly reported in this study, the authors should provide data on the target gene expression in homozygous and heterozygous mutants to confirm that they are indeed true genetic knockouts.\n2. Measuring TH+ area relative to the mushroom body area is a very interesting and effective approach. Additionally, the observation that this index is primarily age-dependent and correlates with motor activity is intriguing. However, it is important to note that the main function of mushroom body in fruit fly is not closely related to the main symptoms of Parkinson’s disease (parkinsonism). Clarifying this point would help readers better understand the relevance of this new readout.\n3. The authors compared parkinsonism phenotypes in Figure 1 using only genetic knockout flies, which is a logical approach. However, in subsequent figures, such as Supplementary Figure 2, where seizure phenotypes are addressed, the authors clustered PD genes using data from both humans and flies. Therefore, it is crucial to consider whether the mutations in PD genes result in loss-of-function or gain-of-function effects. For example, since pathogenic mutations in LRRK2 are known to be gain-of-function, it would be more appropriate to use LRRK2 gain-of-function mutants, such as G2019S mutant, for a more accurate analysis.\n4. The connection between the ERG/sleep phenotypes and MB TH+ area is ambiguous. It would be helpful to include correlation graphs, similar to the correlation graph between SING and TH+ area in Figure 1e, showing the relationship between ERG/sleep and MB TH+ area.\n5. This study further elucidates the genetic interactions among genes within the same subgroup. However, while the authors used chemicals such as Q10 and R55, it is important to examine whether other known genes within the three clusters have functional interactions with the PD-related genes discussed in this manuscript. For example, it would be crucial to demonstrate whether vesicle trafficking-related PD genes, such as Synj1, Rab39, Vps13, or Vps35, functionally interact with other known vesicle trafficking genes. To address this, the authors should consider selecting genes known to be involved in vesicle trafficking, such as sec22, and performing overexpression or knockdown experiments in a Synj1 (or other vesicle trafficking-related PD genes) knockout background.", "reviewer_4": "In this manuscript, Kaempf and co-workers develop a new library of 24 Drosophila models of familial parkinsonism using CRISPR/Cas mediated knockout of homologous genes known to cause parkinsonism in humans. The authors carefully backcrossed their mutant lines for >10 generations to the same control isogenic stock. These genetically well-controlled lines are likely to be a great resource to other related studies. The authors identified three clusters of mutants that overlap in mitochondrial function, retromer/vesicle trafficking and proteostasis/autophagy. The construction of behavior-based clusters with links to common biology of a complex disease is particularly attractive, as it greatly raises the chance that the results can potentially accelerate biomarker discovery and drug development. The results show significant reduction of climbing performance with age as well as reduction of TH positive area for most of the mutants screened. Climbing performance was ameliorated when feeding L-Dopa. Climbing and TH+ area was partially restored after feeding with subgroup specific treatments for mitochondrial dysfunction and endo-lysosomal trafficking.\n\nThe manuscript is well-written and sufficiently rigorous in most places. The studies are conducted with large numbers of animals and the methods are fairly detailed.\n\nA few points of concern are outlined below:\nMajor concerns:\n1. Missing alpha-synuclein staining data: The authors screened the mutants using climbing to measure motor performance and immunostaining with anti-TH (dopaminergic neurons) focusing on MB. Was there alpha-synuclein accumulation? Since these genes have been linked with other neurological disorders, including alpha-synuclein staining would be important for providing more direct conclusion as to these models recapitulating PD specific pathology.\n2. Robustness of the hierarchical classification of mutant sleep patterns: It is not clear why it was necessary to combine NMF features with user-defined behavioral features. How many clusters emerge if only NMF (more unbiased of the two) is used? In using NMF, the use of the “elbow method” is somewhat risky to this referee as this method does not have a rigorous statistical foundation and is widely considered heuristic, at best. The authors may want to consider better alternatives such as those developed by Tibshirani, Walther & Hastie (2001) and/or Calinski and Harabasz (1974). A related question: how sensitive is the clustering to the choice of the number of total features? Lastly, how do we know if the user-defined features are independent of the NMF features since they are used together?\nDid the authors examine locomotor data instead of sleep? Locomotor data may contain more movement details that are “averaged out” in the process of constructing sleep timeseries.\n3. Use of ERG to characterize parkinsonism: It is unclear what is(are) the rationale behind choosing ERGs as a readout for mutants that showed primarily dopaminergic impairment (Figure 3). The major neurotransmitter involved in Drosophila phototransduction is histamine.\nThe ERG signal is a summation of sustained receptor potentials (depolarization) in photoreceptors and the activity of downstream laminar neurons that are receiving synaptic input from the photoreceptors. Therefore, why was only depolarization change measured and ON and OFF transients excluded? Moreover, due to other factors such as pigment cells that affect the initial depolarization response after light stimulation, measurement of receptor potential (depolarization) is commonly measured when it reaches plateau [Belusic, G., ERG in Drosophila, in Electroretinograms. 2011, InTechOpen; Vilinsky I, Johnson KG. Electroretinograms in Drosophila: a robust and genetically accessible electrophysiological system for the undergraduate laboratory. J Undergrad Neurosci Educ. 2012].\nIn addition, why was 4 DAE chosen for ERG readout when other experiments were performed in aged flies given the evidence for normal synaptic innervation in young flies (Figure 1)? The traces shown in Figure 3 do not show significant synaptic dysfunction. The authors conclude that some combinations are alleviating, but this measurement is not tracked with age.\nMoreover, Figure 1 results show significant reduction of climbing performance with age as well as reduction of TH positive area for most of the mutants screened. This data is not consistent with the reduction of depolarization from ERGs in Figure 3. Some mutants such as nutcracker showed significant reduction of TH+ staining and climbing performance but no change in depolarization. As well as others showed no significant reduction in TH positivity or climbing, but significant depolarization. These discrepancies should be discussed in more detail.\n4. Scaled “velocity” in Method details Constructing vectors of temporal behavior: If d=log(pixel distance), using v=10^(d/1000)/0.0042, I get only values of v>240, since 1/0.0042 is the minimum. Is that correct or is there a typo or did I misunderstand? How was this expression for v derived and why was it necessary to re-scale movements? Either way, the 3 speed categories are arbitrary. It is not clear why the authors use rescaled value of velocity (actually speed, if I understand correctly, since velocity is a vector) in Supplementary figure 4 and not the values determined in mm/sec directly from the videos.\nAdditional (less critical) comments:\n1. Fig 1. Panel F second p value seems to have wrong > direction (VAC14).\n2. Fig. 2. How do non-parkinsonism flies, such as pdfKO, cluster with the mutants considered in this study?\n3. In Fig 3c, what are the horizontal dashed lines meant to signify? Could the authors provide more examples of “alleviation” and add examples of “aggravation” as well (pertaining to major concerns# 3)?\n4. How Bayesian analysis is used in ERG analysis is unclear and should be justified better with some intermediate (between raw ERG and final process) results. Why is variance in ERG of particular concern while all other data in this study also come with large variability? Couldn’t single ERG recordings not be used, like single fly behavioral timeseries?\n5. Suppl Fig 2. The authors report seizure -like activity. The methodology used to define a seizure was unclear. Is it based on previously published methods? Are the seizures only occurring after a mechanical trigger (vortex)? Do the flies exhibit any involuntary movements? What about after L-Dopa administration? Supplementary video recordings and an expanded description of the assay would be helpful.\n6. Suppl Fig. 6. By looking at gene expression profiles in the brain, the authors conclude that clustering of gene expression cannot explain clustering of mutants according to behavior. Since PD-important genes are likely expressed widely in the fly, can the authors comment on the possibility of linking genes to behavior by examining expressions outside the brain as well? For instance, the gut (and gut-brain connection) has received particular attention in several Parkinson’s studies.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their thoughtful questions and remarks. We have addressed the comments and below we summarize the most important additions:\n1) Our work revealed different groups of familial Parkinson’s disease and related genetic causes of parkinsonism based on behavioral screening. We conducted several additional analyses to verify the robustness of our clustering. First, we verified that the parameters we included for clustering are independent from one another by conducting a correlation analysis. This shows that the (independent) parameters we used for clustering do not significantly correlate; Second, we removed single parameters from our clustering analysis and find this does not majorly affect the final clustering result; third, we included the pdf mutant. We previously found that parkin, pink1 and pdf mutants show similar sleep phenotypes1 and now find in our unbiased assessment that pdf clusters with pink1, parkin and chchd2 in “group B”. Taken together, these new experiments further support the robustness of our clusters.\n2) We tested if all mutants in each of the groups with a dopaminergic synaptic innervation defect are rescued by Q10 (that boosts mitochondrial function) or R55 (that is a retromer chaperone). We show that the synaptic innervation of dopaminergic neurons of mutants of the group that contains pink1 and parkin (known to affect mitochondria) (group B) are rescued by Q10 but not by R55, and vice versa, that the synaptic innervation of dopaminergic neurons of mutants of the group with iPLA2VIA and vsp35 (known to affect the retromer) (group A1-3) are rescued by R55 and not by Q10. These results are consistent with our conclusion that the genetic space of Parkinson’s disease and related genetic causes can be ‘stratified’ in specific groups each with specific molecular dysfunction.\n3) Some of our mutants showed activity-induced seizure-like behavior. We show this is not rescued by L-DOPA (while the dopaminergic neuron-related phenotypes are rescued).\n4) We included RT-PCR, validating our mutants. We also noted problems and molecular genetic rearrangements in our parkKO-WS mutant that we still do not fully understand. We therefore replaced parkKO-WS throughout the paper by park1/Δ21 (a null mutant; also confirmed by RT-PCR). This includes the initial phenotyping: SING, TH-labeling and L-DOPA responsiveness, as well as the ethoscope recordings, ERG-genetic interactions against the rest of the mutants, clustering and the Q10 and R55 rescue experiments. We find the new parkin mutant -as expected- clusters with pink1 and chchd2 (also thought to affect mitochondria) and it is rescued by Q10 (and L-DOPA), but not by R55.\n5) We expanded the expression analysis of mutants to include single cell sequencing datasets of the entire fly body. Like our analysis of gene expression across brain cell types, we also do not find a correlation between expression pattern and the assignment of mutants to groups in our clusters. In other words: gene expression patterns are likely not sufficient to define the ‘Parkinson’s-disease-and-related-genetic-causes-of-Parkinsonism’ groups.\nBelow we address each reviewer point one-by-one."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revisions, including new experimental data, have resulted in a much improved manuscript. Collectively, a remarkable amount of work is presented, making for an important and valuable contribution. It appears that a technical glitch remains with accessing Supplemental Table 3.", "reviewer_2": "This manuscript provides functional analyses of mutations in Drosophila homologs of 24 genes that have been linked to Parkinson and other motor dysfunctions in humans. They use multiple assays to cluster the phenotypes of these mutants and show that the genes fall into two broad groups (or 3 groups depending on how one sets the cutoffs). Remarkably, orthogonal assays lead to the same basic groupings; even more remarkably, response to pharmacological treatment largely obeys the same groupings. Perhaps most striking, this is true not only of genes that are known to have related functions, but even of genes that seem superficially to have little in common but that co-cluster in the unbiased assays and then are rescued by the same drugs.\n\nThis paper brings an unexpected and long-sought order to what has been a bewildering multiplicity of functions and properties of PD-associated genes. It also offers a whole new way to think about the basic pathogenesis of PD, a meaningful way to segregate pathological genes into functional pathways, and a model for how we can dissect other diseases with complex genetic etiology, including, AD, FTD-ALS, and others. It is of huge significance to the specific field of Parkinson disease, to the larger field of human neurodegenerative disease, and in the broadest sense as a model for dissection of complex genetic processes.\n\nI found the original submission of this manuscript to be quite impressive; the revised version is even better. In my opinion the authors are to be applauded for the clarity and thoroughness with which they addressed the comments of all the reviewers, myself included. I am entirely satisfied by their responses. There were a number of referee requests to expand on details of the clustering methods, which clearly were central to the authors’ conclusions; I found the responses to be sensible and convincing.\n\nOne thing that struck me in reading the comments is that a number of the referee concerns boil down to the question of why assays like ERG and sleep, that seem to be only peripherally related to diagnostic phenotypes of “parkinsonism”, should be appropriate tools for clustering molecular lesions associated with motor dysfunction. However, it seems to me that the success of these measures (as demonstrated by the concordance between clustering by multiple assays and by response to treatment) is one of the most profound results of the entire project. It provides an astonishing demonstration that the disease pathogenesis is rooted, not in the details of dopamine biology per se, but in basic mechanisms of neuronal cell biology. These data provide powerful support for the idea that it is the selective sensitivity of particular neuronal cell types to specific aspects of that basic cell biology that distinguish among the various neurodegenerative diseases (an idea that is also underscored by results in a previous study from this group, Praschberger, 2023). The work by these authors now gives us a way to approach the basis of that selective sensitivity that has been so elusive.", "reviewer_3": "I am fully satisfied with the authors’ responses and additional experiments, and I have no further questions. The authors have appropriately addressed the reviewers’ questions and critiques, and this revision is sufficient for publication.", "reviewer_4": "The investigators have provided a substantially improved revised manuscript. The changes to Figure 3 and its related discussion offer greater clarity regarding both the methodological approach and the conclusions drawn regarding the genetic interaction screen using ERGs. Inclusion of the Calinski-Harabasz index bolsters their earlier claims about data clustering. Revised Supplementary Figures 6 and 8 are helpful. The revised Supplementary Figure 3 along with newly included videos capturing seizure-like behavior is a valuable addition, that enables a more comprehensive assessment of the described phenotype. Moreover, the incorporation of small molecule treatment for all mutants, as detailed in revised Figure 4, provided compelling evidence to support the assigned classifications. I still find description of fly speed calculations unnecessarily complicated but that may be more due to the work of Geissmann et al rather than the current manuscript.\n\nThe revised manuscript is recommended for publication.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their thoughtful questions and remarks. We have addressed the comments and below we summarize the most important additions:\n1) Our work revealed different groups of familial Parkinson’s disease and related genetic causes of parkinsonism based on behavioral screening. We conducted several additional analyses to verify the robustness of our clustering. First, we verified that the parameters we included for clustering are independent from one another by conducting a correlation analysis. This shows that the (independent) parameters we used for clustering do not significantly correlate; Second, we removed single parameters from our clustering analysis and find this does not majorly affect the final clustering result; third, we included the pdf mutant. We previously found that parkin, pink1 and pdf mutants show similar sleep phenotypes1 and now find in our unbiased assessment that pdf clusters with pink1, parkin and chchd2 in “group B”. Taken together, these new experiments further support the robustness of our clusters.\n2) We tested if all mutants in each of the groups with a dopaminergic synaptic innervation defect are rescued by Q10 (that boosts mitochondrial function) or R55 (that is a retromer chaperone). We show that the synaptic innervation of dopaminergic neurons of mutants of the group that contains pink1 and parkin (known to affect mitochondria) (group B) are rescued by Q10 but not by R55, and vice versa, that the synaptic innervation of dopaminergic neurons of mutants of the group with iPLA2VIA and vsp35 (known to affect the retromer) (group A1-3) are rescued by R55 and not by Q10. These results are consistent with our conclusion that the genetic space of Parkinson’s disease and related genetic causes can be ‘stratified’ in specific groups each with specific molecular dysfunction.\n3) Some of our mutants showed activity-induced seizure-like behavior. We show this is not rescued by L-DOPA (while the dopaminergic neuron-related phenotypes are rescued).\n4) We included RT-PCR, validating our mutants. We also noted problems and molecular genetic rearrangements in our parkKO-WS mutant that we still do not fully understand. We therefore replaced parkKO-WS throughout the paper by park1/Δ21 (a null mutant; also confirmed by RT-PCR). This includes the initial phenotyping: SING, TH-labeling and L-DOPA responsiveness, as well as the ethoscope recordings, ERG-genetic interactions against the rest of the mutants, clustering and the Q10 and R55 rescue experiments. We find the new parkin mutant -as expected- clusters with pink1 and chchd2 (also thought to affect mitochondria) and it is rescued by Q10 (and L-DOPA), but not by R55.\n5) We expanded the expression analysis of mutants to include single cell sequencing datasets of the entire fly body. Like our analysis of gene expression across brain cell types, we also do not find a correlation between expression pattern and the assignment of mutants to groups in our clusters. In other words: gene expression patterns are likely not sufficient to define the ‘Parkinson’s-disease-and-related-genetic-causes-of-Parkinsonism’ groups.\nBelow we address each reviewer point one-by-one."}]

Decision letter