Transcriptomic effects of propranolol and primidone converge on molecular pathways relevant to essential tremor

Essential tremor (ET) is one of the most common movement disorders, affecting nearly 5% of individuals over 65 years old. Despite this, few genetic risk loci for ET have been identified. Recent advances in pharmacogenomics have previously been useful to identify disease related molecular targets. Notably, gene expression has proven to be quite successful for the inference of drug response in cell models. We sought to leverage this approach in the context of ET where many patients are responsive to two drugs: propranolol and primidone. In this study, cerebellar DAOY and neural progenitor cells were treated for 5 days with clinical concentrations of propranolol and primidone, after which RNA-sequencing was used to identify convergent differentially expressed genes across treatments. Propranolol was found to affect the expression of genes previously associated with ET and other movement disorders such as TRAPPC11. Pathway enrichment analysis of these convergent drug-targeted genes identified multiple terms related to calcium signaling, endosomal sorting, axon guidance, and neuronal morphology. Furthermore, genes targeted by ET drugs were enriched within cell types having high expression of ET-related genes in both cortical and cerebellar tissues. Altogether, our results highlight potential cellular and molecular mechanisms associated with tremor reduction and identify relevant genetic biomarkers for drug-responsiveness in ET.


INTRODUCTION
Essential tremor (ET) is one of the most common movement disorders 1 affecting around 5% of individuals over 65 years old. The disease causes an 8-12 Hz kinetic tremor that typically affects the upper limbs but can also affect the head, voice, and rarely the lower limbs. Tremor intensity can sometimes increase with age and have a severe impact on activities of daily living. Recent studies aimed at identifying common and rare genetic variants have yielded mixed results, possibly due to clinical heterogeneity thus decreasing power of genetic studies 2 . Only a handful of variants have been identified and even fewer of them were replicated in other studies. Therefore, new approaches are needed, and transcriptomics might yield new insights in the pathophysiology of ET.
Recent studies in psychiatric genetics have successfully used drug effect screens to identify putative disease genes 3,4 . This approach is particularly relevant to diseases that have specific drug-responsive subsets of patients, as is the case with lithium responsive patients in bipolar disorder (BD) 5 . This kind of approach has yet to be used in many drug-responsive neurological disorders such as ET where patients respond to two drugs: propranolol and primidone 6 .
Propranolol and primidone are the most common drug treatments for ET, although both were not initially developed to specifically treat ET. Both are efficient at reducing tremor by about 50% in ET patients 6 . Drug response is variable between patients, with some having a better outcome with either propranolol or primidone. The tremor-reducing effects of propranolol are thought to be due to its dual effect on both muscle spindles in the periphery and CNS effects 7 . Propranolol reduces the excitability of muscles spindles but its effects on neurons in the CNS have yet to be investigated 8 . Evidence demonstrates however that it may act through modifying transcription as shown by its upregulation of SHF transcripts, a gene that was shown to be downregulated in the cerebellum of ET patients 9 . The tremor reducing effects of primidone are thought to be mediated through cation channels in the CNS, but the exact neurons or cells through which this effect is mediated are unknown 8,10 . Therefore, the effects of both drugs on CNS cells have yet to be elucidated.
Studying the effects of tremor-reducing drugs on transcription can inform us on mechanisms that reduce tremors. Furthermore, it is possible that genes that are targeted by both drugs are implicated in ET pathophysiology and could allow for the identification of genes harboring putative ET causing variants. Currently, no ET genotypes have been associated with either propranolol or primidone response. Moreover, few ET rare variants have been validated in multiple cohorts 2 . Therefore, since no robust genetic models for ET exist, we set out to study ET drug responses in wild-type cells that are representative of two brain regions associated with ET: the cortex and cerebellum 11 . Even if these drugs have been repurposed to treat ET and are not specific to the disease, understanding their mechanisms in cells from disease-relevant brain regions might yield valuable information on disease pathophysiology itself.
In this study, we identified convergent transcriptomic targets of primidone and propranolol in cortical neural progenitor cells (NPC) and cerebellar medulloblastoma cells (DAOY). Common cellular pathways affected by both treatments were related to neuronal morphology, axon guidance as well as cell-cell interactions as revealed by co-expression and pathway enrichment analysis. We also found that ET drugs specifically affected the expression of genes intolerant to loss-of-function (LoF) variants, hinting at possible enrichment of such rare LoF variants in ET. Furthermore, with integration of single-cell data, we find that drug-targeted genes are mostly enriched in non-neuronal cell types such as endocytes, astrocytes, and oligodendrocytes in both cortical and cerebellar tissues. Our study identifies new putative ET-and tremor-related genes and informs on the molecular and cellular basis for tremor reduction in ET.

