Association of Rare Genetic Variants in Opioid Receptors with Tourette Syndrome.

Background: Genes involved in Tourette syndrome (TS) remain largely unknown. We aimed to identify genetic factors contributing to TS in a French cohort of 120 individuals using a combination of hypothesis-driven and exome-sequencing approaches.
Methods: We first sequenced exons of SLITRK1-6 and HDC in the TS cohort and subsequently sequenced the exome of 12 individuals harboring rare variants in these genes to find additional rare variants contributing to the disorder under the hypothesis of oligogenic inheritance. We further screened three candidate genes (OPRK1, PCDH10, and NTSR2) preferentially expressed in the basal ganglia, and three additional genes involved in neurotensin and opioid signaling (OPRM1, NTS, and NTSR1), and compared variant frequencies in TS patients and 788 matched control individuals. We also investigated the impact of altering the expression of Oprk1 in zebrafish.
Results: Thirteen ultrarare missense variants of SLITRK1-6 and HDC were identified in 12 patients. Exome sequencing in these patients revealed rare possibly deleterious variants in 3,041 genes, 54 of which were preferentially expressed in the basal ganglia. Comparison of variant frequencies altering selected candidate genes in TS and control individuals revealed an excess of potentially disrupting variants in OPRK1, encoding the opioid kappa receptor, in TS patients. Accordingly, we show that downregulation of the Oprk1 orthologue in zebrafish induces a hyperkinetic phenotype in early development. Discussion: These results support a heterogeneous and complex genetic etiology of TS, possibly involving rare variants altering the opioid pathway in some individuals, which could represent a novel therapeutic target in this disorder.

Introduction

Tourette syndrome (TS) is a neurodevelopmental disorder characterized by multiple involuntary motor and vocal tics that typically begin in childhood and have a waxing and waning course.[1-3] Twin studies have provided compelling evidence that TS is genetically determined,[4-7] but the nature of the involved genetic factors, their mode of inheritance, and the mechanisms by which they act remain largely unknown.

Recent studies have emphasized that the genetic architecture of TS is very likely oligo- (i.e., involving a few gene loci) or poly-genic (i.e., influenced by many genes as well as environmental factors) and overlaps with that of frequently comorbid disorders, such as autism spectrum disorders, obsessive-compulsive disorder (OCD), attention deficit hyperactivity disorder (ADHD), and depression.[8-10] Main genes possibly involved in TS so far include NRXN1, CNTN6, CELSR3, SLITRK1, and HDC.[11-22] SLITRK1 encodes a transmembrane protein modulating neurite outgrowth belonging to a family comprising six paralogues (SLITRK1 to SLITRK6) in mammals.[23,24] HDC encodes histidine decarboxylase, the enzyme catalyzing the synthesis of histamine through decarboxylation of histidine.[11] Slitrk1-, Slirk5-, and Hdc-deficient mice show behavioral disturbances and elevated anxiety-like behaviors, further supporting the involvement of these genes in neuropsychiatric disorders.[25-31]

In this study, we sequenced exons of SLITRK1-6 and HDC genes and analyzed the exome of 12 TS patients with rare variants in these genes to find additional hits contributing to the disorder under the hypothesis of oligogenic inheritance. We focused our analysis on rare, possibly disrupting variants in genes preferentially expressed in the basal ganglia, one of the main brain regions of interest in TS.[32] Screening of candidate genes in our TS cohort and comparison of variant frequencies in control populations identified OPRK1 variants as susceptibility factors for TS.

Methods

Participants

This study included 120 French individuals (92 males and 28 females, M/F sex ratio: 3.3) with TS syndrome recruited from the “Réseau National Gilles de la Tourette” (INSERM RBM 04-18; 2005–2007) and the “Centre de Référence National Maladies Rares Syndrome Gilles de la Tourette” (Pitié-Salpêtrière Hospital, Paris, France). Patient assessment is described in detail elsewhere.[33] Twelve patients with rare variants in SLITRK1-6 and HDC (see supplementary Table 1) from this cohort were identified and further analyzed by exome sequencing (Figure 1). One patient out of 12 (SGT-003-04) had pervasive developmental disorder not otherwise specified (PDD-NOS) in addition to TS. A total of 788 French healthy individuals matched on sex (632 males and 156 females, M/F sex ratio: 4.0) were selected from the same databank (DNA and cell bank of Institut du Cerveau et de la Moelle épinière [ICM]) to compare variant frequencies in NTS, OPRK1, and OPRM1. Genomic DNA was extracted from blood cells using standard phenol–chloroform procedures.

