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 Sleuth12 (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).

Fig. 1: Correlation between DAOYs and NPCs treated with propranolol and primidone.

Volcano plots of propranolol-treated NPCs a and DAOYs b as well as primidone-treated NPCs c and DAOYs d. Blue lines indicate −0.5- and 0.5-Log2FC changes. Red lines indicate q value significance threshold (0.05). e Unweighted Pearson correlations between DEGs z-scores from different conditions of treatment and cell types. f Weighted Pearson correlations between DEGs z-scores from different conditions of treatment and cell types.

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

Table 1 Top 35 convergent DEGs across all conditions.

Across DAOYs and NPCs, 788 significant convergent DEGs were found with propranolol treatment (Supplementary Table 6), 36 convergent DEGs following primidone treatment (Supplementary Table 7) and 265 convergent DEGs across all conditions (Table 1 and Supplementary table 5 for full list). In total, 210 propranolol and 12 primidone-specific convergent DEGs were also found to be convergent across both treatments (e.g., BRD2, MYO1E, ROBO1, etc; Table 1) Propranolol increased expression of TRAPPC11, a trafficking protein previously associated with ET9, in DAOYs (log2FC = 0.98, q val = 5.32E−27) and this gene was also found to be convergently affected across both cell lines (Stouffer’s Z-score = 5.41, q val = 5.87E−06; Supplementary Table 6). Propranolol also decreased expression of G3BP1 in NPCs (Log2FC = −1.22, q val = 3.20E−41) which encodes a protein implicated in stress granule formation and is known to affect axonal mRNA translation as well as nerve regeneration13. This effect was also found to be convergent across both NPCs and DAOYs (Stouffer’s Z-score = −9.07, q val = 7.84E−17). BRD2, a transcription factor previously implicated with epilepsy, was convergently upregulated following propranolol treatment in both cells (DAOY Log2FC = 2.57, q val = 1.17E−230; NPC Log2FC = 0.360, q val = 0.007; Stouffer’s Z-score = 21.13, q val = 4.56E−95). NONO, a gene harboring a splicing variant known to cause X-linked intellectual deficiency with intentional tremor, was found to be upregulated in DAOYs treated with propranolol (Log2FC = 1.23, q val = 3.05E−60)14. Primidone, in NPCs, upregulated VCAM1 (Log2FC = 0.69, q val = 1.77E−08) and this effect was found to be convergent across both cell lines (Stouffer’s Z-score = 5.53, q value = 1.29E−04). VCAM1 is a gene implicated in axonal myelination by oligodendrocytes15. GIPC1 was also found to be convergently downregulated following primidone treatment when leveraging effects in both cell types (DAOY Log2FC = −0.534, q val = 0.004; NPC Log2FC = −0.360, q val = 0.137; Stouffer’s Z-score = −5.46, q val = 1.42E−04). GIPC1 is a known interactor of DRD3 which has previously been associated with ET and Parkinson’s (PD)2,16,17.

Propranolol and primidone act on pathways related to neuronal survival as well as axon guidance

Following the identification of convergent DEGs across treatments, we aimed to identify molecular pathways affected by propranolol and primidone in DAOYs and NPCs. Co-expression enrichment analysis (using GeneNetwork2.018) for convergent DEGs across all conditions showed that Reactome terms related to GPCR signaling (q val = 1.12E−19), axon guidance (q val = 1.68E−08), Semaphorin interactions (q val = 3.24E−13) and VEGF signaling (q val = 2.23E−08) were significantly enriched within the convergent genesets (Supplementary Table 14). Furthermore, Ca2+ signaling (q val = 4.67E−07) and voltage-gated potassium channels (q val = 4.64E−06) were also found to be significantly enriched. Interestingly, GO:cellular components significant terms were mostly related to cell:cell or cell:extracellular matrix interactions as well as axon guidance such as lamellipodium (q val = 4.47E−13), filopodium (q val = 3.54E−11), focal adhesion (q val = 4.70E−11) and growth cone (q val = 1.04E−09)(Supplementary Table 16).

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 neurons19. 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 disease20.

