Transcriptomic Changes Resulting From STK32B Overexpression Identify Pathways Potentially Relevant to Essential Tremor

Objective: Essential tremor (ET) is a common movement disorder that has a high heritability. A number of genetic studies have associated different genes and loci with ET, but few have investigated the biology of any of these genes. STK32B was significantly associated with ET in a large genome-wide association study (GWAS) and was found to be overexpressed in ET cerebellar tissue. The objective of this study is to determine the effects of overexpressed STK32B in cerebellar DAOY cells. Methods: Here, we overexpressed STK32B RNA in human cerebellar DAOY cells and used an RNA-Seq approach to identify differentially expressed genes (DEGs) by comparing the transcriptome profile of these cells to one of the control DAOY cells. Results: Pathway and gene ontology enrichment identified axon guidance, olfactory signaling, and calcium-voltage channels as significant. Additionally, we show that overexpressing STK32B affects transcript levels of previously implicated ET genes such as FUS. Conclusion: Our results investigate the effects of overexpressed STK32B and suggest that it may be involved in relevant ET pathways and genes.


INTRODUCTION
Essential tremor (ET) is one of the most common movement disorders and is typically characterized by a kinetic tremor in the hand or arms (Elble, 2013). The severity tends to increase with age and may involve different regions such as the head, voice, or jaw (Haubenberger and Hallett, 2018). Previous studies investigating the histology of post-mortem tissue have pinpointed the cerebellum as a region of interest for ET (Gironell, 2014). Specifically, abnormalities with Purkinje cell axons and dendrites were found in ET brains (Axelrad et al., 2008). Furthermore, the olivocerebellar circuitry has been implicated in ET pathology and many calcium voltage channels are highly expressed in this circuitry (Schmouth et al., 2014).
The genetic etiology of ET has remained largely elusive and most studies have focused on common or rare variants and on looking for genetic overlap with other disorders (Houle et al., 2017;Liao et al., 2018). Twin studies have shown that ET has a concordance of 69-93% in monozygotic twins and 27-29% in dizygotic twins, which suggests that both genetic and environmental factors drive the onset and development of this complex trait (Tanner et al., 2001). A recent genome-wide association study (GWAS) identified a significant locus in STK32B and found that ET patients overexpressed STK32B in cerebellar tissue by comparison to healthy controls, suggesting potential implications of STK32B in ET (Müller et al., 2016). STK32B is transcribed and translated into YANK2, a serine/threonine kinase, that has not been well characterized. There have been several exome-wide studies that implicated different genetic variants as causes of ET. The first ET-implicated gene found through exome sequencing was the Fused in Sarcoma gene (FUS) (Merner et al., 2012). However, it is unclear whether these "ET" genes interact or have indirect effects on the expression of each other.
To understand the effects of overexpressed STK32B mRNA and identify pathways potentially relevant to ET, we overexpressed this gene in human cerebellar DAOY cells and compared transcriptomic changes in overexpressed cells and empty-vector controls using RNA sequencing. Several interesting pathways such as axon guidance, calcium ion transmembrane transport, and olfactory transduction were significantly enriched after overexpression of STK32B. We also identified previously implicated ET genes whose expression is dysregulated through the overexpression of STK32B, suggesting that overexpressed STK32B may have relevant downstream effects.

Confirming Overexpression of STK32B
The STK32B stable cell lines had higher RNA expression of STK32B, which was detected by reverse-transcriptase qPCR (Supplementary Figure S1). RNA sequencing (RNA-Seq) also confirmed that the stable cell lines have higher STK32B RNA levels compared to controls (Figures 1, 2). The STK32B gene was the top significant differentially expressed gene (P = 7.07E-245, β = 3.1789).