Differential expression following propranolol and primidone treatment
To assess the transcriptomic effect of propranolol and primidone on neuronal and cerebellar cells, NPCs and DAOYs were independently treated with clinically relevant concentrations of both drugs for 5 days. Differential expression was done using the Wald test (WT) in the R package Sleuth 12 (Eq. (1)). Treatment of NPCs with propranolol resulted in 1754 DE genes (Supplementary Table 1) while treatment of NPCs resulted in 1571 DE genes (Supplementary Table 2). Directionality of overall transcriptional effect was widely different between NPCs and DAOYs, with propranolol treatment resulting in mostly underexpression in NPCs and overexpression in DAOYs (Fig. 1a, b). Pearson correlation of propranolol-treated NPCs and DAOYs effectively show a weak negative correlation, indicating transcriptomic effects on different genes (r = −0.35, p val < 2.2E−308, Fig. 1e). This correlation weakens when weighing for the most significant DEGs (r = −0.283, p = 7.1E−214, Fig. 1f). Primidone, on the other hand, had a weak effect on transcription in both DAOYs and NPCs with only 200 DEGs (Supplementary Table 3) and 23 DEGs (Supplementary Table 4) in each, respectively (Fig. 1c, d for volcano plots). In NPCs, propranolol and primidone DEGs were lowly correlated (r = −0.06, p val = 1.6E−11, Fig. 1e) with a weaker weighted correlation (r = −0.021, p val = 2.2E−02, Fig. 1f). Similar weak (weighted and unweighted) correlations are seen between the two drugs in DAOYs (Fig. 1e, f).
ET drug targets converge on genes related to movement disorders and ET Shared effects of propranolol and primidone on specific genes increases the likelihood of these genes being integral to tremor reduction in ET. Therefore, convergence of drug effects on expression was assessed by comparing gene Z-scores from different treatment conditions: convergent drug targets in either DAOYs or NPCs, convergent propranolol or primidone targets in both cell types and convergent targets of both drugs in all cell types. Table 1 shows the  top 35 convergent DEGs across all cells and treatment conditions  (see Supplementary Table 5 for full statistics).
Pathway enrichment analysis of convergent propranolol DEGs (in both cell types) was also performed using g:profileR using genes expressed in both DAOYs and NPCs as background ( Table 2). Pathways known to be affected by propranolol such as HIF-1alpha (q val = 0.001) and regulation of apoptosis (q val = 0.02) were significantly enriched. Much like the co-expression analysis, Reactome terms related to axon guidance were found to be significant, such as RUNX1 transcription (q val = 0.0002), a transcription factor implicated in growth cone guidance of DRG neurons 19 . Interestingly, CaMKK2 signaling pathway was found to be significantly enriched within genes in the propranolol geneset. CAMKK2 encodes a kinase implicated in synapse homeostasis and is also involved in modifying Aβ synaptotoxicity in Alzheimer's disease 20 .
Weighted gene correlation network analysis was also performed to identify co-expression modules associated with combined propranolol/primidone treatment. Module-trait and module correlation heatmaps are shown in Fig. 2. Two modules (cyan and red; corr = 0.74, p val = 0.009; corr = 0.73, p val = 0.01 respectively; Fig. 2a) were found to be significantly associated with treatment in DAOYs and only one module (red; corr = 0.65, p val = 0.03) was significantly associated with NPCs (Fig. 2b). Pathway enrichment analysis of DAOY red module genes found an enrichment of Reactome terms related to RABGAP signaling (q val = 0.009) as well as RUNX1 transcription (q val = 0.02; Table 3). NPC red modules genes were significantly associated with neuronal morphology, axon guidance and neurogenesis ( Table 4).
Correlation of the effects of propranolol and primidone with those of common and rare variants in ET TWAS studies the effect of common SNPs associated with a disease on the expression of genes in different tissues.
We postulated that transcriptomic targets of propranolol and primidone might correlate with the transcriptomic effect of common ET variants. We used TWAS summary statistic from the latest ET GWAS 21 to measure the correlation between TWAS gene Z-scores and convergent drug target Z-scores (across all possible conditions) while controlling for gene length and GC content (Eqs. (2) and (3)). Weak, non-significant correlations between TWAS Z-scores and drug target Z-scores were found across all conditions and all brain tissues (p > 0.05; Fig. 3a). Cerebellar hemispheres and cerebellum tissues, brain regions highly associated with ET, displayed no correlations with convergent drug targets (coeff = −0.0143, p val = 0.549; coeff = −0.000138, p val = 0.994 respectively; Fig. 3b).
We postulated that since propranolol and primidone had a nonsignificant correlation with expression of genes harboring common variants for ET, they might instead act on genes that have rare variants. GnomAD recently published observed/ expected (o/e) LoF scores for all protein coding genes in the genome. These scores inform on the tolerance of genes to rare LoF variants, with genes with a higher frequency of observed to expected LoF variants being more tolerant to mutations. Figure 3c shows the distribution of LoF scores of drug DEGs compared to all protein coding genes passing the initial DE QC. Drug targets displayed a significantly lower o/e score median (n = 256, median = 0.18) than all protein coding genes (n = 11,188, median = 0.36; W = 1,727,520, p val = 1.50E−10) using a Wilcoxon unpaired test. Interestingly, when looking at Log2FC direction (Fig. 3d), upregulated genes (n = 194) had a significantly lower o/e score median (median = 0.15, W = 1,361,482, p val = 2.917E−12) than all protein coding genes whilst no significant difference was found between o/e scores medians of downregulated genes (n = 71) and all protein coding genes (median = 0.35, W = 417,126, p value = 0.3246) using a Wilcoxon unpaired test. Thus, propranolol and primidone increased expression of mutationally constrained genes in cultured DAOYs and NPCs.
Single-cell enrichment of propranolol and primidone-targeted genes Our current understanding of CNS cell types affected in ET is still very limited. Enrichment of disease related genes can indirectly inform on potential cell types implicated in disease pathophysiology 22 . We first sought to assess the enrichment of ET genes discovered through familial linkage, whole-exome, GWAS and transcriptomic studies in cell types of the adult cerebellum and cerebral cortex (Figs. 4 and 5 and Supplementary Tables 17-18). Enrichment Z-scores per cell type for ET genes as well as drug DEGs were calculated based on average normalized expression in singlenucleus cerebellum data from Lake et al. 23 and cortical single-cell Smart-seq data from the Allen Brain Institute. In the cerebellum, ET genes were mostly enriched in astrocytes (enrichment z-score = 3.11, q value = 0.021; Fig. 4a, b). In the cortex, the strongest enrichments of ET genes were found in oligodendrocyte progenitor cells (OPCs; z-score = 3.55) and L3-L5 excitatory neurons with the most significant neuronal cell type being the FEZF2-, DYRKexpressing pyramidal neurons of cortical layer V (z-score = 3.28, q val = 0.0068; Fig. 4c). Significant enrichment was also found in L1 MTG1 astrocytes (z-score = 3.13, q val = 0.0090).