Figure 1

Flowchart Illustrating the Strategy and Results of this Study. The scheme summarized the different experimental steps and main results obtained.

Sanger sequencing

Specific primer pairs were designed to amplify the coding exons and adjacent intron–exon boundaries of SLITRK1-6 and HDC. The same primers were used to sequence the coding regions of these genes in the 120 TS patients. After exome sequencing, six genes were selected and sequenced by Sanger sequencing: OPRK1, OPRM1, PCDH10, NTS, NTSR1, and NTSR2. All primer sequences are provided in supplementary Table 2. Forward and reverse sequence reactions were performed with the Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA, USA). G50-purified sequence products were run on an ABI 3730 automated sequencer (PE Applied Biosystems) and data were analyzed with Seqscape v2.6 software (Applied Biosystems). In addition, 63 variants identified by exome sequencing were confirmed by Sanger sequencing (supplementary Table 3). For genes with multiple expressed isoforms, the longest isoform was used to describe variants. The variant nomenclature is based on the following cDNA reference sequences: OPRM1: NM_001145279 (longest isoform) and NM_000914.4 (main isoform); OPRK1: NM_000912.3; PCDH10: NM_032961.1; NTS: NM_006183.4; NTSR1: NM_002531.2; NTSR2: NM_012344.3.

Exome sequencing

We sequenced the exome of 12 patients with variants in SLITRK1-6 and HDC (Figure 1). Genomic DNA was captured using the SureSelect Human All Exon 50 Mb Kit (Agilent, Santa Clara, CA, USA) followed by paired-end 75 bp massively parallel sequencing on a HiSeq2000-sequencing system (Illumina, San Diego, CA), according to the manufacturer’s instruction and protocols,[34] by Integragen SA (Evry, France). Image analysis and base calling were performed using Illumina Real Time Analysis (RTA) Pipeline version 1.8 with default parameters. Analysis of exome data was based on the Illumina CASAVA1.7 pipeline, using the ELANDv2 alignment algorithm. Variant annotation was realized with an in-house pipeline (Integragen Evry, France), consisting of gene annotation (RefSeq), frequencies of the variants at the homozygous or heterozygous state in databases (HapMap, 1,000 genomes, exome variant server [EVS]) and in exomes sequenced at Integragen, and analysis of the impact of the variants (exonic, intronic, silent, nonsense, etc.). Variant filtering was performed using the Eris interface (Integragen) with the following parameters: (1) minor allele frequencies (MAFs) ≤0.01 (rare heterozygous variants) or 0.001 (homozygous variants) in dbSNP, Hapmap, 1,000 genomes, AVS, and a local exome database; (2) a predictable damaging impact on the gene or the protein (nonsense variants, variants predicted to cause a frameshift, variants altering splice sites, missense variants, and variants introducing in-frame deletions or insertions); (3) coverage ≥10 reads and exclusion of ambiguous genotypes (supplementary Table 3). This list was then filtered to keep genes with variants (same or different) in at least two TS patients (supplementary Table 4). These parameters allowed to detect all 13 variants previously identified in the HDC and SLITRK1-6 genes. For six patients, DNA from both parents was available and their exome was also sequenced to investigate de novo variants. For one family (SGT-003), parents were included in the first analysis. For the five remaining families, the index case and the parents were resequenced independently and these data were used to validate variants identified in the first experiment. Exome sequencing in family SGT-082 revealed an unexplained incompatibility of the paternal DNA.

Microarrays