Quality Control of RNA-Seq Analyses
The QQ-plot for differentially expressed genes (DEGs) did not show stratification or potential biases (Supplementary Figure S2). The M-A plot showed slight enrichment of downregulated genes compared to upregulated genes (Supplementary Figure S3). A principal component (PC) plot showed separation of controls and STK32B cell lines (Figure 3). Additionally, the dendrogram with heatmap shows that controls and STK32B cell lines are distinct (Supplementary Figure S4).
The mean-variance plot shows the shrinkage under the sleuth model (Supplementary Figure S5).

Gene Network Analyses
To identify which cluster of DEGs drove the significant pathways, gene network analysis based on co-expression was done using and identified 14 different clusters within the DEGs (Supplementary Table S1). Analysis of the top 3 largest clusters for pathway enrichments found several similar pathways in different databases including axon guidance (P < 1.09E-12) in cluster 1 and calcium signaling pathway in cluster 3 (2.06E-05).

DISCUSSION
The genetic etiology of ET is complex and likely explained by a combination of copy number variants, rare variants, gene-gene interactions and common variant drivers. One of the most significantly enriched pathways we identified was olfactory signaling and transduction. Previous reports have shown conflicting evidence about olfactory loss in ET (Applegate and Louis, 2005;Shah et al., 2008). It may be that ET patients with dysregulated olfactory signaling have overexpressed STK32B, which may contribute to the subset of ET patients with olfactory loss.
Interestingly, FUS was amongst the genes dysregulated when STK32B was overexpressed. In the exome familial ET study that identified FUS, we also found reduced mRNA levels for FUS (Merner et al., 2012). Similarly, the expression of FUS is lower when STK32B is overexpressed, suggesting an indirect relationship between the two genes. Additionally, two calcium voltage-gated channel genes, CACNA1C and CACNA1A were overexpressed. Enrichment analyses found the GO term calcium ion transmembrane transport to be highly significant in GO Biological Processes, suggesting that STK32B may play a role upstream of these genes. In the olivocerebellar circuitry, a system implicated in ET, is enriched for both of these genes (Mori et al., 1991;Schmouth et al., 2014). CACNA1A has been shown to be predominantly expressed in Purkinje cells, a cell type relevant to ET (Mori et al., 1991). Another enriched pathway was axon guidance. A previous study identified the TENM4 missense variant segregating within ET families (Hor et al., 2015). TENM4 is a regulator of axon guidance and myelination and Tenm4 knockout mice showed an ET phenotype, suggesting that axon guidance is important in ET (Hor et al., 2015). Similarly, overexpressed STK32B could be dysregulating important pathways relevant to axon guidance.
By sub-stratifying the overexpressed genes into different clusters that are co-expressed, several interesting pathways involving the cardiovascular system were identified. GO terms and pathways such as right ventricular cardiomyopathy in KEGG and angiogenesis in GO Biological Processes were found among the significant pathways in Table 3. In certain ET patients, beta-blockers can reduce tremor magnitude and frequency. Beta-blockers lower blood pressure and are used to treat irregular heart rhythm and other cardiovascular phenotypes. A common treatment for ET is propranolol, a beta-blocker. This was a drug developed for cardiovascular health that affects adrenergic activity, but still reduces tremor in ET individuals. This could suggest that overexpressed STK32B affects cardiovascular health, which may in turn affect or lead to the ET phenotype or that STK32B may be pleiotropic and affect the nervous and cardiovascular system differently.
Due to the heterogeneity of ET, it is likely that not all ET-affected individuals have overexpression of STK32B. One limitation of this study is that DAOY cells are a cancerous cell line and may not be the ideal model for understanding the transcriptome of Purkinje cells. Future studies on the effects of STK32B overexpression could use induced pluripotent stem cells differentiated into cells with Purkinje cells properties; such cells would likely have less biological noise. However, our findings support the idea that STK32B is a gene of interest for a subset of ET cases, and further investigations into how STK32B may interact with other genes are warranted.