DISCUSSION
Understanding the cellular and molecular mechanisms behind drug treatments can inform on disease pathophysiology. In this study, we sought to investigate the transcriptomic effects of first line treatments for ET in cerebellar DAOY cells as well as NPCs, to gain insight on potential disease related genes. Overall, weak Pearson correlations were observed between the same treatment in different cells indicating that the drugs might have completely different effects on genes and pathways in the cortex compared to the cerebellum. Nonetheless, 265 genes were found to be convergent across both treatment and cell types. Indeed, we found that propranolol and primidone affected expression of multiple genes related to movement disorders and ET. Notably, TRAPPC11, whose expression was previously shown to be altered in ET cerebellar cortex and is also involved in protein trafficking 9 .
Other genes related to endosomal trafficking were found to be differentially expressed after propranolol treatment, such as MYO1E and SYNJ1. Convergent DEGs also displayed an enrichment of genes related to the ESCRT complex, known to be a pillar of endosomal trafficking in neurons. These findings potentially increase the likelihood of endosomal trafficking being altered in ET and possibly partly restored through transcriptomic effects of propranolol.
Axon guidance was previously associated with ET in several studies 2,9,11,24,25 . Bulk-RNA-sequencing of cerebellar cortex and dentate nucleus of ET patients showed a significant enrichment of axon guidance genes 9 . Hallmark axon guidance genes such as ROBO1 (z-score = 5.87, q val = 1.88E−06) and NEO1 (z-score = 4.01, q val = 5.04E−03) were both found to have increased expression following drug treatment. NEO1 (and its paralog DCC), which binds netrin-1, is implicated in cell-cell adhesions, mostly between axons and oligodendrocytes, as well as cell-extracellular matrix adhesions. Netrin-1 also acts on dendrite arborization, increasing connections in excitatory synapses 26 . Interestingly, NEO1 protein remains expressed in Purkinje cells of the adult cerebellum (GTEx V8). Thus, the post-developmental role of axon guidance signaling pathways is to maintain adhesions and important synaptic connections between cells. This might be an important process by which ET tremorolytic drugs diminish tremor. These findings on axon guidance are concordant with other Reactome/GO-terms found to be enriched amongst DEGs, most notably semaphorin interactions, cadherin binding, and actin cytoskeleton reorganization. Purkinje cell axons in ET patients have shown accumulations of disordered neurofilaments ("axonal torpedoes") leading to abnormal axonal morphologies 11 . This process is thought to either be part of a neurodegenerative cascade or a response to neurodegeneration. Moreover, decreased neuronal density was observed in multiple brain regions of ET patients, most notably the inferior cerebellar peduncles through which afferent axons from the brainstem nuclei pass in order to reach the cerebellar cortex 27 . Our findings therefore provide additional support for the involvement of axon guidance molecules in ET pathophysiology.
We also identified the CaMKK2 signaling pathway as significantly enriched in propranolol DEGs in DAOYs and NPCs. CaMKK2 exacerbates amyloid-b synaptotoxicity in Alzheimer's disease through Tau protein phosphorylation by AMPK 20 . This pathway is sensitive to cellular calcium intake, which was shown to be affected at the transcriptome level by both propranolol and primidone. Both Tau protein and amyloid-beta abnormalities have been observed in ET cerebellar tissues, with multiple findings pointing toward protein aggregation being a hallmark of the disease 28,29 . Propranolol affecting transcription of genes implicated in both CAMKK2 and Ca 2+ signaling pathways might imply that ET drugs could reduce aggregate-induced neurotoxicity.
Convergent drug DEGs did not correlate with transcriptomic effects of common ET variants (TWAS DEGs). Moreover, propranolol and primidone DEGs displayed weak non-significant correlations with gene expression in the cerebellum of ET patients, the principal brain region affected in this disorder 1 . There are several possible explanations for these results. The relatively underpowered state (for a common disease) of the current ET GWAS might not capture the effects of common variation on transcription, in part explaining the absence of correlation with drug DEGs. Moreover, the lack of good cell models for cerebellar neurons as well as the neurodevelopmental state of NPCs also impair adequate comparisons between TWAS statistics and drug DEGs presented in this study.
Convergent drug DEGs are significantly more likely to be genes predicted to be intolerant to LoF variants. Mutationally constrained genes are more likely to be essential for cell homeostasis  and survival and thus more likely to be implicated in disease when affected by LoF mutations 30 . Given that both ET drugs converged on these genes in multiple cell types increases the likelihood that these genes harbor rare variants associated with ET. Upregulated DEGs were found to be significantly less tolerant than all protein coding genes while downregulated DEGs were as tolerant as all protein coding genes. These genes could be good candidates for future targeted sequencing, especially within propranolol and primidone responsive cohorts. Identifying cell types affected in ET remains difficult. Several conflicting studies have tried to identify specific pathological morphologies in post-mortem cerebellum of ET patients, most notably in Purkinje cells, yet no defining histopathological markers have been found 11 . Here we sought to identify the relevant ET cell types by assessing the enrichment of variant-harboring ET genes within single cells in cerebellar and cortical tissues. Expression of ET genes were mostly enriched within L3-L5 excitatory neurons in the cerebral cortex, more specifically FEZF2 L5 glutamatergic pyramidal neurons 31 . These neurons originate in the primary motor cortex (M1) and form the corticospinal tract that projects to lower motor neurons, which controls conscious movements. These neurons are influenced by multiple cortico-cortical pathways but also input from the cerebellothalamic tract, crucial for movement coordination. The primary motor cortex has previously been shown to be important for tremor generation in ET as subdural stimulation of M1 can reduce tremor intensity in patients 32 . Moreover, propranolol-targeted genes were mostly enriched in VIP-expressing inhibitory neurons of L1-L3. These neurons are known to inhibit motor neurons through different cortical pathways 33 . The enrichment of ET genes within M1 pyramidal neurons coupled with the enrichment of ET drug genes in motor neuron-inhibiting cells does suggest new potential cellular mechanisms through which tremor generation (and/or reduction) occurs in ET.
In the cerebellum, both ET genes and convergent drug DEGs were significantly enriched within astrocytes in the cerebellum. This somewhat contradicts previous histopathological findings postulating that Purkinje cells were the defining cell type in ET pathophysiology. Not much is known about the role of astrocytes in ET but based on other neurodegenerative diseases, it could be argued that they may play an important role in the onset or development of the disease 11 . Oligodendrocytes, whose dysfunction contributes to numerous other neurological diseases, also showed an enrichment of propranolol and primidone-targeted genes. Both astrocytes and oligodendrocytes might be targeted by ET drugs to reduce tremor since non-neuronal cell types are known to be involved in neurodegeneration in numerous diseases 34 . The lack of single-cell data on ET tissues is a limitation in the study of this disease but our results highlight a possible role for non-neuronal cells in the cerebellum in ET.
This study has a number of limitations. Propranolol and primidone are known to act on cell excitability and this effect was postulated as being important for tremor reduction in ET. Given that DAOYs and NPCs are non-excitable, it is very hard to assess the electrophysiological effects of these drugs in these cells. Moreover, the electrophysiological effects of drugs on cells are known to influence transcription 35 . This might explain why primidone had such a mild effect on transcription in both DAOYs and NPCs. The lack of transcriptomic effects of primidone might also be related to the low expression of certain TRP channels by both NPCs and DAOYs which are predicted to be affected by the drug in the context of ET. Contrary to propranolol which acts on Fig. 3 Effects of ET drugs on common and rare variants. a Correlation heatmap of ET TWAS gene Z-scores in different brain tissues and drug effect gene Z-scores from different meta-analysis conditions. Values indicate Z-score regression coefficient from linear model. b Correlation plot of TWAS gene Z-scores from cerebellar tissue and convergent primidone and propranolol gene Z-scores across DAOYs and NPCs. c Line histogram displaying the distribution of O/E LOEUF scores from all protein coding genes (blue) and convergent DEGs (red) following drug treatment. O/E scores were directly transformed to percentages (ex. 0.25 as 25%) with scores over 10 counted as 100%. d Violin plots of O/E LOEUF scores for upregulated DEGs (yellow), downregulated DEGs (red) and non-significant DEGs (green). Boxplot elements: center line = median; upper and lower hinges = 1st and 3rd quartiles respectively; whiskers = mean ± interquartile range × 1.5.
GPCRs that mediate multiple effects on gene expression, primidone might convey its tremor reducing actions by shortterm effects on cell excitability. These effects would evade our transcriptomic screen performed over the span of multiple days. Moreover, cells used in this study do not represent the complete range of cell types in the cortex and cerebellum. NPCs do not completely replicate neuronal expression and do have a more neurodevelopmental transcriptomic state. DAOYs, on the other hand, are derived from cancerous cells and do have dysregulated expression of genes related to cell division and cell growth. Moreover, these cells do not accurately replicate the disease-state present in cerebellar and cortical neurons and other cells in ET patients. These cells are therefore not specific cellular ET models but are nonetheless helpful in understanding the effects of ET drugs in an ET-related cellular context. This study only serves as an ET drug effect screen and remains a steppingstone for more indepth functional studies, leveraging better ET models such as patient derived iPSCs.
Our study identifies multiple cellular and molecular pathways implicated in ET pathophysiology and tremor reduction by both propranolol and primidone. Our findings also suggest a role for genes harboring potentially rare, deleterious variants associated with ET. Targeted sequencing of these convergent drug genes in case-control cohorts could help to confirm or infirm this hypothesis. These genes could also be used as biomarkers for propranolol treatment in responsive ET patients. Our results also identify several cell types involved in ET in both cerebellar and cortical tissues, as well as cells potentially affected by propranolol and primidone through which tremor might be reduced. We believe that this relatively novel paradigm to study pharmacogenomics could be leveraged to repurpose other drugs to treat ET. Moreover, this approach could be used in other diseases to understand the biological effects of drugs with unknown mechanisms of action. Future studies will be needed to further identify the transcriptomic and electrophysiological effects of both drugs in ET, possibly using more representative neuronal models such as iPSC-derived Purkinje cells, non-neuronal cell types as well as motor neurons.