Copy number variants (CNVs) were independently assessed in 11/12 index cases with HDC and SLITRK1-6 variants using OmniExpress24 microarrays (Illumina). Automated Illumina microarray experiments were performed according to the manufacturer’s specifications. Image acquisition was performed using an iScan System (Illumina, San Diego, CA, USA). Image analysis and automated CNV calling were performed using GenomeStudio v2011.1 and CNVPartition v3.1.6 with the default confidence threshold of 35. CNVs identified in TS patients were compared to those present in the Database of Genomic Variants (DGV), using Cartagenia Bench (Agilent). Variants absent from DGV, present in less than three individuals, or with an overlap <80% with CNVs present in DGV were listed (supplementary Table 5).

Pyrosequencing

To compare the frequencies of rare variants in NTS, OPRM1, and OPRK1, exons and intron–exon junctions of these genes were screened in 120 patients with TS and 788 healthy controls using amplicon-based massive parallel pyrosequencing with a Roche 454 GS Junior. A first polymerase chain reaction (PCR) was performed to amplify the targeted exons by PCR (primers in supplementary Table 2). A second PCR was performed to incorporate multiplex identifier and 454 adaptors, and an emulsion PCR was carried out according to the emPCR Amplification Method Manual (Roche, Basel, Switzerland). All variants previously identified by Sanger sequencing were also detected by pyrosequencing in the TS cohort.

Bioinformatic analyses

Variants identified in TS individuals were annotated using Alamut 2.7.1 (Interactive Biosoftware, Rouen, France). This software allows automatic gene variant annotation including missense predictions using Align GVGD, SIFT, Mutation Taster and Polyphen-2. Rare variants correspond to variant with MAF<0.01; ultrarare variants have an MAF<0.001. We used the statistical overrepresentation tests from PANTHER 14.1[35] with the “GO biological process complete” option to look for molecular pathways and functions enriched in gene lists with rare variants identified by exome sequencing in TS patients (supplementary Tables 6 and 7).