Cell Culture and Plasmid Construction
The DAOY cell line was cultured in Eagle's Minimum Essential Medium (EMEM) with 10% fetal bovine serum (FBS) at 37 • C and 5% CO 2 with the Glico Pen-Strep-Glutamine cocktail at 1×. Cells were passaged every 2 days at 80-90% confluence and at the same time. The cDNA of STK32B (NM_001306082.1) was inserted into a pcDNA3.1(+) vector containing a neomycin resistance gene, used as a selectable marker. This transcript was picked because it had the highest expression in the cerebellum based on GTEx 53 v7. The plasmid was transformed and expanded in XL10-Gold Escherichia coli strain (Stratagene). The cDNA of the plasmid was sent for Sanger sequencing to validate the gene sequence.

Transfection and Deriving Stable Cell Lines
The cells were transfected for 48 h with 8 ug of the vector of using the jetPRIME transfection reagent (Polyplus). A green fluorescent protein (GFP) control vector was used to assess and estimate transfection efficiency and optimize transfection parameters. Stable cell lines were established by subjecting transfected cells to G418 antibiotic at 1 µg/mL for 10 days and maintained at a concentration of 400 µg/mL for an additional 2 weeks. A kill curve was done to determine the ideal antibiotic selection. In parallel, empty pcDNA3.1(+) vector control lines were also grown and underwent antibiotic selection.

RNA Extraction and Sequencing
Total RNA was extracted using the RNeasy Mini Kit (Qiagen). The total RNA concentration was measured using the Synergy H4 microplate reader. RNA for the top three overexpressed STK32B cell lines and three empty-vector controls was sent to Macrogen Inc., for sequencing. RT-qPCR was used to confirm overexpression in the cell lines compared to controls. Preparation of the cDNA library was done with the TruSeq Stranded Total RNA Kit (Ilumina) with Ribo-Zero depletion. Sequencing was done on the NovaSeq 6000 at 150 bp paired-end reads with a total of 200M reads. Samples were randomized for cell lysis, RNA extraction, Ribo-Zero depletion, library preparation and sequencing to account for potential batch effects.

Data Processing and Quality Control
The FASTQ files were pseudo-aligned with Salmon using the Ensembl v94 annotation of the human genome (Patro et al., 2017). The data is available on the public repository of NCBI, Gene Expression Omnibus (GEO), with the GEO accession as GSE150393. It can be accessed with the following URL: https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150393. The parameters of Salmon included the following: 200 bootstraps with mapping validation, GC bias correction, sequencing bias correction, 4 range factorization bins, a minimum score fraction of 0.95 and VB optimization with a VB prior of 1e-5. The estimated counts were analyzed using sleuth to identify DEGs (Pimentel et al., 2017). A Wald test was used to get β values [log 2 (x + 0.5)]. Additionally, transcript-level aggregation was done against the Ensembl v94 annotation to determine gene-level differential expression. P-values were corrected using the Benjamini-Hochberg procedure to account for false discovery rate (FDR). The FDR threshold was set at 0.05. A PC plot and heatmap was done to identify clusters.

Pathway and Gene Set Enrichment Analysis
Pathway and gene set enrichment for relevant databases was done using Gene Network 2.0 (Deelen et al., 2019). To investigate genes with larger fold changes, only DEGs with a β > | 0.5| were analyzed for pathway enrichments. Network clustering was done by grouping genes based on co-expression using Gene Network 2.0. Pathway analyses were also conducted with Gene Network 2.0 for the clusters using a Wilcoxon test. We included the top five significant pathways across different databases.

DATA AVAILABILITY STATEMENT
The data is available on the public repository of NCBI, Gene Expression Omnibus (GEO), with the GEO accession as GSE150393. It can be accessed with the following URL: https: //www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150393.

ETHICS STATEMENT
The review board at the McGill University Health Center Research Ethics Board (MUHC REB) approved the study protocols (reference number IRB00010120).

AUTHOR CONTRIBUTIONS
CL conducted the experiments and analyses and drafted the manuscript. FS, VV, DR, FA, SD, GH, QH, and HC helped with experiments. AL and DS helped with bioinformatic analyses. PD and GR oversaw the study and drafting of the manuscript. All authors contributed to the article and approved the submitted version.