Cell culture and drug treatment
The NPC line is from a healthy individual who is a Caucasian male of 62 years old and provided written informed consent. iPSCs were derived from fibroblasts and all the quality control criteria for validation of the iPSC line reprogramming and integrity have been performed. iPSC colonies were then cultured on Matrigel-coated dishes (BD Biosciences) using mTeSR1 medium (StemCell Technologies). Embryoid bodies were formed by mechanical dissociation of iPSC colonies to induce neural induction, using collagenase, and plating onto low-adherence dishes in STEMdiff™ Neural Induction Medium + SMADi (StemCell Technologies). Embryoid bodies were growing and maintained for 20 days in STEMdiff™ Neural Induction Medium + SMADi (StemCell Technologies). To obtain NPC, embryoid bodies were plated onto polyornithine/laminin (Sigma)-coated dishes in DMEM/F12 plus N2 and B27. Rosettes were manually collected and dissociated with accutase (Chemicon) after 1 week and plated onto Polyornithine/laminin-coated dishes in NPC media (DMEM/F12, 1× N2, 1× B27 (Invitrogen), 1 μg ml −1 laminin and 20 ng ml −1 FGF2 (Invitrogen)). The NPCs were then stained for the neural precursor marker (Nestin) and the pluripotency marker (Sox2) to confirm the NPC state. DAOYs (ATCC) cells were cultured as previously described 9 in Eagle's Minimum Essential Media (ATCC) supplemented with 10% FBS and 5% penicillin/streptomycin (Gibco). DAOYs and NPCs were treated for 5 days with 20 ng/ml of propranolol or 5 μg/ml of primidone (n = 3 per treatment/cell line). H2Oor DMSO (0.023%)-treated cells were used as controls for propranolol and primidone, respectively. Drug concentrations were chosen based on previous studies that tested efficient tremor-reducing serum levels of propranolol and primidone in ET patients 8,36 . A kill curve was used to determine lethal drug concentrations for DAOY cells and NPCs in culture (Supplementary Tables 10-12 and Supplementary Figs. 1 and 2). hexamer cDNA generation (New England Biolabs). Samples were sequenced on the Illumina NovaSeq6000 platform (150 bp paired-end reads, 150 M reads). FASTQ files were pseudo-aligned to the Ensembl v102 annotation of the human genome using Salmon v1.4.0 37 . Gene-level differential expression analysis was done using the R package Sleuth 12 . Only genes with a minimum of 10 scaled reads per base in 90% of samples were kept to filter out low-count genes. Cell types and treatments were analyzed separately using the WT. The full model for the WT was:

RNA-sequencing and differential expression analysis
Differentially expressed genes ðDEGÞ $ plate þ buffer þ treatment (1) MA plots and p value histograms displayed expected distributions (Supplementary Figs. 3 and 4). Meta-analysis of gene Z-scores was performed to analyze convergent DEG across cell types and treatments. Briefly, Z-scores for each gene were calculated and then summed across different combinations of cell types and treatments using Stouffer's Z method 38 .
Multiple analyses were performed notably propranolol specific effect across cell types (labeled "prop"; Supplementary Table 4), primidone effect across cell types ("prim"; Supplementary Table 5), convergent propranolol and primidone effect in each cell type ("daoy" and "npc"; Supplementary Tables 8  and 9 respectively) and convergent primidone and propranolol effects across both cell types ("all"; Supplementary Table 4). False discovery rate was controlled for using the Benjamini-Hochberg procedure (q value threshold <0.05). TRAPPC11 Log2FC was validated in both propranolol-treated DAOYs and NPCs (Supplementary Table 13). At least three DEGs with highest Log2FC per condition were validated using TaqMan qPCR probes as well as 4 top convergent DEGs with 3 out of 4 Log2FCs in the same direction (ROBO1, FAT4, ZMIZ1 and MAP1B) (Supplementary Table 13).
WGCNA, co-expression and pathway enrichment WGCNA was done using the R package 39 . DAOY and NPC sequencing results were analyzed separately, merging both primidone and propranolol treatments in the analysis. Normalized TPM values obtained from Sleuth ("sleuth_to_matrix") were used for the analysis. To filter out noisy lowcount genes, only genes with a minimum of 10 TPM in 47% of samples were kept, for a final list of 8549 genes in DAOYs and 9260 genes for NPC. Two outlier samples ("DAOY_PRIM_03" and "NPC_PRIM_02") were removed from the analysis based on sample clustering dendrogram for a final 22 samples in our WGCNA analysis. Fisher's exact test was used to calculate gene-module p values. Co-expression analysis was performed using GeneNetwork2.0 18 . Pathway enrichment analysis was done using the gprofileR R package 40 . Briefly, gene-lists were made from convergent DEGs across multiple conditions (both drugs in DAOYs or NPCs, propranolol or primidone in both cells, both drugs in both cells). Custom background used in gprofiler comprised genes expressed in either DAOYs, NPCs or both when pertinent. The g:SCS algorithm was used for multiple testing correction (q value threshold <0.1).