Median Reads Per Kilobase Million (RPKM) values for each human gene GTEx Analysis V6p were downloaded from the Genotype-Tissue Expression (GTEx) portal (https://www.gtexportal.org/). Genes preferentially expressed in the adult basal ganglia or substantia nigra (maximum RPKM value in caudate, nucleus accumbens, putamen, or substantia nigra>maximum RPKM value in other tissues/brain regions) were listed and ranked according to their expression in basal ganglia/substantia nigra compared to other tissues. The final list includes 491 transcripts with RPKM>0.5 (supplementary Table 8). We then listed the 64 variants present in the 12 individuals in basal ganglia-enriched genes (supplementary Table 9). Gene annotation (gene function, gene expression, related phenotypes, or disorders in mice and humans) was performed using GeneALaCarte (http://www.genecards.org/cgi-bin/BatchQueries/Batch.pl) and MGI Batch Query (http://www.informatics.jax.org/batch).

Missense and Loss-of-function (LoF) variants present in PCDH10, NTS, NTSR1, NTSR2, OPRK1, and OPRM1 present in the gnomAD[36] r2.0.1 database (http://gnomad.broadinstitute.org/) were downloaded. In silico pathogenicity predictions calculated by PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/) were systematically retrieved for missense variants present in gnomAD using the batch query tool (http://genetics.bwh.harvard.edu/pph2/bgi.shtml). Possibly or probably damaging predictions by Polyphen-2, which usually are the most stringent of all four predictions tools, were used to filter possibly damaging variants and compare variant frequencies (supplementary Table 11). Variant frequencies were compared using exact Fisher’s tests.

Zebrafish studies

Adult and larval zebrafish (Danio rerio) were maintained in the animal platform of Inserm Unit U1127, ICM in Paris, according to the National and European Guidelines for Animal Welfare. Experiments were performed on wild type embryos from AB and TL strains. Zebrafish larvae were raised in embryo water: 0.6 g/L aquarium salt (Instant Ocean, Blacksburg, VA) in reverse osmosis water + 0.01 mg/L methylene blue in a 28°C incubator. Microinjections of 0.5 nL of solution were carried out at one-cell stage of zebrafish embryos. Antisense Morpholino oligonucleotides (AMOs) (GeneTools Philomath, OR) were used to knock down the expression of the sole orthologue of the opioid receptor kappa 1 (Oprk1: NM_182886) gene in zebrafish. The sequence of the ATG-targeting AMO was the following AAATCTGCACCACGTTGCTGTCCA that would bind to the following genomic sequence in the zebrafish genome: ACAGGCATCACTATTGTTTTCAAACCTTG-GGATTTTCGTCTTAAACTACCTGTTTAAAATTCACAACAAGACG[(ATG)GACAGCAACGTGGTGCAGATTT]. The mismatch control for this AMO was the following (small letters scrambled to disrupt binding) AAATaTGaACCACaTTGCTaTCaA. To ensure the AMO specificity, we also utilized a splice AMO targeting the exon2–intron2 junction with the following sequence: AGGACAATGCTGTTTCTTACCTGAT. Mismatch control for the splice AMO was the following: AGcACAATcCTcTTTCTTACgTcAT. All AMOs were used at a concentration of 0.6 mM in the injection solution for a total amount of 2.5 ng of AMO injected per egg. For rescue experiments, the human wild-type OPRK1 cDNA construct was obtained from Dharmacon, Lafayette, CO, USA, cloned into a pCS2 plasmid, and injected alongside the AMO at a concentration of 75 ng/μL. Total RNA from embryos was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA). RNA was quantified using the Nanodrop 8000 (Thermo Scientific, Waltham, MA). cDNA was synthesized using Transcriptor Universal cDNA Master Mix (Roche) with the following primer pairs to detect the expression and splicing of the OPRK1 transcript: forward: 5′-CATCACTGCCGTTTACTCCGT and reverse 5′-TCATTGGTGTCCGGAAGTCC. To record global activity, 28 hpf zebrafish embryos in their chorion were placed on a plastic mesh, submerged in embryo water, and imaged at the frequency of 30 fps as previously described.[37] Activity parameters were quantified using an automatized Zebralab system (ViewPoint, France). Analysis of variance (ANOVA) was performed using the SigmaStat software (California) and data were reported as significant when p<0.05.

Results

To study the contribution of variants in HDC and SLITRK genes in the French TS cohort, coding exons of SLITRK1, SLITRK5, SLITRK6 (chromosome [chr] 13), SLITRK3 (chr3), SLITRK2, SLITRK4 (chrX), and HDC (chr15) were screened in 120 TS patients by Sanger sequencing (Figure 1). Seventeen rare heterozygous variants (MAF<0.01) were identified in 19 unrelated index probands (supplementary Table 1). Among these, 12 individuals carried 13 ultrarare (MAF<0.001) variants modifying the coding sequence (supplementary Table 1, supplementary Figure 1). All variants but one (HDC: c.1492A>G, p.Arg498Gly), which was present in an individual who also had a variant in SLITRK6 (c.1232C>G, p.Thr411Arg), alter amino acids conserved in mammals (supplementary Figure 2). One of the three SLITRK1 variants (c.1252 A>T, p.Thr418Ser) was previously described in patients with OCD.[38] Seven variants for which inheritance could be assessed were inherited from asymptomatic parents, indicating that these variants, if involved in the etiology of TS, are insufficient on their own to lead to the disorder.

To find additional hits contributing to the disorders under the assumption of oligogenic inheritance, we sequenced the exome of these 12 patients (Figure 1). We identified a total of 3,755 rare variants (mean number of variants = 408 per patient; range: 342–595) altering coding sequences or splice sites in 3,041 genes (supplementary Table 3). A set of 63 variants in selected candidate genes was validated by Sanger sequencing in the affected subjects, confirming the presence of all variants but one (supplementary Table 3). At least two different variants were identified in two or more individuals in 424 genes (supplementary Table 4). One individual (SGT-040-10) for whom DNA from both parents was available had two de novo variants absent or present at a very low frequency in gnomAD: a de novo duplication (arr[hg19] 17q12(37,042,006-37,276,683) × 3) on chr 17q12 including three genes (FBXO47, LASP1, PLXDC1) and a de novo missense variant (NM_005739.3: c.2327T>C, p.Ile776Thr) in the RASGRP1 gene (supplementary Table 5).

We used PANTHER to determine whether specific biological processes were enriched in the 3,041 and 424 gene lists, respectively. Interestingly, both lists were enriched in genes involved in cell differentiation, cell adhesion, cell projection organization, neurogenesis, and neuron differentiation among other cellular processes (supplementary Tables 6 and 7). Genes involved in these pathways include NRP2 (three individuals), ROBO4 (two individuals), SLIT2, NTN5, RELN, RGMA, plexins (PLXNA1, PLXNB1), semaphorins (SEMA4F, SEMA6B: two individuals each, SEMA3B, SEMA3C, SEMA4B, SEMA5A, SEMA5B, SEMA6D), ephrins (EPHA1, EPHA6: two individuals each, EPHA8, EPHA10, EPHB4), cadherins (CDH4: three individuals, CDH20), protocadherins (PCDH10), cadherin EGF LAG Seven-Pass G-type receptors (CELSR3: two individuals; CELSR2), and fat-like cadherins (FAT4: four individuals; FAT3: three individuals; FAT1; two individuals; FAT2).

We then examined variants in genes preferentially expressed in the basal ganglia (supplementary Table 8), corresponding to a total of 64 variants in 54 genes (supplementary Table 9). Among these variants, two in PCDH10/protocadherin 10 (c.3076C>T; p.Arg1026*) and OPRK1/opioid receptor kappa-1 (c.787C>T; p.Arg263*) lead to a theoretical protein truncation. A missense variant in the neurotensin receptor-2 (NTSR2, c.397C>G, p.Arg133Gly) altering a highly conserved amino acid (until C. elegans) was predicted to be damaging by all four prediction tools used.

Opioid and neurotensin signaling pathways being relevant to TS etiology, we examined variants in genes encoding ligand or receptors independently of their preferential expression in the brain. We observed rare variants possibly altering other genes involved in opioid (c.743C>T, p.Thr248Ile in OPRM1/opioid receptor Mu-1) and neurotensin (c.423G>C, p.Lys141Asn in NTS/neurotensin, and c.533C>T, p.Thr178Ile in NTSR1/neurotensin receptor-1) pathways in three additional patients.

We then screened NTS, NTSR1, NTSR2, OPRK1, OPRM1, and PCDH10 by Sanger sequencing for additional variants in our TS cohort. Additional rare variants were identified in all six genes, including PCDH10 missense variants (c.1526G>A, p.Gly509Asp and c.2411T>C, p.Ile804Thr) in three individuals, one frameshift (c.348del, p.Gln117Asnfs*5) and one missense variant (c.190C>A, p.Leu64Ile) in NTS in two individuals, rare variants in NTSR1 and NTSR2 in 11 and 3 individuals each, one nonsense (c.1107C>A, p.Tyr369*) and three missense variants in OPRK1 (four individuals), and eight rare variants including a nonsense variant (c.829C>T, p.Arg277*) in OPRM1 (supplementary Table 10).

Comparison of rare variant frequencies between TS cases and individuals from the gnomAD database suggested that variants in OPRK1 may be more frequent in TS patients than in control individuals (TS: 4.2%, controls: 0.8%; p = 0.003, Fisher’s exact test, Table 1, supplementary Tables 10 and 11). To further confirm these findings, we reanalyzed TS cases and a control population comprising 788 individuals matched on sex from the same databank for variants in OPRK1 using pyrosequencing (supplementary Table 10). We also included NTS and OPRM1 in this screening based on the identification of truncating variants in these genes. Comparison of variant frequencies in TS cases and matched controls confirmed that there was a significant excess of rare variants in OPRK1 (TS: 4.2%, controls: 0.5%; p = 0.003, Fisher’s exact test), and to a lesser extend in OPRM1 (TS: 6.7%, controls: 2.3%; p = 0.048, Fisher’s exact test) in TS patients, while the difference for NTS (TS: 2.5%, controls: 1.2%; p = 0.397, Fisher’s exact test) was not significant (Table 1). Variants of OPRK1 and OPRM1 in which inheritance could be assessed were all inherited from an asymptomatic parent. In two families including affected siblings, the variant was also present in the affected sib (supplementary Table 3, supplementary Figure 3). These results suggest that variants in opioid receptors could represent susceptibility factors to TS.

Table 1

Comparison of the Frequencies of Rare Variants in Six Selected Genes in TS and Control Populations

Gene	Type of variants	TS cohort (120 individuals)	gnomAD (138,632 individuals)	p	Control individuals (788 individuals)	p
OPRK1*	Missense/LoF (MAF<1%)	5 (4.2%)	1099 (0.8%)	0.003	4 (0.5%)	0.003
	Damaging/LoF (MAF<1%)	4 (3.3%)	541 (0.4%)	0.001	3 (0.4%)	0.007
OPRM1**	Missense/LoF (MAF<1%)	8 (6.7%)	8,085 (5.8%)	0.694	22 (2.3%)	0.048
	Damaging/LoF (MAF<1%)	5 (4.2%)	3,786 (2.7%)	0.268	10 (1.2%)	0.037
NTS	Missense/LoF (MAF<1%)	3 (2.5%)	998 (0.7%)	0.057	10 (1.2%)	0.397
	Damaging/LoF (MAF<1%)	2 (1.7%)	474 (0.3%)	0.064	2 (0.25%)	0.087
NTSR1	Missense/LoF (MAF<1%)	12 (10.0%)	14,333 (10.3%)	1.000	ND
	Damaging/LoF (MAF<1%)	6 (5.0%)	8,248 (5.9%)	0.847	ND
NTSR2	Missense/LoF (MAF<1%)	4 (3.3%)	3,142 (2.3%)	0.353	ND
	Damaging/LoF (MAF<1%)	4 (3.3%)	1,772 (1.3%)	0.069	ND
PCDH10	Missense/LoF (MAF<1%)	4 (3.3%)	2,494 (1.8%)	0.171	ND
	Damaging/LoF (MAF<1%)	3 (2.5%)	1,174 (0.8%)	0.083	ND

The details of variants included in the comparisons appear in supplementary Tables 10 and 11. The number of variants/individuals calculated from gnomAD (supplementary Table 11) assumed that rare variants correspond to one individual each. Variant frequencies were compared with Fisher’s exact tests.

Calculations for OPRK1 excluded Asp374Asp, which is frequent in the African population (longest isoform NM_001318497).

For OPRM1, only variants in exons 3, 4, and 5 of isoform NM_001145279, which are present in at least 11/19 isoforms, were taken into account.

We next tested the functional consequences of altering the expression of the TS candidate orthologous genes during vertebrate development. For this, we focused on the expression of oprk1 in zebrafish, a simple vertebrate model that exhibits stereotyped movements during early development.[39] The zebrafish genome contains a single orthologue of the human OPRK1 gene (ensembl: oprk1 ENSDARG00000006894) with 67% identity at the protein level. Oprk1 expression is enriched in brain structures of the diencephalon and mesencephalon in adult zebrafish,[40] and is expressed during the development of the embryo.[41] To investigate the function of Oprk1 during early development, we performed a knockdown of its expression via microinjections of morpholino-modified antisense oligonucleotides (AMO) in the one-cell stage embryo. We used two different strategies to knock down the expression of oprk1, targeting either the ATG starting codon or the splice junction exons 2/intron 2. Injection of 2.5 ng of either AMOs did not result in gross morphological defects at 28 hpf (Figure 2A) or at later stages of development (not shown) when compared with non-injected or with mismatch AMO-injected controls. To determine the effectiveness of the splice-modifying AMO, we performed reverse transcription PCR using primers flanking the exon 2–exon 3 junction. Microfluidic capillary electrophoresis of the PCR products revealed a defective splicing event in fish injected with the splice-blocking AMO (Figure 2B). Further, a lower amount of PCR product was detected in the ATG-targeting AMO injected fish, indicating a degradation of the oprk1 transcript (Figure 2B).

Figure 2

Knockdown of the Zebrafish OPRK1 Orthologue Results in Hyperactivity during Development. (A) Activity heatmap illustrating spontaneous movement of 28 hpf embryos inside the chorion during a 10 seconds period. Higher levels of activity are evident in embryos injected with oprk1 splice-blocking AMO (spOPRK1) or ATG-targeting AMO (OPRK1 KD) when compared to mismatch AMO-injected embryos (Mismatch control). Expression of human OPRK1 transcript (hOPRK1) rescues the hyperactivity phenotype. (B) RT-PCR product of amplification of a region spanning exons 2 and 3 of the zebrafish oprk1 transcript. Injection of the splice-blocking AMO (spOPRK1) results in a defective splicing in this region. (C) Quantification of the total activity of 28 hpf zebrafish embryos. (D) Quantification of the frequency of movement bursts in the 28 hpf embryos.

To determine the effect of oprk1 knockdown on the motor activity, we analyzed the movement of zebrafish embryos at 28 hpf. This developmental stage marks the appearance of the first stereotyped movement, consisting of bidirectional periodic coils inside the chorion. We quantified this type of movement using an automated video capture and analysis system. As illustrated in Figure 2A, knockdown of oprk1 was associated with an overall increase in activity in zebrafish. Indeed, injection of either the ATG-targeting or the splice blocking AMO resulted in significantly higher total movement (Figure 2C) and frequency of coiling bursts (Figure 2D). To confirm that the hyperkinetic phenotype is specifically associated with a loss of function of oprk1, we coexpressed human OPRK1 together with the ATG-blocking AMO. As shown in Figure 2, overexpression of OPRK1 functionally rescued the knockdown of the zebrafish orthologue, as both the total activity (Figure 2C) and the frequency of coiling movements (Figure 2D) were restored to control levels. Investigation of the global activity of zebrafish at later stages of development (48 hpf–5 dpf) revealed a gradual disappearance of the hyperactivity phenotype, with no significant differences from controls (not shown). Thus, in zebrafish, we observed an association of oprk1 loss of function with a transient hyperkinetic phenotype at early stages of development.

Discussion

Genetics of TS is both complex and heterogeneous, making the identification of interacting variants extremely challenging. In an attempt to circumvent this difficulty, we focused on candidate genes enriched in brain structures relevant to TS pathophysiology, that is, the basal ganglia. Although straightforward, this approach would miss variants in genes potentially involved in TS-related networks dysfunction but rather mainly expressed in other brain regions or tissues.

The main finding of this study is an excess of possibly disrupting variants in OPRK1, the gene encoding the kappa opioid receptor (KOR) in patients with TS. We also detected a possible, weaker association of TS with rare variants in OPRM1, although these results have to be interpreted with caution, notably because OPRM1 is a very polymorphic gene, encoding many different isoforms. We also show that knockdown of the oprk1 orthologue in zebrafish is associated with a hyperkinetic phenotype that is confined to the early stages of development. Overexpression of the human OPRK1 transcript was able to rescue the loss of function of the zebrafish orthologue, highlighting the high functional conservation of the opioid receptor kappa during evolution and the relevance of its investigation in lower vertebrates. Altogether, these results further support the hypothesis that alteration of opioid receptor kappa in vivo could lead to an early transient hyperactivity phenotype mimicking TS. Extending this study to the functional analysis of variants identified in TS patients may also provide additional support for the potential role of OPRK1 in TS pathogenesis.

The endogenous opioid system is composed of three different opioid receptor subtypes (mu, kappa, and delta), which have different function and brain distribution.[42] These receptors are activated upon binding with endogenous peptides including dynorphins, enkephalins, endorphins, endomorphins and nociceptin, or exogenous drugs such as morphine or heroin. The opioid endogenous system, and more specifically the KOR system,[43] is involved in both motor control and psychiatric conditions such as anxiety, depression, schizophrenia, and addiction.[44] Sex differences in KOR function[45] might contribute to sex-dependent clinical manifestations or to imbalanced sex ratios existing in TS and many other neuropsychiatric conditions such as susceptibility to develop addictions and depression.

The implication of the opioid system in pathophysiology of TS has been hypothesized for more than 30 years, first based on treatment studies using opiate agonists and/or antagonists.[46-50] Importantly, a postmortem report of a single patient with TS showed total absence of dynorphin-like staining in the dorsal part of the globus pallidus externus and very faint staining in the ventral pallidum.[51] This case was particularly interesting since it pinpointed toward underactivity of a specific opioid system, the dynorphin/kappa receptor system. Very recently and in line with our own findings, downregulation of OPRK1 and OPRM1 opioid receptor gene expression in the postmortem striatum of nine patients with TS was reported.[52] Frequent variants in OPRK1 (rs702764) and OPRM1 (rs1799971) have also recently been associated with impulse control disorders in Parkinson disease,[53,54] a finding in agreement with genetic variants in opioid receptor contributing to impulsivity disorders, a comorbidity frequently associated with TS.[55]

With relevance to TS, KORs are located on dopaminergic afferents where they reduce release and increase the reuptake of dopamine,[56] implying that inactivation of KORs (by OPRK1 variants, for instance) may increase cerebral dopamine release, especially in regions such as the striatum.[57] This would be in line with the predominant theory that TS pathophysiology is associated with a dopaminergic hyperinnervation.[58] Little is known about the specific function of oprk1 during zebrafish development; however, a role in the differentiation of Pax2-positive neurons has been suggested,[41] with implications for the correct specification of the dopaminergic system.

Contrary to Oprk1 knockdown in zebrafish, Oprk1-deficient mice did not show evidence of hyperactivity.[59] Oprk1 knockdown in zebrafish and Oprk1 knockout in mice cannot be strictly compared as compensation mechanisms are very likely different when Oprk1 is downregulated or completely absent during development in these distinct species. Furthermore, locomotor activity is usually tested in adult mice and an early transient hyperactive phenotype could have been easily missed. Furthermore, locomotor phenotypes and behaviors also depend on the genetic background, as observed for Oprd1-deficient mice, which mainly showed impulsivity on a C57Bl/6J background.[60,61] However, a recent study has shown that triple knockout of all three opioid receptors in mice is associated with hyperactivity (motor restlessness) and dopamine abnormalities.[62]

The number of patients with opioid receptor variants in our TS cohort remains low (3.3–4.2% compared to 0.5–0.8% in control individuals). It is also important to note that a series of 120 individuals is small for case-control association studies, and that TS and control populations have not been matched for ancestry based on genotypes. Therefore, it would be important to replicate these findings in larger TS and control cohorts. The impact of population stratification on frequencies of variants in opioid receptors is illustrated by the OPRM1 NM_001145279:c.397A>G, p.Asn133Asp (NM_000914.4:c.118A>G, p.Asn40Asp, rs1799971) variant, a common functional polymorphism that was shown to act as a susceptibility factor for addiction and to alter the response treatments targeting opioid receptors.[63] This variant was identified in 31 TS patients (28 heterozygotes and 3 homozygotes), a distribution similar to that observed in our control population as well as in the European (non-Finnish) population of gnomAD, suggesting that this polymorphism does not play a major role in TS etiology. However, we observed that the frequency of the G minor allele significantly varies across populations (much higher in Asian populations while lower in Africans; supplementary Table 12). This finding might have important consequences when considering opioid receptors as possible therapeutic targets in TS or other psychiatric disorders.

KOR agonists – assuming that residual receptor function is maintained – might indeed be novel therapeutic options to be investigated in TS. However, classic KOR agonists such as U50,488 or U69,593 used experimentally to treat drug addiction have shown dysphoric and psychomimetic properties. Therefore, novel molecules based on Salvinorin A, the active compound of Salvia divinorum, are currently being developed for the treatment of substance abuse (particularly cocaine) but might also prove to be useful for the treatment of tics.[64]

Finally, this work has also identified rare variants in several individuals in additional candidate genes including PCDH10, encoding a protocadherin with preferential expression in basal ganglia, possibly contributing to the establishment of synaptic connections in this brain region; ROBO4, a gene of the ROBO/SLIT axon guidance pathway previously associated with TS[65]; NRP2, encoding neuropilin 2, a protein possibly controlling axon guidance of dopaminergic neurons[66]; NOS1, encoding the neuronal nitric oxide synthase, a gene also downregulated in the striatum of TS patients[52] and possibly playing a role in controlling motor behavior;[67] and SHANK1, in which truncating variants were previously associated with male-limited Asperger syndrome.[68] Two individuals also had variants in CELSR3, a gene that was recently recognized as a candidate gene for TS.[22] The possible contribution of these genes to TS etiology needs to be confirmed in further studies, as are the nature and number of the variants acting in concert necessary to lead to TS.