Table 2 Pathway enrichment for convergent propranolol DEGs in both DAOYs and NPCs.

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

Fig. 2: Co-expression gene modules for convergent propranolol and primidone targets.

a Module-treatment (propranolol/primidone) and -buffer (H2O/DMSO; control) correlation heatmaps for DAOYs. b Module-treatment (propranolol/primidone) and -buffer (H2O/DMSO; control) correlation heatmaps for NPCs. Value indicates correlation between gene-trait and gene-module associations with p value in parenthesis. c Module dendrograms with module membership correlation heatmaps for DAOYs. d Module dendrograms with module membership correlation heatmaps for NPCs.

Table 3 Pathway enrichment analysis of red gene module for drug treatment in DAOYs.Table 4 Pathway enrichment analysis of red gene module for drug treatment in NPCs.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 GWAS21 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).

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.

We postulated that since propranolol and primidone had a non-significant 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 pathophysiology22. 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 single-nucleus 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-, DYRK-expressing 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).

Fig. 4: Single-cell enrichment of ET genes in cerebellar and cortical tissues.

a Single-cell enrichment Z-score heatmap of ET-related genes in adult cerebellar tissue. Rows represent ET genes; columns represent cell types of the cerebellum (Purk1 SORC3+ Purkinje cells, Purk2 SORC3− Purkinje cells, Ast Astrocytes, OPC Oligodendrocyte progenitor cells, Oli Oligodendrocytes, Mic Microglia, End Endocytes, Gran Granule cells, Per Pericytes). b Ridge plots displaying distribution of average expression counts of ET-related genes in different cell types of the adult cerebellum. c Z-score expression heatmap of ET genes in single-cell types of the adult cortex. Rows represent ET genes; columns represent cortical cell types (Exc Excitatory, Inh Inhibitory, Astro Astrocytes, End Endocytes, Peri Pericytes, VLMC vascular and leptomeningeal cells, OPC Oligodencrocyte progenitor cells, Oligo Oligodendrocytes, L# cortical layer).

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.

Next, we assessed the enrichment of propranolol and primidone DEGs identified in this study in cortical and cerebellar single-cell data using a one sample Z-test (Eq. (4); Fig. 5 and Supplementary Tables 17 and 18). In cerebellum single-nucleus data, convergent propranolol DEGs were mostly enriched in endocytes (z-score = 3.38, q val = 0.014) and microglia (z-score = 3.36, q val = 0.014) whilst convergent propranolol/primidone DEGs in all cell types were mostly enriched in oligodendrocytes (z-score = 2.90, q val = 0.034; Fig. 5e). Interestingly, convergent propranolol/primidone DEGs in DAOYs, a cell type specific to the cerebellum, had enriched expression in astrocytes (z-score = 2.74, q val = 0.047), much like the enrichment of ET genes in cerebellar astrocytes (Fig. 4a). In cortical tissue, convergent drug DEGs were mostly significantly enriched in non-neuronal cell types (Fig. 5d), notably oligodendrocytes (z-score = 5.09, q val = 3.65E−07), astrocytes (z-score = 4.92, q val = 1.00E−04) and endocytes (z-score = 3.95, q val = 1.70E−03). Unsurprisingly, given the use of propranolol to lower blood pressure, convergent propranolol DEGs were mostly enriched in endocytes (z-score = 6.18, q val = 4.48−07) and vascular and leptomeningeal cells (z-score = 4.77, q val = 1.52E−04). Of note, propranolol DEGs were also enriched in L1-L3 inhibitory neurons, notably vasoactive intestinal peptide (VIP) expressing inhibitory neurons (Fig. 5d and see Supplementary Tables 17 and 18 for statistics).