Correlation with ET TWAS summary statistics
ET transcriptome wide-association studies (TWAS) summary statistics were obtained from Liao et al. 21 . A linear model was used to measure the strength of association between gene-level drug Z-scores and TWAS Zscores, controlling for gene length and gene GC content ("lm" function in R). Weighted Z-scores were also used to account for significance of effect. The formula used were: And for the weighted Z-score analysis, given by: Both an unweighted and weighted Z-score test was used to assess the correlation with TWAS summary statistics. The weighted test allows to weigh Z-scores by their significance, akin to subsetting for only those with FDR < 0.05 but retains direction and size effect information (upregulated (positive Z-scores) or downregulated (negative Z-scores) gene expression). For example, a small non-significant Z-score (Z = 1), will result in a small Zscore after weighing (e.g., Z * (absolute Z) = 1*1 = 1). A large Z-score (Z = 4) will result in a larger significant Z-score after weighing (e.g., Z * (absolute Z) = 4 *4 = 16). This permits a more holistic approach to testing for correlations between DEGs across TWAS and the drug screen as opposed to subsetting DEGs based on direction and significance.
Association p values were corrected for multiple testing using Benjamini-Hochberg (q value threshold <0.05).

Single-cell enrichment analysis
A one sample Z-test was used to test enrichment of drug-targeted genes as described previously 22 . An ET gene-set was curated from genes associated with ET from linkage, whole-exome, GWAS and transcriptomic studies 2,9 . Drug gene-sets were made from convergent DEGs (FDR < 0.05) across different conditions (DAOY, NPC, propranolol, primidone, all conditions). Adult cerebellum single-nucleus RNA-sequencing data was obtained from Lake et al. (GEO accession: GSE97930) 23 . Average cell counts per cell type were obtained using Seurat v4.0.1 41 . Trimmed means per cell type from adult cortex single-cell RNA-sequencing were obtained from the Allen Brain Atlas Smart-seq multiple cortical regions dataset 42 . To account for drop-out rates and reduce zero-inflation of the single-cell count matrices, low average count genes were filtered out in both cerebellum (<0.5 counts in 7/10 cell types) and cortex (<1 count in 85/121 cell types). Single sample Z-tests were used to obtain cell type specific enrichment Zscores: Z À score ¼ Mean geneset counts À Mean cell type expression counts Geneset ¼ standard deviation ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Number of genes in geneset p (4)

