Whole‐exome sequencing for variant discovery in blepharospasm

Abstract Background Blepharospasm (BSP) is a type of focal dystonia characterized by involuntary orbicularis oculi spasms that are usually bilateral, synchronous, and symmetrical. Despite strong evidence for genetic contributions to BSP, progress in the field has been constrained by small cohorts, incomplete penetrance, and late age of onset. Although several genetic etiologies for dystonia have been identified through whole‐exome sequencing (WES), none of these are characteristically associated with BSP as a singular or predominant manifestation. Methods We performed WES on 31 subjects from 21 independent pedigrees with BSP. The strongest candidate sequence variants derived from in silico analyses were confirmed with bidirectional Sanger sequencing and subjected to cosegregation analysis. Results Cosegregating deleterious variants (GRCH37/hg19) in CACNA1A (NM_001127222.1: c.7261_7262delinsGT, p.Pro2421Val), REEP4 (NM_025232.3: c.109C>T, p.Arg37Trp), TOR2A (NM_130459.3: c.568C>T, p.Arg190Cys), and ATP2A3 (NM_005173.3: c.1966C>T, p.Arg656Cys) were identified in four independent multigenerational pedigrees. Deleterious variants in HS1BP3 (NM_022460.3: c.94C>A, p.Gly32Cys) and GNA14 (NM_004297.3: c.989_990del, p.Thr330ArgfsTer67) were identified in a father and son with segmental cranio‐cervical dystonia first manifest as BSP. Deleterious variants in DNAH17,TRPV4,CAPN11,VPS13C,UNC13B,SPTBN4,MYOD1, and MRPL15 were found in two or more independent pedigrees. To our knowledge, none of these genes have previously been associated with isolated BSP, although other CACNA1A mutations have been associated with both positive and negative motor disorders including ataxia, episodic ataxia, hemiplegic migraine, and dystonia. Conclusions Our WES datasets provide a platform for future studies of BSP genetics which will demand careful consideration of incomplete penetrance, pleiotropy, population stratification, and oligogenic inheritance patterns.