Loss-of-function analysis
The distribution of mutational constraint scores for drug DEGs was assessed using pLoF o/e ratio scores obtained from gnomAD 30 . pLof scores for convergent genes across all conditions with q val <0.05 were compared all protein coding genes passing QC from the Sleuth differential expression analysis. To account for coding sequence length and gene GC percentage, propensity score matching with replacement was used (matchIT package in R 43 ) to measure pLoF score distribution differences between DE drug genes and all protein coding genes included in the meta-analysis. Nearest neighbor matching with the maximum number of matches (ratio = 1:43) between non-DEGs and DEGs was used. A Wilcoxon unpaired test was done on the matched data. The same methods were used to assess pLoF score differences of upregulated (match ratio = 1:57) and downregulated (match ratio = 1:178) DEGs with all protein coding genes.

Statistical tests
No statistical methods were used to determine sample sizes prior to experiments. A randomized layout was used when treating both cell types with propranolol and primidone as well as during RNA-sequencing. Statistical tests were performed in R version 4.0.2 and Rstudio version 1.4.11061. Differential expression analysis was done using the WT in Sleuth. Fisher's exact test was used for WGCNA as well as pathway enrichment analysis using g:profileR. A Wilcoxon unpaired test was used to compare Fig. 5 Single-cell enrichment of drug DEGs in cerebellar and cortical tissues. a Single-cell enrichment Z-score heatmap of convergent propranolol/primidone DEGs in adult cerebellar tissue. Rows represent DEGs; columns indicate cell types; legend color scheme is based on enrichment z-score direction. b Violin plot of average expression per cerebellar cell type of convergent propranolol/primidone DEGs. Boxplot elements: center line = median; upper and lower hinges = 1st and 3rd quartiles respectively; whiskers = mean ± interquartile range × 1.5. c Single-cell enrichment Z-score heatmap of convergent propranolol/primidone DEGs in adult cortical tissue. Rows represent cell types; columns indicate DEGs; legend color scheme is based on enrichment Z-score direction. d Enrichment Z-score heatmap of DEGs gene-sets from different conditions (see below for abbreviations) in single-cell data from adult cortex. Rows represent cell types; columns represent condition gene-sets. e Enrichment Z-score heatmap of DEGs gene-sets from different conditions in single-nucleus sequencing data from adult cerebellar tissue. Rows indicate condition gene-sets; columns indicate cerebellar cell types. ET ET-related genes, prop convergent propranolol DEGs in both cell types, prim convergent primidone DEGs in both cell types, DAOY convergent propranolol and primidone DEGs in DAOY cells only, NPC convergent propranolol and primidone DEGs in NPCs only, all convergent propranolol and primidone DEGs in both cell type.