| INTRODUCTION
Dystonia is defined as a movement disorder characterized by sustained or intermittent muscle contractions causing abnormal, often repetitive, movements, postures, or both (Albanese et al., 2013). In general, adult-or late-onset dystonia without evidence of overt degeneration or structural lesions of the nervous system is referred to as isolated dystonia and can be inherited in an autosomal-dominant fashion with reduced penetrance. The most common forms of focal dystonia are cervical dystonia and blepharospasm (BSP). Blepharospasm (BSP) (OMIM: 606798) is characterized by involuntary orbicularis oculi spasms that are usually bilateral, synchronous, and symmetrical (Defazio et al., 2015). Review of BSP epidemiological data provides prevalence estimates ranging from 16 to 133 per million (Defazio, Abbruzzese, Livrea, & Berardelli, 2004). BSP is significantly more common in females (>2F:1M) with a mean age of onset at approximately 55 years (O'Riordan et al., 2004). In comparison to cervical and laryngeal dystonia, BSP is more likely to spread to other body parts (Weiss et al., 2006). Most commonly, BSP spreads to contiguous craniocervical segments (lower face, masticatory muscles, and neck). The term segmental craniocervical dystonia is applied to the combination of BSP and dystonia of other head and neck muscles (LeDoux, 2009). Herein, BSP-plus (BSP+) will be used to denote subjects with BSP who exhibit subsequent spread to other anatomical segments (LeDoux, 2009;Waln & LeDoux, 2011). Sensory tricks or geste antagonistes are highly specific to dystonia, reported in a high percentage of patients with BSP, and can facilitate the diagnosis of BSP (Defazio, Hallett, Jinnah, & Berardelli, 2013). However, without valid genetic biomarkers, the diagnosis of BSP can be difficult, even for experienced clinicians (Defazio et al., 2013).

| Ethical compliance
All human studies were conducted in accordance with the Declaration of Helsinki with formal approval from the University of Tennessee Health Science Center Institutional Review Board (IRB; 01-07346-FB, 05-08331-XP, and 14-03320-XP) and ethics committees of all participating centers. All subjects gave written informed consent for genetic analyses and disclosure of medical information.

| Subjects
Subjects in this study were examined by at least one neurologist with subspecialty expertise in movement disorders. Subjects were asked to perform specific tasks, including holding their eyes open, opening and closing their eyes gently, opening and closing their eyes forcefully, along with additional verbal and postural maneuvers designed to capture masticatory, laryngeal or cervical involvement. A clinical diagnosis of definite BSP was given to subjects that exhibited increased blinking and stereotyped, bilateral and synchronous orbicularis oculi spasms inducing narrowing/closure of the eyelids (Defazio et al., 2013). Subjects with isolated episodes of increased eyelid blinking were given a diagnosis of possible BSP. Each affected or possibly affected family member was queried for the presence of sensory tricks. WES was completed on a total of 31 subjects from 21 pedigrees from the United States, Poland, and Italy (Table 1). Prior to WES, pathogenic variants in THAP1, GNAL (OMIM 139312) and Exon 5 of TOR1A were excluded as previously described Vemula et al., 2013;Xiao et al., 2009Xiao et al., , 2010. Two pedigrees were African-American and 19 pedigrees were Caucasian of European descent. The results of WES on the proband of African-American pedigree 10908 were previously reported (Xiao, Thompson, Vemula, & LeDoux, 2016) and deposited in Sequence Read Archive (SRX1790848).

| Whole-exome sequencing
The concentration and quality of genomic DNA (gDNA) extracted from peripheral blood were examined with a NanoDrop â ND-1000 (Thermo Scientific), the Qubit â dsDNA BR Assay Kit (Thermo Scientific) and agarose gel electrophoresis. DNA was then forwarded to Otogenetics or Beijing Genomics Institute (BGI) for additional in-house quality control assessments prior to WES.
For WES at Otogenetics, 3 lg of genomic DNA (gDNA) was sheared to yield 100-450 bp fragments. Insolution whole-exome capture and massively parallel sequencing was performed using the Agilent SureSelect XT All Exon Kit 51 Mb. Enriched DNA fragments were sequenced on Illumina's HiSeq 2500 platform as pairedend 100-125 base-pair reads. On average, over 95% of exons were covered at >209. The percentage of exome coverage was based on exons targeted by the 51 Mb All Exon v4 Kit which incorporates Consensus Coding Sequence (CCDS), NCBI Reference Sequence (RefSeq) and GENCODE annotations.
For WES at BGI, the gDNA samples were fragmented by Covaris, and, after two rounds of bead purification, the resulting gDNA fragments were mainly distributed between 200 and 400 bp. Then, AdA 5 0 -and 3 0 -adaptors were ligated to the 5 0 -and 3 0 -ends of the fragments, respectively. The AdA adaptor-ligated fragments were amplified by PCR, and the PCR products were used for exon capture. A 58.95 Mb region was targeted for capture. The captured exon fragments were purified by DynabeadsM-280 streptavidin bead purification and were further amplified by another round of PCR. Then, the PCR products were circularized and the resulting double strand (ds) circles digested with Ecop15. Among these digested fragments, small fragments were collected after bead purification. Similar to the AdA adaptor ligation, AdB adapters were ligated to both ends of the purified fragments and the fragments were then used for single strand (ss) circularization. The resulting ss circles were the final library products used on the CG Black Bird sequencing platform. Finally, highthroughput sequencing was performed for each captured library.

| Read mapping
Sequence reads (FASTQ) from Illumina (Otogenetics) were mapped to the human reference genome (NCBI build 37.1) with NextGENe â (SoftGenetics). Using the consolidation and elongation functions of NextGENe, instrument sequencing errors were reduced and sequence reads were lengthened prior to variant analysis. The condensation tool polished the data for adequate coverage by clustering similar reads with a unique anchor sequence. Using this process, short reads were lengthened and reads with errors were filtered or corrected. To maximize the probability of detecting causal variants, all base changes occurring in ≥4 reads in any individual sample were classified as variants for downstream analyses. An Overall Mutation score of 5 was used as a cutoff to filter read errors and reduce the effects of allelic imbalances. The Overall Mutation score is generated via a proprietary algorithm (SoftGenetics) to provide an empirical estimation of the likelihood that a given variant call is genuine and not an artifact of sequencing or alignment errors. This score is based on the concept of Phred scores, where quality scores are logarithmically linked to error probabilities. With NextGENe â software, intergenic and deep intronic (≥12 nt from splice sites) variants were eliminated prior to downstream in silico analyses. Complete Genomics (BGI) developed high-speed mapping software capable of aligning read data to reference sequences. Using GRCh37 as the reference, the mapping is tolerant of small variations from a reference sequence, such as those caused by individual genomic variation, read errors, or unread bases. To support assembly of larger variations, including large-scale structural changes or regions of dense variation, each arm of a DNA Nanoball (DNB) is mapped separately, with mate pairing constraints applied after alignment. Initially, mapping reads to the human reference genome is a constrained process that does not allow for insertions and deletions. All mate-pair constraint-satisfying paired-end mappings are used to detect small variants. DNBs are then filtered and individual reads are optimized. Optimization collects reads likely to lie in regions of interest, using mate alignment information and performs local de novo assemblies.
REVEL, MetaLR and CADD scores were used to prioritize nonsynonymous missense variants for additional scrutiny whereas CADD and ExAC Probability of Loss-of-Function (LoF) intolerance (pLI) scores were used to prioritize nonsense SNVs and frameshift INDELs. Muta-tionTaster was also used for analysis of small INDELs which are not scored by REVEL or MetaLR. Each category of variant (nonsynonymous, synonymous, splice-site, nonsense, frameshift, other INDELs, and other SNVs) was ranked by in silico scores of deleteriousness. Population frequencies for the highest scoring variants were additionally assessed with genome Aggregation Database (gno-mAD), NHLBI Exome Sequencing Project (ESP) Exome Variant Server (EVS) with particular attention to racial subcategories. All NCBI databases were queried with gene symbols and the names of encoded proteins. Particular attention was paid to data contained in PubMed, ClinVar, OMIM, and BioSystems. OMIM was searched for allelic disorders/phenotypes. MARRVEL and its link outs were used to explore available data related to animal models of homologs, genomic structural variants (DGV and DECI-PHER), gene expression (GTex), and protein expression (ProteinAtlas). Gene expression was also analyzed with Allen Brain Atlas and BioGPS. Candidate genes were eliminated if not expressed in at least one "motor" region of the brain (striatum, cerebellum or frontal motor cortex). UniProt was used to access protein-protein interactions, sites of known or predicted posttranslational modifications and known or putative protein functions. Multiple sequence alignments were performed with Clustal Omega. A subset of candidate pathogenic variants was confirmed with bidirectional Sanger sequencing to exclude next generation sequencing read errors. After Sanger confirmation, cosegregation was assessed in individual pedigrees.

| Copy number variant analysis
CNVkit (Talevich, Shain, Botton, & Bastian, 2016), a Python library and command-line software toolkit to infer and visualize copy number variants (CNVs) from targeted DNA sequencing data, was used to detect CNVs in WES data generated by Otogenetics on the Illumina platform. CNVkit was designed for use on hybrid capture sequencing data where off-target reads are present and can be used to improve copy number estimates. CNVkit normalizes read counts to a pooled reference and corrects for three main sources of bias: GC content, target footprint size, and repetitive sequences. For this purpose, Otogenetics provided us with WES data from 15 random subjects of unknown race and unknown geographic region of origin sequenced as part of unrelated projects using the Agilent SureSelect XT All Exon Kit 51 Mb for exome capture and sequenced on Illumina's HiSeq 2500 platform.
CNVkit reports log2 copy ratios. Assuming pure samples and germline mutations, the log2 ratio should be À1.0 for a deletion mutation and infinity if both alleles are deleted. The log2 ratio is 0.585 for duplications and 1.0 for triplications. The relationship between the estimated copy number and the true copy number depends on a number of factors including read depth and number of probes covering a region of interest.

| Sanger sequencing
PCR was performed using 40 ng of peripheral blood gDNA along with 200 nmol/L of each primer (Table S1) in a 10-ll reaction volume with HotStarTaq â Plus DNA polymerase from Qiagen. The following cycling conditions were employed: 95°C for 15 min; 35 cycles at 95°C for 10 s, 58°C for 30 s, and 72°C for 30 s.

| PCR validation of copy number variants
Quantitative PCR (qPCR) was used for initial assessment of a random selection of predicted CNVs identified with CNVkit. Primers and probes for qPCR were designed with Roche's Universal Probe Library to cover (Table S1). qPCR was performed using 20 ng of template DNA and 200 nmol/L of each primer in a 10-ll reaction volume with the LightCycler TM 480 system and Universal Taqman â probes (Roche). The following cycling conditions were employed: 95°C for 5 min; 45 cycles at 95°C for 10 s, 58°C for 30 s, and 72°C for 12 s. Copy numbers were calculated against an endogenous control, HLCS, holocarboxylase synthetase. All assays were carried out in triplicate and means were used for calculating fold changes.
Digital PCR (dPCR) was then used for confirmation of select deletion and duplication CNVs identified with CNVkit. Literature mining as described for SNVs and small lNDELs was used to select genes with deletion log2 scores of À0.75 to À1.25 and covered by ≥4 probes, or genes with duplication log2 scores of 0.385 to 0.835 and covered by ≥4 probes. Primers and probes (FAM dyelabeled) were designed via Roche's Universal Probe Library to encompass the estimated deletion regions (Table S1). The TaqMan copy number reference assay (Applied Biosystems 4403326) contained RNase P-specific F I G U R E 1 Blepharospasm (BSP) and BSP+ Pedigrees. Pedigrees with two or more affected individuals. Arrows, probands. Arrowheads, other family members analyzed with whole-exome sequencing. White symbol, unaffected. Black symbols, BSP, BSP+ or other anatomical distribution of dystonia. Gray symbols, possibly affected forward and reverse primers and VIC dye-labeled TAMRA hydrolysis probe. RNase P, a single copy gene, is used as the reference for this work (Qin, Jones, & Ramakrishnan, 2008).
Reaction mixtures (4.0 ll) containing TaqMan geneexpression master mix (Life Technologies), 20X GE sample loading reagent (Fluidigm 85000746), 20X gene-specific assays, 20X TaqMan copy number reference assay (Applied Biosystems) and 1.2 ll target gDNA (20 ng/ll) was pipetted into each loading inlet of a 48.770 dPCR array (Fluidigm). The BioMark IFC controller MX (Fluidigm) was used to uniformly partition the reaction from the loading inlet into the 770 9 0.84 nl chambers and dPCR was performed with the Fluidigm BioMark System for Genetic Analysis. The Fluidigm dPCR software was used to count gene copy numbers. The quality thresholds were manually set specific to each assay, but consistent across all panels of the same assay. The CNV calculation is based on "relative copy number" so that apparent differences in gene copy numbers in different samples are not distorted by differences in sample amounts. The relative copy number of a gene (per genome) is expressed as the ratio of the copy number of a target gene to the copy number of a single copy reference gene in the sample. By using assays for the two genes (the gene of interest and the reference gene) with two fluorescent dyes on the same Digital Array IFC, we are able to simultaneously quantitate both genes in the same DNA sample. The ratio of these two genes is the relative copy number of the gene of interest.

| BSP and BSP+ pedigrees
Whole-exome sequencing was completed on 31 subjects from 21 distinct pedigrees with either concordant or discordant BSP and BSP+ phenotypes ( population (1.49E-02) with a much lower allele frequency of (6.76E-04) in non-Finnish Europeans and quite rare in other racial populations. The identified amino acid substitution is located in the C-terminal, intracellular domain of the encoded voltage-dependent P/Q-type calcium channel subunit a-1A, which is conserved among mammals (Figure 3). We did not screen other variants for cosegregation given previously established associations between CAC-NA1A and dystonia. Five SNVs had CADD_phred scores >15 and REVEL scores >0.5 but none had a MetaLR score >0.75, REVEL score >0.75 and CADD_phred score >30. A frameshift INDEL in MMP28 with a CADD_phred score of 34 is reported in ExAC and gnomAD. Four nonsense SNVs had CADD_phred scores >30 but two are reported in ExAC and gnomAD and none seem biologically plausible candidates.

| REEP4 missense variant
was identified in seven subjects with BSP+ or BSP and one asymptomatic female family member from a three-generation African-American pedigree ( Figure 4, Tables 1, 5, 8 and S2; Data S1). This variant is present at very low frequency in gno-mAD and predicted to be deleterious by in silico analysis, including CADD (phred score = 34), REVEL (0.767), MetaLR (0.960), and MutationTaster2 (disease causing, probability value: 1.0). In gnomAD, this variant is not present in 15,290 African alleles. The p.Arg37Trp variant alters an amino acid that is highly conserved among vertebrates as shown by the multiple pairwise alignments generated with Clustal Omega (Figure 4).

| DNAH17 variants found in pedigree and isolated subject with BSP
Deleterious variants in DNAH17 were identified in two brothers with BSP and one isolated case of BSP (Figure 7; Tables 1, 6, 8 and S2, Data S1). Both variants are present at low frequency in ExAC and gnomAD. DNAH17 encodes dynein axonemal heavy chain 17. The FANTOM5 dataset reports expression of DNAH17 in testes and brain (hippocampus, caudate and cerebellum; Kawaji, Kasukawa, Forrest, Carninci, & Hayashizaki, 2017). DNAH17 has not yet been linked to any other neurological or non-neurological disease. A roundworm homolog (dhc-1) of human DNAH17 is involved in cytokinesis, microtubule-based movement, mitotic spindle organization, meiotic nuclear division and nervous system development (MARRVEL).

| Copy number variants
CNVkit called from 11 to 217 CNVs per shared exome. Assessing randomly selected CNVs with qPCR showed high discordancy (Table S3), particularly for variants that did not have log2 ratios near À1.0. We then focused on CNVs with log2 ratios compatible with a single-copy gain (~0.585) or single-copy loss (À1.0) using dPCR. Deletions in LILRA3 were confirmed in three unrelated subjects with BSP (Table 7). LILRA3 (OMIM 604818) deletions are common in the general population and may increase risk for HIV infection and some autoimmune disorders (Ahrenstorf et al., 2017;Du et al., 2015). A deletion in BTNL3 (OMIM 606192) and duplications in SLC2A14 (OMIM 611039), SLC2A3 (OMIM 138170), TOP3B (OMIM 603582), and UNK (616375) were identified in single exomes (Tables 7 and 8). UNK is expressed at high levels in brain (Allen Brain Atlas, BioGPS, and The Human Protein Atlas) and plays an important role in the development of neuronal morphology. Two UNK duplications are reported in ExAC. To date, UNK has not been linked to any medical disorder (OMIM). Copy number analysis of GOLGA8A (Chr15) was compromised by the presence of pseudogenes and a homolog with very close sequence similarity on Chr15.
3.9 | Other candidate genes found in two or more pedigrees  (Granneman et al., 2003). A missense variant in UBR4 (p.Arg5091His) was found to segregate with episodic ataxia in a large Irish pedigree (Conroy et al., 2014). UBR4 is expressed at high levels in cerebellar Purkinje cells (Allen Brain Atlas), interacts with calmodulin, colocalizes with ITPR1, and may be involved in Purkinje cell calcium homeostasis (Conroy et al., 2014). ARHGEF19 shows significant expression in cerebellar Purkinje cells (Allen Brain Atlas) and zebrafish arhgef19 is involved in neural tube closure (Miles et al., 2017).

| DISCUSSION
The molecular and cellular mechanisms underlying BSP and other anatomical distributions of isolated dystonia remain fragmentary. Accordingly, treatments for BSP are entirely symptomatic (Pirio Richardson et al., 2017). Most commonly, BSP patients are treated with injections of botulinum toxin although, in some series, almost 50% report minimal improvement, no improvement or worsening of BSP after injections of botulinum toxins (Fernandez et al., 2014). Identification of genetic etiologies for BSP may permit development of targeted disease-modifying therapeutics. In this study, we used exome sequencing to explore genetic contributions to BSP and provide a foundation for future case-control studies of this important focal dystonia. Although we do provide data suggesting potential roles for CACNA1A, REEP4, TOR2A, ATP2A3, HS1BP3/GNA14, DNAH17, TRPV4, CAPN11, VPS13C, UNC13B, SPTBN4, MYOD1, and MRPL15 in the pathogenesis of BSP, the limitations of our work should be bordered. First, we did not identify a common cosegregating genetic etiology in more than one pedigree. This points to the likely genetic heterogeneity of BSP but also suggests that one or more variants identified herein cosegregated by chance alone. Unfortunately, none of our pedigrees were powered to generate LOD (logarithm [base 10] of odds) scores >3 thereby precluding the usage of linkage analysis for validation of  cosegregating variants. Second, several of the candidate variants identified with WES are reported in population databases (ExAC and gnomAD) with MAFs near the minimal population prevalence of BSP. On the other hand, noted MAFs are significantly lower than the maximal population prevalence of BSP with corrections for the markedly reduced penetrance characteristic of isolated dystonia. Furthermore, BSP and premonitory increased blinking may be much more common in the general population than commonly accepted (Conte et al., 2017). Thirdly, our genetically heterogeneous cohort included Polish, Italian, Caucasian-American and African-American pedigrees, possibly reducing the probability of detecting variants shared among pedigrees and singletons. Accordingly, follow-up case-control analysis of individual variants identified herein will require careful attention to population stratification and large sample sizes to confidently determine if variants in candidate genes are enriched in BSP. Fourth, our prioritization of variants was predominantly driven by in silico predictions of deleteriousness and many potentially pathogenic candidate variants were not confirmed with Sanger sequencing or subjected to cosegregation analysis. Fifth, WES will miss most repeat expansions and does not access the mitochondrial genome. In this regard, repeat expansions are a common cause of late-onset neurological disease and mitochondrial mutations may include dystonia as part of a more expansive neurological phenotype (LeDoux, 2012). Furthermore, our approach to CNV analysis was largely insensate to smaller structural variants such as single exonic deletions. Despite these limitations, our findings are compatible with common themes in dystonia research (calcium signaling, Purkinje cells, and dopaminergic signaling), point out potential genetic common ground with PD and ET, suggest a role for oligogenic inheritance in BSP, and provide motivation for treating a subset of BSP patients with acetazolamide.
CACNA1A is highly expressed in the cerebellum, particularly the Purkinje cell layer. Mutations in several genes related to calcium signaling and homeostasis and expressed in Purkinje cells have been causally associated with dystonia in humans and mice (LeDoux, 2011). In fact, virtually all genes associated with dystonia in spontaneous mutants (tottering, stargazer, ophisthotonus, ducky, lethargic, waddles, and wriggle) are involved in Purkinje cell Ca 2+ signaling (Canca1a,Cacng2,Itpr1,Cacna2d2,Cacnb4,and Pmca2). In humans, autosomal-recessive mutations in HPCA (OMIM 142622) cause childhood-onset dystonia and the encoded protein, hippocalcin, is robustly expressed in Purkinje cells and serves as a Ca 2+ sensor (Charlesworth et al., 2015;Tzingounis, Kobayashi, Takamatsu, & Nicoll, 2007). SVs in CACNA1A have been associated with a variety of neurological disorders including episodic ataxia type 2, familial hemiplegic migraine, spinocerebellar ataxia type 6 (SCA6), and various anatomical distributions of dystonia such as benign paroxysmal torticollis of infancy and BSP (Naik, Pohl, Malik, Siddiqui, & Josifova, 2011;Sethi & Jankovic, 2002;Shin, Douglass, Milunsky, & Rosman, 2016 , 2008). A notable percentage of patients with dystonia due to mutations in CACNA1A show significant improvement with acetazolamide (Spacey, 1993;Spacey et al., 2005). Unfortunately, our pedigree was lost to follow-up and none of the affected family members were treated with acetazolamide. The a-1 subunit of P/Q type, voltage-dependent, calcium channel harbors a polyglutamine expansion in its C-terminal intracellular domain and the novel missense variant p.Pro2421Val identified in our pedigree with BSP is near this expansion (Figure 3). In contrast, the previously described BSP-variant was likely associated with nonsensemediated decay and haploinsufficiency (Spacey et al., 2005). Mutations linked to familial hemiplegic migraine appear to operate via gain-of-function mechanisms whereas the SCA6 polyglutamine repeat and loss-of-function mutations may lead to neuronal cell death (Cain & Snutch, 2011). In this context, it is worthy to note that reduced Purkinje cell density was found in two individuals with BSP and cervical dystonia (Prudente et al., 2013). REEP4 is a microtubule-binding endoplasmic reticulum and nuclear envelope protein (Schlaitz, Thompson, Wong, Yates, & Heald, 2013). Depletion of REEP4 from HeLa cells is associated with defective cell division and proliferation of intranuclear membranes derived from the nuclear envelope (Schlaitz et al., 2013). Similarly, omega-shaped nuclear blebs have been used as a phenotypic measure of torsinA (encoded by TOR1A) dysfunction (Laudermilch et al., 2016). In Xenopus, loss of REEP4 causes defects of nervous system development and paralysis of embryos (Argasinska et al., 2009). Mutations in REEP1 (OMIM 609139) and REEP2 (OMIM 609347) are associated with spastic paraplegia (SPG) types 31 (SPG31), and 72 (SPG72). Although dystonia is not a clinical feature typically reported in SPG31 and SPG72 cases, dystonia is not uncommon in several other SPGs, including SPG7, SPG15, SPG26, SPG35, and SPG47 (van Gassen et al., 2012;Klebe, Stevanin, & Depienne, 2015).
A DGAG deletion in Exon 5 of TOR1A was the first SV to be linked to isolated dystonia (Ozelius et al., 1997). TorsinA interacts with LAP1, a transmembrane protein ubiquitously expressed in the inner nuclear membrane. Recessive mutations of TOR1AIP1 (OMIM 614512) which encodes LAP1 are associated with severe early-onset generalized dystonia and progressive cerebellar atrophy (Dorboz et al., 2014). Another torsinA interacting protein, torsin family 2 member A (encoded by TOR2A) was found to harbor a missense variant in one of our pedigrees with BSP. Similar to the DGAG mutation in TOR1A, the penetrance of the p.Arg190Cys missense variant identified in our pedigree was less than 50%. TOR2A is a member of the human torsin gene family (Laudermilch et al., 2016;Ozelius et al., 1999 show relatively high expression in cerebellar Purkinje cells (Allen Brain Atlas). A nonsynonymous SNV in ATP2A3 (NM_005173.3: c.1966C>T, p.Arg656Cys) was found in five definitelyaffected subjects from a discordant pedigree with BSP from Italy. However, this variant was not detected in one possibly affected family member with writer's cramp. This could be either a phenocopy or evidence against the causality of ATP2A3. Furthermore, the p.Arg656Cys variant is present at notably high frequency in gnomAD (183/276,114 alleles, no homozygotes, 0.13% of 138,057 subjects). BSP is the most common focal dystonia in Italy with a crude prevalence rate of 133 per million or 0.013%. Even with a penetrance of <20%, this suggests that p.Arg656Cys may not be pathogenic or, at least, pathogenic in isolation, requiring digenic inheritance of another pathogenic variant. On the other hand, p.Arg656Cys is predicted to be highly deleterious, may contribute to other anatomical distributions of dystonia, and, like ATP1A3, could be involved in the etiopathogenesis of other neurological disorders such as Parkinson disease, Alzheimer disease, and brain tumors (Kawalia et al., 2017;Korosec, Glavac, Volavsek, & Ravnik-Glavac, 2009;Matak et al., 2016). In this regard, ATP2A3 shows striking expression in cerebellar Purkinje cells and dopaminergic neurons of the substantia nigra pars compacta (Allen Brain Atlas). ATP2A3 encodes a sarcoplasmic/endoplasmic reticulum Ca 2+ ATPase and disorders of Purkinje cell (LeDoux, 2011) and dopaminergic (Surmeier, Halliday, & Simuni, 2017) calcium homeostasis have been linked to dystonia and Parkinson disease, respectively.
A small pedigree ( Figure 6) with BSP+ and Parkinsonism harboring variants in HS1BP3 and GNA14 highlights the distinct possibility of oligogenic inheritance in BSP and other anatomical distributions of dystonia. In particular, all of the exomes sequenced in this study harbored more than one potentially pathogenic variant. Since most of our pedigrees were small and moderate numbers of variants showed in silico evidence of deleteriousness, we did not assess cosegregration for all of the identified candidate variants. However, we determined that both GNA14 and HS1PB3 were attractive candidate genes. Guanine nucleotide-binding protein subunit alpha-14 (encoded by GNA14) interacts with dynein, axonemal, light chain 4 (UniProt) which is expressed at high levels in sperm and brain. GNA14 appears to play a key role in the genetic architecture underlying normal gray matter density (Chen et al., 2015) and a GNA14 deletion mutation has been reported in a patient with early-onset Alzheimer disease (Lazarczyk et al., 2017). HS1BP3 shows moderate expression in brain (The Human Protein Atlas), and, in cerebellum, appears at highest levels in Purkinje cells (Allen Brain Atlas). Multipoint linkage analysis in four large pedigrees with ET identified a critical region between loci D2S2150 and D2S220 on Chr 2p which includes HS1BP3 (Higgins, Loveless, Jankovic, & Patel, 1998). The p.A265G HCLS1-binding protein 3 (HS1BP3) variant encoded by HS1BP3 is in linkage disequilibrium with ET but is unlikely to be causal since it is present at high frequency in the general population (Shatunov et al., 2005). It remains unknown if other coding or noncoding variants in HS1BP3 are causally related to the pathogenesis of ET. HS1BP3 negatively regulates autophagy (Holland et al., 2016), a cellular pathway closely tied to several neurodegenerative disorders including PD (Nash, Schmukler, Trudler, Pinkas-Kramarski, & Frenkel, 2017). In this regard, ET and PD may be related to adultonset dystonia through common genetics (De Rosa et al., 2016;Dubinsky, Gray, & Koller, 1993;Hedera et al., 2010;LeDoux et al., 2016;Louis et al., 2012;Straniero et al., 2017).
Oligenic inheritance is caused by mutations in two or more proteins with a functional relationship through direct interactions, membership in a pathway, or coexpression in a specific cell type. Given that functional groups of genes tend to colocalize within chromosomes (Thevenin, Ein-Dor, Ozery-Flato, & Shamir, 2014), the possibility of oligogenic inheritance of variants found within a locus defined by linkage analysis cannot be ignored. Our focused analyses of the DYT13 and DYT21 loci provide genes and variants for cosegregation analysis in these previously detailed dystonia pedigrees and suggest that digenic or higher-order oligogenic inheritance of variants within a disease-associated locus may be causal in some pedigrees and isolated cases with BSP. In this context, cosegregating variants in CIZ1 and SETX were linked to cervical dystonia in a large American pedigree .
Blepharospasm exerts important effects on health-related quality of life (Hall et al., 2006). Many patients with BSP experience annoying dry eye symptoms and photophobia (Hallett, Evinger, Jankovic, Stacy, & Workshop, 2008). Oral medications such as anticholinergics and benzodiazepines are mildly beneficial in some subjects. Many patients with BSP show moderate benefit from injections of botulinum toxin. However, injections are expensive, painful and may be denied by third-party payers. Although deep brain stimulation has been used to treat some individuals with BSP+ phenotypes, responses have been mixed (Reese et al., 2011). Major advances in the treatment of BSP demand a deeper understanding of its genetic etiopathogenesis. Our work provides a platform for follow-up case-control analyses of identified variants, evaluation of digenic and higher-order oliogenic etiologies for BSP (Deltas, 2017), and generation of animal models to help assess the pathogenicity of identified variants. Future work will demand attention to the effects of genetic background, oligogenic inheritance, pleiotropy, confounds of phenocopies, and the limitations of WES. TIAN ET AL.