A Point Mutation in a lincRNA Upstream of GDNF Is Associated to a Canine Insensitivity to Pain: A Spontaneous Model for Human Sensory Neuropathies

Abstract

Human Hereditary Sensory Autonomic Neuropathies (HSANs) are characterized by insensitivity to pain, sometimes combined with self-mutilation. Strikingly, several sporting dog breeds are particularly affected by such neuropathies. Clinical signs appear in young puppies and consist of acral analgesia, with or without sudden intense licking, biting and severe self-mutilation of the feet, whereas proprioception, motor abilities and spinal reflexes remain intact. Through a Genome Wide Association Study (GWAS) with 24 affected and 30 unaffected sporting dogs using the Canine HD 170K SNP array (Illumina), we identified a 1.8 Mb homozygous locus on canine chromosome 4 (adj. p-val = 2.5x10-6). Targeted high-throughput sequencing of this locus in 4 affected and 4 unaffected dogs identified 478 variants. Only one variant perfectly segregated with the expected recessive inheritance in 300 sporting dogs of known clinical status, while it was never present in 900 unaffected dogs from 130 other breeds. This variant, located 90 kb upstream of the GDNF gene, a highly relevant neurotrophic factor candidate gene, lies in a long intergenic non-coding RNAs (lincRNA), GDNF-AS. Using human comparative genomic analysis, we observed that the canine variant maps onto an enhancer element. Quantitative RT-PCR of dorsal root ganglia RNAs of affected dogs showed a significant decrease of both GDNF mRNA and GDNF-AS expression levels (respectively 60% and 80%), as compared to unaffected dogs. We thus performed gel shift assays (EMSA) that reveal that the canine variant significantly alters the binding of regulatory elements. Altogether, these results allowed the identification in dogs of GDNF as a relevant candidate for human HSAN and insensitivity to pain, but also shed light on the regulation of GDNF transcription. Finally, such results allow proposing these sporting dog breeds as natural models for clinical trials with a double benefit for human and veterinary medicine.

About the Authors:

Jocelyn Plassais

* E-mail: [email protected] (CA); [email protected] (JP)

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Laetitia Lagoutte

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Solenne Correard

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Manon Paradis

Affiliation: Department of Clinical Sciences, Faculté de Médecine Vétérinaire, University of Montreal, Montreal, Québec, Canada

Eric Guaguère

Affiliation: Clinique Vétérinaire Saint Bernard, Lomme, France

Benoit Hédan

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Alix Pommier

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

ORCID http://orcid.org/0000-0003-0588-8099

Nadine Botherel

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Marie-Christine Cadiergues

Affiliation: UDEAR, Université de Toulouse, Ecole Nationale Vétérinaire de Toulouse, INSERM, Toulouse, France

Philippe Pilorge

Affiliation: Clinique Vétérinaire de Saint-Cyr, Rennes, France

David Silversides

Affiliation: Department of Clinical Sciences, Faculté de Médecine Vétérinaire, University of Montreal, Montreal, Québec, Canada

Maud Bizot

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Mark Samuels

Affiliation: Department of Biochemistry and Molecular Medicine, CHU Sainte-Justine, University of Montreal, Montreal, Québec, Canada

Carme Arnan

Affiliations Centre for Genomic Regulation (CRG), Barcelona, Spain, Universitat Pompeu Fabra (UPF), Barcelona, Spain, Institut Hospital del Mar d'Investigations Mèdiques (IMIM), Barcelona, Spain

Rory Johnson

Affiliations Centre for Genomic Regulation (CRG), Barcelona, Spain, Universitat Pompeu Fabra (UPF), Barcelona, Spain, Institut Hospital del Mar d'Investigations Mèdiques (IMIM), Barcelona, Spain

ORCID http://orcid.org/0000-0003-4607-2782

Christophe Hitte

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Gilles Salbert

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Agnès Méreau

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Pascale Quignon

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Thomas Derrien

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, France

Catherine André

* E-mail: [email protected] (CA); [email protected] (JP)

Affiliations CNRS, UMR 6290, Institut de Génétique et Développement de Rennes, Rennes, France, Université Rennes 1, UEB, Biosit, Faculté de Médecine, Rennes, FranceAbstract

Human Hereditary Sensory Autonomic Neuropathies (HSANs) are characterized by insensitivity to pain, sometimes combined with self-mutilation. Strikingly, several sporting dog breeds are particularly affected by such neuropathies. Clinical signs appear in young puppies and consist of acral analgesia, with or without sudden intense licking, biting and severe self-mutilation of the feet, whereas proprioception, motor abilities and spinal reflexes remain intact. Through a Genome Wide Association Study (GWAS) with 24 affected and 30 unaffected sporting dogs using the Canine HD 170K SNP array (Illumina), we identified a 1.8 Mb homozygous locus on canine chromosome 4 (adj. p-val = 2.5x10-6). Targeted high-throughput sequencing of this locus in 4 affected and 4 unaffected dogs identified 478 variants. Only one variant perfectly segregated with the expected recessive inheritance in 300 sporting dogs of known clinical status, while it was never present in 900 unaffected dogs from 130 other breeds. This variant, located 90 kb upstream of the GDNF gene, a highly relevant neurotrophic factor candidate gene, lies in a long intergenic non-coding RNAs (lincRNA), GDNF-AS. Using human comparative genomic analysis, we observed that the canine variant maps onto an enhancer element. Quantitative RT-PCR of dorsal root ganglia RNAs of affected dogs showed a significant decrease of both GDNF mRNA and GDNF-AS expression levels (respectively 60% and 80%), as compared to unaffected dogs. We thus performed gel shift assays (EMSA) that reveal that the canine variant significantly alters the binding of regulatory elements. Altogether, these results allowed the identification in dogs of GDNF as a relevant candidate for human HSAN and insensitivity to pain, but also shed light on the regulation of GDNF transcription. Finally, such results allow proposing these sporting dog breeds as natural models for clinical trials with a double benefit for human and veterinary medicine.

Author Summary

In this study, we present a canine neuropathy characterized by insensitivity to pain in the feet, sometimes combined with self-mutilation described in four sporting breeds. This particular phenotype has the clinical hallmarks of human Hereditary Sensory Autonomic Neuropathies (HSAN). As we hypothesized that a monogenic recessive disorder was shared between these breeds, we performed a Genome Wide Association Study (GWAS) to search for the genetic causes and found one homozygous chromosomal region in affected dogs. High-throughput sequencing of this region allowed the identification of a point mutation upstream to the GDNF gene and located in the last exon of a long non-coding RNA, GDNF-AS. We confirmed the perfect association of this variant with the disease using more than 900 unaffected dogs that do not present with this mutation. Functional analyses (qRT-PCR, EMSA) confirmed that the mutation alters the binding of regulatory complex, leading to a significant decrease of both GDNF and GDNF-AS mRNA expression levels. This work in canine spontaneous forms of human neuropathies allowed the identification of a novel gene GDNF and its regulation mechanism, not yet described in human HSAN, opening the field of clinical trials to benefit both canine and human medicine.

Citation: Plassais J, Lagoutte L, Correard S, Paradis M, Guaguère E, Hédan B, et al. (2016) A Point Mutation in a lincRNA Upstream of GDNF Is Associated to a Canine Insensitivity to Pain: A Spontaneous Model for Human Sensory Neuropathies. PLoS Genet 12(12): e1006482. doi:10.1371/journal.pgen.1006482

Editor: Adam R. Boyko, Cornell University College of Veterinary Medicine, UNITED STATES

Received: May 20, 2016; Accepted: November 15, 2016; Published: December 29, 2016

Copyright: © 2016 Plassais et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Raw data are available on the following confidential URL Genotyping: http://tools.genouest.org/data/tderrien/data/neuropath/data.map.txt http://tools.genouest.org/data/tderrien/data/neuropath/data.ped.txt Sequences produced by Targeted Sequencing: http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS01_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS01_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS02_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS02_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS03_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS03_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS04_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS04_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS05_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS05_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS06_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS06_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS07_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS07_R2_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS08_R1_fastq.txt.gz http://tools.genouest.org/data/tderrien/data/neuropath/Sample_NGSAMS08_R2_fastq.txt.gz

Funding: This study was supported by the CNRS (Centre National de la Recherche Scientifique), the Brittany Region (France) (PhD funding for JP), the European Commission (FP7-LUPA, GA-201370), the Rosembloom Family and the Companion animal health fund from the Faculté de Médecine Vétérinaire, Université de Montréal, and the CRB-Anim infrastructure, ANR-11-INBS-0003, funded by the French National Research Agency in the frame of the 'Investing for the Future' program.The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Inherited peripheral neuropathies are neurodegenerative diseases of the peripheral nervous system (PNS), including the Charcot-Marie-Tooth diseases (CMT) and Hereditary Sensory and Autonomic Neuropathies (HSANs) [1]. They are categorized based on the degree of involvement of motor, sensory and/or autonomic nerve fibers. HSANs are rare peripheral neuropathies constituting a clinically and genetically heterogeneous group which phenotypes range from pure sensory involvement to variable levels of motor and autonomic disturbances [2]. The HSAN types can be subdivided into two groups based on their mode of inheritance, HSAN I being autosomal dominant disorders and HSAN II to V being recessive. The main symptoms correspond to a progressive degeneration of sensory and autonomic neurons causing insensitivity to pain and temperature in feet and hands, leading to ulcerative mutilations [3-5]. These lesions could result in severe tissue infections and/or osteomyelitis and may lead to amputation of the affected limbs. To date, the pathophysiologic mechanisms known to be implicated in HSAN include sphingolipid metabolism with mutations in SPTCL1 or SPTCL2 (serine palmitoyl-transferase genes), vesicular transport (RAB7, IKBKAP) or neurotrophic factors and their tyrosine kinase receptors (NTRK1, NGF-B) [6-8]. Although 12 genes have been associated with HSANs until now, they only explain one third of the patients affected with HSANs, highlighting the fact that the genetic causes of HSANs remain only partially understood in humans [2] and that natural models of these diseases are needed.

Interestingly, while neurological diseases such as HSANs are rare in humans, some dogs breeds spontaneously develop such diseases with high frequencies [9]. Indeed, recent breeding practices have created genetically isolated populations with large homozygous haplotypes and linkage disequilibrium (LD) blocks [10-12]. With this particular genetic structure of population, the dog model allowed the identification of new genes involved in rare diseases [13,14], and thus constitutes an excellent model to study sensory neuropathies. Indeed, a large number of inherited motor and sensory neuropathies were reported in the dog species [9,15]. Since 1960's, symptoms similar to human HSANs were described in four sporting breeds, with the first cases reported in German short-haired Pointers (GP) in 1964 [16,17], then in English Pointers (EP), English Springer Spaniel (ESS) and more recently in French Spaniel (FS) [16,18-21]. For these four breeds, the hallmarks of the disease are loss of pain sensation and feet ulceration. Symptoms naturally occur in puppies approximately four months old when they begin to lick and bite their paws. Affected dogs present an acral insensitivity to pain with, in the majority of cases, severe self-mutilations of the feet including claw loss, painless fractures, and digit amputation (Fig 1A). No sign of autonomic disturbances are detected and motor abilities and proprioception are normal. The pathologic process thus seems to only affect sensory neurons, primarily their development, followed by progressive postnatal degeneration of unmyelinated fibers [18-20]. In addition, pedigree analyses in these breeds revealed a recessive autosomal mode of inheritance such as in human type II-V HSANs [2,21].

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 1. Association study for the acral mutilation syndrome in sporting dogs.

(A) Clinical phenotype in the French spaniel breed (left picture). Dogs present insensitivity to pain in feet and can sometimes show severe self-mutilations with the absence of a toe. (B) Manhattan plot of -log10 transformed p-values by canine chromosome highlighting a strong signal on chromosome 4 (Wald test; p ≤ 10-16). C) Manhattan plot of -log10 empirical p-values (EMP) obtained by permutations test confirming the signal on chromosome 4 (Permutation = 100,000).

http://dx.doi.org/10.1371/journal.pgen.1006482.g001

We thus investigated these affected breeds as a potential opportunity to identify new genes for unexplained cases of human HSANs. We first started the genetic study with the French Spaniel breed and identified one locus strongly associated with the disease, on canine chromosome 4, and refined the locus by intersecting the haplotypes of the other affected breeds. Using this strategy, we highlighted a common homozygous haplotype of 1.8 Mb (chr4:70,6-72,4Mb) shared by all affected dogs in the four sporting breeds. Next, we performed targeted DNA sequencing of the locus in four affected and four unaffected dogs. We found 478 variants which were individually removed by Allele-Specific PCR except one located upstream of the GDNF gene (Glial cell-Derived Neurotrophic Factor), a neurotrophic factor involved in neuronal development and adult neuronal survival [22-25]. More precisely, we showed that this variant is localized in the last exon of a canine long-non-coding RNA (GDNF-AS) that is transcribed in the opposite direction of GDNF. Comparative genomic analysis with the human orthologous genomic region led us to hypothesize that the mutation disrupts a regulatory region, such as an enhancer, controlling GDNF expression. We then performed quantitative PCR analyses and showed lower expression levels of both GDNF and GDNF-AS RNAs in affected dogs compared to unaffected dogs. By gel shift assays, we further showed that the mutation modified the affinity of a regulatory complex. Such findings highlight the GDNF/GDNF-AS partnership as a subtle regulatory mechanism and relevant candidates to search for human HSAN mutations and potentially to develop targeted therapies.

Results

Genome-wide association study and GDNF sequencing

We first focused the genetic analysis on the French Spaniel breed. The collection of 173 blood samples with pedigree information allowed us to design a large pedigree and to confirm the monogenic recessive inheritance of the disease in this breed as previously described [21]. Affected dogs presented with severe self-mutilations with insensitivity to pain only in the feet (Fig 1A). Neurological examinations (proprioception, motor abilities, spinal reflex exams) confirmed that pain perception became progressively normal above the knees and was not altered in the rest of the body. One affected dog had never shown self-mutilation marks leading us to consider insensitivity to pain as the first clinical sign for this disease. As described in Paradis et al, no signs of autonomic disturbances were detected by neurological exams [21] and motor abilities and proprioception were normal in all dogs.

We first performed a genome-wide association study (GWAS) using a linear mixed-model method accounting for population stratification and relatedness, with 49 French spaniels (21 affected showing insensitivity to pain in their feet and 28 unaffected dogs) (Fig 1B). Following quality control of the genotyping data, 123,579 SNPs were retained for genetic mapping. This analysis identified a 3 Mb locus on canine chromosome 4 strongly associated with the disease, with significant p-values (Wald test; p-value = 10-16; Table 1) and confirmed the association using a permutation test (p100,000 = 10-5) (Fig 1C). Since all affected French spaniels present a homozygous haplotype in this 3 Mb locus, we took advantage of the particular genetic structure of the dog population. Based on the hypothesis that sporting breeds share the same disease and founder mutation, we genotyped additional affected dogs in three sporting dog breeds: German Short-haired Pointer (GP), English Pointer (EP) and English Springer Spaniel (ESS). These breeds revealed a common homozygous haplotype block with French Spaniel (FS) confirming the founder effect between these four breeds. Recombination breakpoints in the FS and ESS reduced the critical interval to 1.8 Mb (chr4:70,6-72,4Mb) (Fig 2). This interval contains a dozen genes annotated in the dog and human orthologous genomic region including the GDNF gene (Glial cell-Derived Neurotrophic Factor), which appeared as a new and good candidate for HSAN (Table 1). Indeed, GDNF codes a neurotrophic factor involved in neuronal development and adult neuronal survival [22-25] and was never associated with HSAN. However, Sanger sequencing did not reveal any mutation neither in the three exons of the transcripts, nor in the promoter sequence of GDNF.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 2. UCSC screenshot of the dog chromosome 4 locus (CanFam3) found by GWAS and its orthologous genomic region in human (GRCh38).

The upper panel represents the canine locus with the first track showing homozygous haplotype blocks found in four sporting breeds (shaded greys) and the bottom horizontal bar corresponding to the common 1.8 Mb locus shared by the four breeds (black). Next tracks show the 3 Mb region targeted for high throughput sequencing (in red) and the resulting 478 variants identified by GATK program (dark green). The AMS variant is represented as a vertical bar (light green) and highlighted by a green arrow. Finally, the non-dog RefSeq genes (longest isoform) are represented in blue. The bottom panel corresponds to the human orthologous region (inverted orientation) of the canine targeted interval as given by the LiftOver tools available on the UCSC genome browser [26,27]. Human RefSeq genes are annotated in blue and the predicted orthologous genomic position of the AMS variant corresponds to the green arrow.

http://dx.doi.org/10.1371/journal.pgen.1006482.g002

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Table 1. Best-SNPs obtained by GWAS using 21 affected vs 28 unaffected French spaniels.

http://dx.doi.org/10.1371/journal.pgen.1006482.t001

Targeted re-sequencing of the associated locus and validation of the mutation

To identify the disease-associated variants (InDel and SNP), we chose to enlarge the homozygous candidate region to 3 Mb (70.3 Mb to 73.3 Mb on CanFam3 assembly) and performed a targeted sequencing (Fig 2). Four affected and four unaffected sporting dogs were thus sequenced with high coverage and depth (coverage 25X: 99%; 10 X: 100%; mean depth: 374X; S1 Table). Within this 3 Mb interval, 478 variants were identified with the expected distribution between affected and unaffected individuals (Tables 2 and 3). After QC and filtering for known SNPs from the dbSNP database (dbSNP build 139), 156 variants still remained (Table 4 and S2 Table). To exclude non-disease variants due to polymorphisms, we then chose to genotype dogs from a large number of unaffected breeds making the hypothesis that the rare allele present in affected dogs should not be found in unaffected breeds (S3 Table). Using AS-PCR experiment and Sanger sequencing, we found only one SNP perfectly segregating with the disease and located in an intergenic region ~90 kb upstream of the GDNF gene (chr4.g.70,875,561C>T). To validate this variant, we sequenced over 900 unaffected dogs from 130 different breeds and observed that the mutated allele was never detected in these dogs, thus excluding a polymorphism. We also checked the allelic distribution of this variant in a panel of 250 sporting dogs of known clinical status and we observed a perfect correlation between the clinical status of dogs and their genetic status, confirming the strong association of the variant with the disease in the analyzed sporting breeds. Interestingly, heterozygous dogs did not present any clinical signs and only one dog was homozygous mutated with only insensitivity to pain (no acral mutilations). Even if it was lying in a not yet described genomic region, we considered this variant as potentially causal and further named it "AMS variant" for Acral Mutilation Syndrome variant.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Table 2. Summary of the SNPs found in the targeted sequencing.

Only variants that are homozygous for a rare allele in all affected dogs and not homozygous in controls were considered to perfectly segregate. Nb SNPs: Number of SNPs; Hom: Nb Homozygous variants; HTZ: Nb Heterozygous variants.

http://dx.doi.org/10.1371/journal.pgen.1006482.t002

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Table 3. Summary of the indels found in the targeted sequencing.

Only variants that are homozygous for a rare allele in all affected dogs and not homozygous in controls were considered to perfectly segregate. Hom: Nb Homozygous variants; HTZ: Nb Heterozygous variants.

http://dx.doi.org/10.1371/journal.pgen.1006482.t003

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Table 4. Summary of the variants obtained in the targeted sequencing before and after filtering.

http://dx.doi.org/10.1371/journal.pgen.1006482.t004

Genomic analyses of the chromosome 4 locus

To unravel the genomic context surrounding the GDNF locus, we mapped the canine 3 Mb locus targeted for sequencing (chr4:70,3-73,3 Mb) onto the human genome (GRCh38) using the LiftOver tools available on the UCSC genome browser [26,27] and the AutoGRAPH website [28] (Fig 2). The locus is well conserved between dog and human with 63.3% sequence identity with human region (chr5:34,9 Mb-38,5 Mb) and 53.1% reciprocally. At the syntonic level, most of the canine orthologous protein-coding genes harbor the same order but in the opposite orientation of transcription as compared to humans. Remarkably in human, a long non-coding RNA gene (GDNF-AS) is annotated upstream of GDNF [29], transcribed in a divergent orientation, and close to the orthologous position of the canine AMS variant (Fig 2).

To refine the annotation of the canine locus, we used 33 RNA sequencing data produced from 25 different tissues of seven breeds (S4 Table). The sequencing data were provided by the Broad Institute and as part of a collaborative work of the LUPA consortium [30-32]. This extended annotation confirmed the presence of a multi-exonic RNA transcribed in the vicinity and in a divergent orientation of GDNF as in human (Fig 3). Next, to assess the protein-coding capabilities of each of the 12 isoforms of the new gene, we used three complementary programs (CPC [33], CPAT [34], PLEK [35]) and confirmed that all of the new transcripts are classified as probable non-coding RNAs (See Methods and S1 File). Remarkably, the AMS variant is located in the last exon of the canine lincRNA, which is shared by all of the 12 isoforms (Fig 3). The lincRNA, that we named canine GDNF-AS, is weakly expressed in the 33 RNA-seq (S5 Table) and is mostly detected in the olfactory bulb (RPKM = 0.849) and spleen (RPKM = 0.732). We also observed that the highest level of GDNF mRNA expression level is in spleen (RPKM = 4.6).

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 3. UCSC screenshots of the candidate GDNF/GDNF-AS regions in dog (CanFam3 assembly) and in human (GRCh38 assembly).

The upper panel represents the canine GDNF/GDNF-AS locus (chr4:70,8-71 Mb) with a subset of the variants identified by targeted sequencing in black and the AMS variant highlighted in green, the GDNF mRNA isoforms annotated by the BROAD institute [32] (in light blue), the mapping of the 12 canine GDNF-AS isoforms identified in our RNA-seq analysis (in black), and the RepeatMasker track of SINEs and LINEs in the locus (shaded gray). The bottom panel represents the human orthologous locus (chr5:37,8-37,9 Mb) with tracks corresponding (from the top to bottom) to the alignment of the canine GDNF-AS isoform (TCONS_00182466) using BLAT [26] (in grey and red) with the orthologous position of the AMS variant (in green), the RefSeq genes (in blue) and the H3K4me1 and H3K27ac chromatin signals annotated by the ENCODE project [36] with NHEK cell line in purple. The dotted arrows illustrate the inversion of the syntonic block between dog and human species.

http://dx.doi.org/10.1371/journal.pgen.1006482.g003

It has been shown that lincRNAs transcribed from a bi-directional promoter could regulate the level of expression of the corresponding neighbor mRNAs [37-39]. We therefore computed the correlation of expression using the 33 expression data values (RPKM) between the canine GDNF-AS and the four closest protein-coding genes i.e. GDNF, WDR70, NUP155 and C5orf42 (S1 Fig). We found that GDNF-AS expression is highly correlated with GDNF (Pearson correlation = 0.67, p-value = 1.8 10-5) but not with the other protein-coding genes (Pearson correlations = -0.10, -0.03 and 0.08 for WDR70, NUP155 and C5orf42, respectively). These results, together with the close proximity of GDNF and GDNF-AS transcription start sites (TSSs), reinforced the relationship between both genes and pointed toward the presence of a bi-directional promoter [37,40,41].

We next analyzed the sequence conservation of the dog genomic region encompassing GDNF/GDNF-AS (chr4:70,8-71 Mb) and showed that it matches the human genome (chr5:37,8-37,9 Mb) with 84.2% of nucleic acid identity (Fig 3). Interestingly, the sequence surrounding the AMS variant on the human genome overlaps with regulatory regions defined by ENCODE data [36,42] using typical histone marks such as monomethylation of histone H3 at lysine 4 (H3K4me1) and acetylation of histone H3 at lysine 27 (H3K27ac) (Fig 3). These typical histone marks led us to hypothesize [43] that the AMS variant is located in a region which could control GDNF and GDNF-AS expression.

Functional analyses

To confirm relationship between GDNF and GDNF-AS, we performed quantitative real-time PCR analyses on Dorsal Root Ganglia (DRG). Indeed, DRG are interesting, as it was previously recognized as harboring the predominant lesions in these breeds [18,19] and contain sensory neurons where neurotrophic factors are highly expressed [22,23,44-47]. We thus extracted RNAs from DRG of two affected and two unaffected dogs and characterized the expression levels of GDNF and GDNF-AS together with the neighboring genes WDR70, NUP155, and SLC1A3. QRT-PCR revealed that the GDNF mRNA level was significantly decreased (60%) in affected dogs (t-test: p-value *<0.0001). We also showed that the expression levels of the other genes of the locus were not altered (Fig 4). Finally, we also observed a strong decrease of the GDNF-AS expression level (80%) in this tissue in agreement with the hypothesis that GDNF-AS and GDNF share a bi-directional promoter highlighted by correlated expression shown by RNA-seq and genomic analyses (Fig 3, S1 Fig and S5 Table).

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 4. qRT-PCR results in Dorsal Root Ganglia using the [delta][delta]Ct methodology (cases in grey, controls in white).

Measures obtained from two DRG by individual were pooled, and we also pooled affected dogs and controls. We used PPIB (Peptidylprolyl Isomerase B) as control and all genes of the locus (GDNF-AS, GDNF, WDR70, NUP155, SLC1A3) were tested 3 times with 3 replicates for each (t-test: p-value *<0.0001).

http://dx.doi.org/10.1371/journal.pgen.1006482.g004

The comparative genomic analysis indicates that the region containing the AMS mutation corresponds to an enhancer element (Fig 3). Mutations in regulatory elements are already described in cancer [48,49] or in developmental disorders such as Van der Woude syndrome [50] or isolated pancreatic agenesis [51]. To test the hypothesis of a mutated regulatory element, we constructed reporter systems using a pTAL-Luc plasmid in which the wild-type and mutated sequences were linked to a firefly luciferase reporter gene with a promoter (S2 Fig). While we detected a weak luciferase signal (p-value < 0.05) showing potential enhancer activity with these clones in HeLa cells, no significant difference between cases and controls were observed. However, the distance between the AMS variant and the promoter is ~2kb in these constructs which is a short distance compared to the ~90kb of the genomic region. We thus hypothesized that a long-distance regulation could exist, explaining the qRT-PCR results. The AMS variant region should thus contain regulatory elements sequences.

To test this hypothesis, we performed Electrophoretic Mobility Shift Assay (EMSA). We first used increasing amounts of nuclear extract from HeLa cells with duplex of wild-type or mutated sequences. We showed a high-affinity specific interaction of a Nuclear Extract (NE) complex with the double-stranded oligonucleotide containing the wild-type sequence but a low-affinity with the sequence containing the AMS variant (Fig 5A). In addition, the competition experiment showed that the complex binds preferentially on the wild-type sequence rather than on the mutated sequence. This result indicates that the affinity of the complexes formed on the two sequences is different, (Fig 5B). The same results were obtained when using human neuroblastoma cell line NE (SY5Y). These data indicate that the mutation prevents the binding of a nuclear complex supporting the hypothesis that the decrease expression of both GDNF and GDNF-AS observed in qRT-PCR experiments (Fig 4) is due to a lack of activity of the mutated enhancer.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 5. The AMS variant affect binding of nuclear complex.

EMSAs were performed with wild type (WT) or mutated (Mut) duplex and nuclear extract (NE). (A) Detection of a mobility shift after incubation of 40 fmoles of radiolabelled WT-duplex (lane 1 to 4) or Mut-duplex (lane 5 to 8) with increasing amounts (0,3-0,6-1,2-2,4 Ml) of HeLa NE. (B) 40 fmoles of radiolabelled WT-duplex were incubated without (lane 9 and 13) or with 1Ml of HeLa NE (lane 10 to 12) or 5 Ml of SY5Y (lane 14 to 16) NE and in the presence of 2 pmol of competitor WT-duplex (lane 11 and 15) or Mut-duplex (lane 12 and 16). Black bar: specific complex associated with the radioactive duplex, asterisk: non specific complex.

http://dx.doi.org/10.1371/journal.pgen.1006482.g005

Discussion

This study presents the identification of the mutation involved in a severe form of hereditary sensory autosomal neuropathy (HSAN) shared by four sporting breeds having a common history. The identification of a 1.8 Mb homozygous region on chromosome 4 by GWAS, with only one variant perfectly segregating with the disease, illustrates the advantages of the canine model in genetics. Since this variant is located 90 kb upstream of the neurotrophic factor GDNF, we hypothesized a regulatory effect of this variant on GDNF. We then showed a strong correlation between the presence of the mutation and the decreased expression levels of both GDNF and its divergent long non-coding RNA GDNF-AS. Finally, we demonstrated the impact of the mutation on the binding of a complex confirming the hypothesis of a regulatory element. With this work, we provide relevant arguments to explore GDNF/GDNF-AS partnership in human patients with unexplained HSAN.

We started this project focusing our study on dogs with self-mutilations using a precise clinical questionnaire and we quickly detected that all dogs with self-mutilations also presented insensitivity to pain, not always reported. This important clinical sign led us to improve the genetic analysis. Indeed, we found additional affected dogs without self-mutilations but related to cases with self-mutilations. These dogs showed the "affected" homozygous haplotype, while owners did not detect the insensitivity to pain. This observation reflects the difficulty to diagnose insensitivity to pain in dogs, which contributed to the spreading of this severe disorder in the four related breeds [18,19,21]. This feature led us to consider self-mutilation, probably triggered by small fractures of toes or other injuries [21], as a consequence of the insensitivity to pain.

With a high coverage targeted sequencing of the locus (365x), we identified one variant perfectly associated with the phenotype, and absent in more than 900 unaffected dogs. This variant located upstream of the GDNF gene appeared as an excellent candidate. Indeed, GDNF (glial cell-derived neurotrophic factor) promotes the axon development and the survival of sensory neurons at different stages of the development in mice and chicken model [22-25]. Moreover, it has been shown in vitro that after injury, nearly 100% of the sensory neurons are rescued by the production or injection of GDNF [24,46,52,53]. We thus expected that this AMS variant could disturb the expression of GDNF leading to the decrease of the number of sensory neurons and then their death, as previously observed in affected dogs [18-20]. This feature could happen during the first months of life, which would correlate well with the early age of onset of this insensitivity to pain.

The genomic analyses of the human and dog orthologous regions revealed the presence of long non-coding RNA transcripts in both species. One of them, GDNF-AS, described in human [29], contains the same structure as the canine lincRNAs identified in this study and more likely share a bi-directional promoter with the GDNF gene. While the GDNF mRNA sequence is well conserved between human, mouse and dog, the canine GDNF-AS transcript sequences diverge significantly from the human sequences at the nucleotide level (only 14.3% identity using pairwise alignment with longest isoforms) and no GDNF-AS are yet annotated in mouse. However, the whole locus is well conserved at the syntonic and nucleotide levels between human and dog. Regarding shared lncRNAs functions between species in spite of low sequence conservation, several well-described lncRNAs (e.g. Xist, TUNA and Hotair...) exhibit conserved biological functions even though nucleotide conservation is only limited to small patches and their exon/intron structure is not maintained [54].

An interesting point is the fact that the AMS variant is located in the last exon of GDNF-AS, which is the only exon shared by all of the 12 isoforms annotated by our RNA-seq data. Strikingly, another neurotrophic factor, BDNF (Brain Derived Neurotrophic Factor precursor), playing similar functions as GDNF, also presents a similar genomic feature than GDNF, with the presence of a lincRNA transcribed in the opposite direction. The parallel between BDNF and GDNF is interesting because they both stimulate the growth and differentiation of new neurons and support the survival of existing neurons in central and peripheral nervous systems [23,52,55]. Moreover, BDNF interacts with the TrkB receptor encoded by the NTRK2 gene, which has a paralog NTRK1, already annotated as mutated in HSAN type IV [5,56,57]. In addition, it has been shown that the inhibition of BDNF-AS increases BDNF levels in vivo [38,58]. Modarresi et al. recently showed that the knock-out of GDNF-AS increases the expression level of GDNF in human HEK293T cells [59]. We could thus hypothesize that the GDNF and GDNF-AS partnership could fine-tune the expression of GDNF and that the AMS variant disturbs this relationship at a specific developmental stage [60].

In embryonic rats, dorsal root ganglia (DRG) neurons were described as potential target of GDNF [61,62]. In addition, a recent paper in human and rat confirmed that DRG are the main tissue where genes involved in sensory neuropathies are highly expressed [47]. Thus, the qRT-PCR analyses in canine DRG brought a relevant result: the AMS variant seems to affect a regulatory element. Indeed, while the AMS variant is located in the last exon of GDNF-AS, both GDNF and GDNF-AS transcripts are weakly expressed in affected dogs as compared to unaffected dogs (Fig 4). This important observation supports the implication of GDNF in the disease, and suggests that too low levels of GDNF in the peripheral nervous system prevent the maintenance of the integrity of sensory adult neurons [22-25]. Moreover, this observation led us to consider the hypothesis of a mutated regulatory element that could impair the activation of the promoter shared by GDNF and its lincRNA GDNF-AS.

While we detected an enhancer activity, we did not find any difference of activity between the wild-type and the mutated sequences in the luciferase reporter assay. This is probably due to the fact that GDNF regulation is modulated by a long distance interaction between the region containing the AMS variant and the promoter. As the luciferase reporter assay did not validate this hypothesis, we checked if the AMS variant could inhibit the enhancer activity by preventing the binding of a transcription complex. The EMSA experiment clearly showed that a complex binding the wild-type sequence has less affinity for the mutated sequence. The fact that the mutation of a single nucleotide, in a 20-nucleotides sequence, inhibits the fixation of regulatory elements comforts the causal effect of the variant. It will be interesting to further investigate the nature of the complex as we noticed that the AMS variant introduces a change of one nucleotide in the conserved binding motif recognized by NEUROD1 (S3 Fig). This transcription factor plays an important role during the neurogenesis by binding to regulatory elements of neuronal genes that are developmentally silenced by epigenetic mechanisms [65]. In this context, decreased binding of neurogenic transcription factors to the mutated enhancer in dog could explain the weak expression of GDNF/GDNF-AS partnership in neuropathies.

While the complete knock-out GDNF-/- mouse model is not viable since they die at 1-1.5 days after birth [63,64], the AMS variant is not lethal in dogs. In affected dogs, we showed that GDNF and GDNF-AS RNAs are still produced but at very low levels. With the functional analyses, we can conclude that the AMS variant affects the regulation of GDNF expression, which could explain the non-lethal effect of the mutation in homozygous dogs who only present insensitivity to pain. However, heterozygous dogs are not affected similar to GDNF+/- mice. Indeed GDNF+/- mice were normal, viable and indistinguishable within controls by visual inspection [64]. Comparison with the mouse model revealed another difference: in dogs, no sign of autonomic disturbances was detected by neurological exams [21] and no metanephric kidney or gastrointestinal tract disorders were observed contrary to what is reported in GDNF-/- (knock-out) mice [63,64]. These results confirm that GDNF acts in a dose-dependent manner as suggested by GDNF-/- mice model [63].

In conclusion, considering the parallel with BDNF function, our RNA expression analyses, and dog/human comparative genomic analyses, we hypothesize that the GDNF/GDNF-AS partnership represents a very precise regulation mechanism for the mRNA expression of GDNF. Indeed, since we showed that the mutation (in the last exon of GDNF-AS) dramatically diminishes both the expression of GDNF/GDNF-AS, by probably disrupting a regulatory element, the AMS variant most probably leads to a decrease of the GDNF expression, preventing a correct growth and maintenance of the small fibers from the knees down to the paws. To date, no mutation in the GDNF gene or in its regulatory regions has been described in human HSAN patients. This work illustrates the power of the canine model in genetics and shows the potential implication of GDNF in such neurological diseases, revealing the importance to explore this gene and its regulatory elements as serious candidates in human HSAN patients.

Methods

Ethics statement and sample collection

All dogs were client-owned and no harmful procedures were performed, so there was no animal experimentation according to the European legal definition (Subject 5b and 5f of Article1, Chapter I of the Directive 2010/63/UE of the European Parliament and of the Council). Blood and biopsies were obtained as part of routine clinical procedures for diagnostic purposes, approved by the CNRS ethical board (France), and were sent by veterinarians and French Veterinary Schools. The biological samples were obtained from the 'Cani-DNA_CRB', which is part of the CRB-Anim infrastructure, ANR-11-INBS-0003 (http://dog-genetics.genouest.org). The work with dog samples was approved by the CNRS ethical board, France (35-238-13) for UMR6290.

Samples

Genomic DNA was extracted from 2 mL blood samples collected in EDTA, using the NucleoSpin Blood L kit (Macherey-Nagel) according to the manufacturer's instructions. Epidemiological and clinical data were collected using a dedicated questionnaire for each affected dog. Dogs over three years old without any symptoms were considered unaffected. Complementary neurological diagnoses (proprioception, motor abilities, spinal reflex exams, electromyography) were done by neurologists, in the Veterinary Schools of Lyon (France) and St-Hyacinthe (Canada).

Genome-wide association study

Using the Illumina Canine HD 173k (BeadChip), genotyping was performed at the Centre National de Génotypage (CNG; Evry, France) in the frame of European LUPA project in 62 dogs: 49 French spaniels (21 affected and 28 unaffected dogs). This design contained 2 subpopulations: one from France and one from Canada, each including 3 small families (parents + 2 sibling including one affected and one unaffected), all other dogs were unrelated to one another at the grandparent level. To reduce the 3 Mb locus, additional dogs were genotyped including 5 German Short-haired Pointers (1 affected dog, his unaffected brother and father, and 2 other unrelated unaffected dogs), 4 unrelated English Pointers (1 affected and 3 unaffected dogs) and 4 unrelated English Springer Spaniels (1 affected, and 3 unaffected dogs). Using the PLINK software (v1.06-1.07) [66] for SNP filtering (minor allele frequency <0.01; call rate by marker and individual>75%), a dataset of 123,579 SNPs was obtained and analyzed. Genome-wide association study was performed using the software GEMMA v0.94.1 (Genome-wide Efficient Mixed-Model Association) [67,68]. The linear mixed-model method accounts for population stratification and relatedness between dogs, especially for small families.

Targeted capture by next-generation sequencing and allele-specific PCR

Genomic library sample preparation was performed using the Illumina paired-end library sample preparation kit (Illumina Inc., San Diego, CA, USA). Sample preparation was carried out by Integragen (Evry, France) according to the manufacturer's instructions (Agilent), using 4Mg of genomic DNA. We used two unrelated affected dogs in French Spaniel (FS), one affected German Shorthaired Pointer (GP) and one affected English Springer Spaniel (S1 Table). For controls, we used father of one the two affected FS + one unrelated, one brother of the affected GP and one unrelated ESS. All these dogs were previously genotyped on the Illumina Canine HD 173k. For targeted sequencing, solution-based capture was performed using the Agilent SureSelect Target Enrichment System Kit. Repeated elements were not sequenced due to the limit of this technology, this is why we extended the size of the targeted locus to 3 Mb (chr4: 70.3-73.3 Mb, CanFam3). A custom panel of 75 bases cRNA primers was designed to sequence the 3 Mb locus based on the canine genome (CanFam3 assembly) in the candidate region on chromosome 4. The targeted regions were covered by approximately 25,000 probes designed for 2x coverage (i.e., each base was covered by two different probes). Targets were pulled down via streptavidin magnetic beads, purified, and enriched through 13 cycles of PCR amplification. Samples were paired-end sequenced on an Illumina HiSeq2000. Reads mapping was performed using Bowtie 2 [69] and variants were called using GATK 3.5 [70]. Candidate variants were then genotyped on 100 breeds, each breed corresponding to a pool of four DNA samples of unrelated dogs (S3 Table). We also used small French Spaniel families previously genotyped, to create 12 different pools used as references to detect the presence of the rare allele. Genotyping was performed using an Allele Specific PCR method (AS-PCR) (Integragen) developed from the methodology described by Nazarenko and Myakishev [71,72]. SNP genotyping was performed on the Biomark (Fluidigm) with a microfluidic multiplex 96.96 dynamic array chip. We excluded all variants where rare alleles were found in minimum three different breeds (minimum one dog by breed).

Identification of the candidate variant in dogs

All primers used were designed using Primer3plus (www.bioinformatics.nl/primer3plus/) [73]. We re-sequenced the GDNF gene and the variant on 16 affected and 16 unaffected French Spaniel using the Type-it Microsatellite PCR kit (Qiagen) on C1000 and S1000 thermocyclers (Bio-Rad). PCR products were cleaned using the Illustra Exostar-1 Step reagent (Dutscher) and sequenced by Sanger method using the BigDye Terminator v3.1 cycle sequencing kit (Thermo Fisher Scientific). Electrophoresis of the products was realized on a 370 ABI sequencer. Sequences were analyzed using the PhredPhrap and Consed software's [74-76]. The presence or absence of the variant was validated on 250 French Spaniel, German-Short-Haired Pointer, English Pointer and English Springer Spaniel (24 affected and 226 unaffected), and 900 dogs (including the 400 dogs used for AS-PCR) from 130 unaffected breeds among the 10 FCI (Fédération Cynologique Internationale) groups.

RNA sequencing analyses and GDNF-AS annotation

To detect the GDNF-AS long non-coding RNA (lincRNA), we analyzed 33 RNA sequencing data produced from 25 different tissues of 7 breeds (S4 Table) as a part of the LUPA consortium and the Broad Institute annotation efforts [31]. Each RNA sequencing data includes between 60 and 150 millions paired-end and strand-specific reads respectively, that were analyzed by Bowtie/Tophat2 and CUFFLINKS/CUFFMERGE v2.2.1 [69,77] revealing 245,276 transcripts belonging to 81,363 genes. The protein-coding capabilities of all novel transcripts (i.e. not annotated in Ensembl) was measured using three complementary programs CPC, CPAT, PLEK [33-35] (S1 File). We quantified RNA expression at the gene-level using the HTSEQ-count program version 0.6.1 [78] where each gene is considered as the union of its exons. This gene-count based measure was then transformed in RPKM using the edgeR software version 3.14.0 [79] in order to normalize gene expression by both its effective length and the total number of mapped reads (thus the sequencing depth) in each of the 33 samples [80]. Therefore, each gene is assigned a vector of 33 points corresponding to the normalized expression of the gene in the 33 RNA-seq samples. Using the cor.test function of the R software (https://www.r-project.org), we computed all pairwise Pearson correlations for the 5 representative genes in the AMS locus i.e. GDNF-AS, GDNF, WDR70, NUP155, and C5orf42 (S1 Fig) based on the 33 points.

Candidate gene expression analysis

RNA was extracted from tissues, using the NucleoSpin RNA kit (Macherey-Nagel) according to the manufacturer's instructions. We used two affected French Spaniel (brother and sister) and two unrelated dogs for controls, all already genotyped on the Illumina Canine HD 173k. Reverse transcription was performed on 1 Mg of total RNA using the High-capacity cDNA Reverse Transcription kit (Thermo Fisher Scientific), according to the manufacturer's instructions. The total Dorsal Root Ganglia cDNA of the canine GDNF and GDNF-AS were amplified and sequenced using two primer pairs (S6 Table). qRT-PCR was performed on 1:80 diluted cDNA samples after pre-amplification with the TaqMan PCR master mix (Thermo Fisher Scientific) on the 7900HT Fast Real-Time PCR System (Applied Biosystems) using standard procedures. We used pre-designed primers for GDNF (Thermo Fisher Scientific Reference: Cf03986046_g1), WRD70 (TFS Ref: Cf02651565_m1), NUP155 (TFS Ref: Cf02644933_m1) and SLC1A3 (TFS Ref: Cf02702629_m1) and we also specifically designed other probes for GDNF-AS (S6 Table). Canine PPIB (Peptidylprolyl Isomerase B) was used as the reference gene (Thermo Fisher Scientific Reference: Cf02629556_m1). Each sample was measured in triplicate and each qRT-PCR was carried out three times by different experimenters. Relative amounts of the transcript were determined using the [delta][delta]Ct method [81]. Expression analyses were carried out on two tissues (DRG L6, DRG L7) and their results were pooled. We also pooled results from both affected dogs and both controls to increase the number of points. Using the t.test function of the R software (https://www.r-project.org), we determined the significance of the variation of expression level for all genes.

Cloning and luciferase assays

Clones were created to be centered on the AMS variant and correspond to a fragment of 2359 bp containing two SINE, two LINE, one LTR, and eight SNVs. PCR was performed from the genomic DNA of two dogs: one affected dog carrying the mutation and rare alleles for 8 SNVs, and one unaffected dog with the same haplotype as the dog genome reference. The PCR fragments were cloned into the pTAL Luciferase plasmid (Clontech Laboratories, Inc.) by DNA ligation or by homologous recombination using the Gibson assembly MasterMix (New England, Biolabs). All primer sequences used to prepare these constructs are given in S6 Table. All constructs were verified by DNA sequencing. Hela cells were grown in high-glucose Dulbecco's Modified Eagle's Medium (DMEM) containing 10% fetal calf serum (GIBCO) according to the manufacturer's instructions. Cells were transfected in triplicates with 250 ng of reporter plasmid and 50 ng of Renilla control vector (pRL-Null from Promega) using JetPEI (Polyplus Transfection) in 24-well plates. After 24 h, Firefly and Renilla luciferase activities were determined using the Dual-Luciferase Reporter Assay System (Promega) on a Veritas Microplate Luminometer (Turner Biosystems). Enhancer assays of selected regulatory regions were run as previously described [82]. Firefly luciferase activities of individual transfections were normalized against Renilla luciferase activities. As positive control (control +), we used a home-made construction with enhancer activity previously described by Sérandour et al [82].

Cell culture and nuclear extracts

HeLa nuclear extract were purchased from the Computer Cell Culture Center S.A. (Belgium) and prepared according to Dignam et al [83]. SY5Y cells were grown in 50% Dulbecco's Modified Eagle Medium (DMEM) and 50% Nutrient Mixture F-12 with 10% of fetal calf serum and 1% of antibiotics and maintained at 37°C and 5% CO2. Cells were harvested using 1 ml of Lysis buffer containing protease inhibitors, 10 mM Tris-HCl (pH 7.5), 0,5% NP40, 2 mM MgCl2, 3 mM CaCl2 and 10% Glycerol and gently spun down. The supernatant was removed, and the cells were resuspended in the lysis buffer and gently spun down. The supernatant was removed and the cells were resuspended in 500ul of Extraction buffer containing protease inhibitors, 20 mM Hepes, 0,4 M KCL, 1,5 mM MgCl2, 2 mM DTT, 0,5 mM PMSF and 20% Glycerol. The cells were left on ice for 30 min to swell and were then vortexed vigorously for 10 s once during the time on ice. The cells were centrifuged at 12,000 rpm, 15 min at 4°C. The supernatant was kept and the protein concentration was subsequently estimated by the BCA assay.

Electrophoretic Mobility Shift Assays

Probe labeling and duplex formation: 10 pmoles of primers wild type sequence (WT: 5'- TGTTGTCTTTGCTGCTGTCATGATGG-3') and mutated sequence (Mut: 5'- TGTTGTCTTTGCTACTGTCATGATGG-3') were 5'-end labeled using T4-PNK (Promega M4101) and 20MCi of 32P-g ATP. Labeled oligonucleotides were purified on a Sephadex G-25 column. Unlabeled complementary oligonucleotides were then annealed to form a WT-duplex or Mut-duplex in hybridization buffer (10mM Tris pH7.4; 6mM MgCl2; 50mM NaCl; 6mM B-Mercapto-Ethanol) by incubation for 1 min at 95°C and gradually cool down to room temperature. Complete annealing of the radioactive probe was checked by electrophoresis on 20% native acrylamide gel. For complex formation, 40fmoles (15000cpm) of duplex was incubated with various amount of nuclear extracts in a total volume of 10Ml containing 25 mM HEPES (pH 7.5), 150 mM KCl, 5 mM dithiothreitol (DTT), 10% glycerol, 1 Mg of poly (dI-dC). Incubation was performed for 30 min on ice prior to loading. For competition experiments, the mixture was preincubated for 30 min at 4°C with 2 pmoles of unradiolabeled competitor duplex, before addition of the radiolabeled duplex. Complexes were then fractionated onto a 5% polyacrylamide (38:2) gel with 2.5% (v/v) glycerol, 0.5 mM EDTA and 22.25 mM Tris-borate (pH 8.3) by electrophoresis at 4°C in 0.5 mM EDTA, 0.25%glycerol and 22.25 mM Tris-borate (pH 8.3). Free and bound duplexes were visualized with a phosphoimager using the ImageJ software.

Supporting Information

S1 File. Fasta files of the canine GDNF-AS sequences of the 12 isoforms.

(TXT)

S1 Fig. Heatmap of the Pearson correlations between GDNF-AS, GDNF, WDR70, NUP155 and C5orf42 expressions.

(PDF)

S2 Fig. Functional analyses of the candidate regulatory element.

(PDF)

S3 Fig. Canine sequence with a conserved NEUROD1 motif binding site

(PDF)

S1 Table. Description of the targeted sequencing using Agilent technologies.

(XLSX)

S2 Table. List of the 156 variants detected by targeting sequencing.

(XLSX)

S3 Table. Lists of the different breeds used to check candidate variants.

(XLSX)

S4 Table. Lists of the different tissues used for the RNA sequencing.

(XLSX)

S5 Table. Expression values (RPKM) for the Acral Mutilation Syndrome locus and candidate receptors.

(XLSX)

S6 Table. Primers used for Sanger sequencing, qRT-PCR and cloning.

(XLSX)

Acknowledgments

The authors are grateful to referring veterinarians and to all dog owners and breeders who donated samples, pedigree data and follow-up of their dogs. We particularly thank Dr.'s G. Mogicato (ENV Toulouse), G. Marignac and L. Tiret (ENV Alfort), C. Escriou and D. Pin (ENV Lyon), P. Prelaud (Clinique Advetia, Paris), A. Muller (Clinique St Bernard, Lomme), J. Fontaine (Vet School of Liege) and F. Degorce (LAPVSO, Toulouse) for getting us samples as well as C. Fuster, G. Tournier and O. Blot for pedigree and health data. We thank A.S. Lequarré and M.L Ramos-Pamplona, as coordinators for the European FP7 LUPA project, A. Boland and D. Zelenika for the genotyping performed at CNG, Evry, France. We thank the GenOuest Bioinformatics core facility (www.genouest.org) for storing sequencing data, for hosting the Cani-DNA website and for the use of the bioinformatic cluster to analyze the data. We specially acknowledge Emmanuel Martin and Integragen for their expertise in targeted sequencing and AS-PCR experiments, C. Demeautis and G. Hatte for their advices in cloning experiments, Maud Rimbault for the help on the shared haplotype analysis and Luc Paillard for useful discussions and suggestions about EMSA experiments.

Author Contributions

1. Conceptualization: JP CA.

2. Data curation: JP CH TD CA.

3. Formal analysis: JP LL SC BH NB AM PQ TD CA.

4. Funding acquisition: MP CA.

5. Investigation: JP LL SC MP EG BH AP NB MCC PP MB CA RJ CH GS AM PQ TD CA.

6. Methodology: JP LL SC NB AM PQ TD CA.

7. Project administration: JP CA.

8. Resources: MP EG NB MCC PP DS MB MS CA RJ CH GS AM CA.

9. Supervision: CA.

10. Validation: JP LL SC BH NB CH AM PQ TD CA.

11. Visualization: JP LL SC NB AM PQ TD CA.

12. Writing - original draft: JP LL SC MP EG MS RJ CH GS PQ TD CA.

13. Writing - review & editing: JP LL SC NB AM PQ TD CA.

References

Jerath NU, Shy ME. Hereditary motor and sensory neuropathies: Understanding molecular pathogenesis could lead to future treatment strategies. Biochim Biophys Acta. 2015;1852: 667-678. doi: 10.1016/j.bbadis.2014.07.031. pmid:25108281

Rotthier A, Baets J, Timmerman V, Janssens K. Mechanisms of disease in hereditary sensory and autonomic neuropathies. Nat Rev Neurol. Nature Publishing Group; 2012;8: 73-85. doi: 10.1038/nrneurol.2011.227. pmid:22270030

Axelrod FB, Gold-von Simson G. Hereditary sensory and autonomic neuropathies: types II, III, and IV. Orphanet Journal of Rare Diseases. BioMed Central Ltd; 2007;2: 39. doi: 10.1186/1750-1172-2-39. pmid:17915006

Auer-Grumbach M. Hereditary sensory and autonomic neuropathies. Handb Clin Neurol. 2013;115: 893-906. doi: 10.1016/B978-0-444-52902-2.00050-3. pmid:23931820

Indo Y. Neurobiology of pain, interoception and emotional response: lessons from nerve growth factor-dependent neurons. Eur J Neurosci. 2014;39: 375-391. doi: 10.1111/ejn.12448. pmid:24494678

Auer-Grumbach M. Hereditary sensory neuropathy type I. Orphanet Journal of Rare Diseases. BioMed Central Ltd; 2008;3: 7. doi: 10.1186/1750-1172-3-7. pmid:18348718

Murphy SM, Ernst D, Wei Y, Laurà M, Liu Y-T, Polke J, et al. Hereditary sensory and autonomic neuropathy type 1 (HSANI) caused by a novel mutation in SPTLC2. Neurology. Lippincott Williams & Wilkins; 2013;80: 2106-2111. doi: 10.1212/WNL.0b013e318295d789. pmid:23658386

Fridman V, Bundy B, Reilly MM, Pareyson D, Bacon C, Burns J, et al. CMT subtypes and disease burden in patients enrolled in the Inherited Neuropathies Consortium natural history study: a cross-sectional analysis. J Neurol Neurosurg Psychiatr. BMJ Publishing Group Ltd; 2015;86: 873-878. doi: 10.1136/jnnp-2014-308826. pmid:25430934

Granger N. Canine inherited motor and sensory neuropathies: an updated classification in 22 breeds and comparison to Charcot-Marie-Tooth disease. Vet J. 2011;188: 274-285. doi: 10.1016/j.tvjl.2010.06.003. pmid:20638305

Sutter NB, Ostrander EA. Dog star rising: the canine genetic system. Nat Rev Genet. Nature Publishing Group; 2004;5: 900-910. doi: 10.1038/nrg1492. pmid:15573122

Lindblad-Toh K, Wade CM, Mikkelsen TS, Karlsson EK, Jaffe DB, Kamal M, et al. Genome sequence, comparative analysis and haplotype structure of the domestic dog. Nature. Nature Publishing Group; 2005;438: 803-819. doi: 10.1038/nature04338. pmid:16341006

Vaysse A, Ratnakumar A, Derrien T, Axelsson E, Rosengren Pielberg G, Sigurdsson S, et al. Identification of genomic regions associated with phenotypic variation between dog breeds using selection mapping. Akey JM, editor. PLoS Genet. Public Library of Science; 2011;7: e1002316. doi: 10.1371/journal.pgen.1002316. pmid:22022279

Grall A, Guaguère E, Planchais S, Grond S, Bourrat E, Hausser I, et al. PNPLA1 mutations cause autosomal recessive congenital ichthyosis in golden retriever dogs and humans. Nat Genet. 2012;44: 140-147. doi: 10.1038/ng.1056. pmid:22246504

Merveille A-C, Davis EE, Becker-Heck A, Legendre M, Amirav I, Bataille G, et al. CCDC39 is required for assembly of inner dynein arms and the dynein regulatory complex and for normal ciliary motility in humans and dogs. Nat Genet. 2011;43: 72-78. doi: 10.1038/ng.726. pmid:21131972

Bardagí M, Montoliu P, Ferrer L, Fondevila D, Pumarola M. Acral mutilation syndrome in a miniature pinscher. J Comp Pathol. 2011;144: 235-238. doi: 10.1016/j.jcpa.2010.08.014. pmid:20961556

Pivník L. [Comparative problems of acrodystrophic neuropathies in man and dogs]. Schweiz Arch Neurol Neurochir Psychiatr. 1973;112: 365-371. pmid:4725277

Sanda A, Pivník L. Die zehernekrose bei kurghaarigen. Vorstehlunden Kleintierproxis. Vorstehlunden Kleintierproxis; 1964;9: 76

Cummings JF, de Lahunta A, Winn SS. Acral mutilation and nociceptive loss in English pointer dogs. A canine sensory neuropathy. Acta Neuropathol. 1981;53: 119-127. pmid:6259871

Cummings JF, de Lahunta A, Braund KG, Mitchell WJ. Hereditary sensory neuropathy. Nociceptive loss and acral mutilation in pointer dogs: canine hereditary sensory neuropathy. Am J Pathol. American Society for Investigative Pathology; 1983;112: 136-138. pmid:6574711

Cummings JF, de Lahunta A, Simpson ST, McDonald JM. Reduced substance P-like immunoreactivity in hereditary sensory neuropathy of pointer dogs. Acta Neuropathol. 1984;63: 33-40. pmid:6203326

Paradis M, de Jaham C, Pagé N, Sauve F, Helie P. Acral mutilation and analgesia in 13 French spaniels. Vet Dermatol. Blackwell Science Ltd; 2005;16: 87-93. doi: 10.1111/j.1365-3164.2005.00443.x. pmid:15842538

Molliver DC, Wright DE, Leitner ML, Parsadanian AS, Doster K, Wen D, et al. IB4-binding DRG neurons switch from NGF to GDNF dependence in early postnatal life. Neuron. 1997;19: 849-861. pmid:9354331

Buj-Bello A, Buchman VL, Horton A, Rosenthal A, Davies AM. GDNF is an age-specific survival factor for sensory and autonomic neurons. Neuron. 1995;15: 821-828. pmid:7576631

Markus A, Patel TD, Snider WD. Neurotrophic factors and axonal growth. Curr Opin Neurobiol. 2002;12: 523-531. pmid:12367631

Lee C-S, Tee LY, Dusenbery S, Takata T, Golden JP, Pierchala BA, et al. Neurotrophin and GDNF family ligands promote survival and alter excitotoxic vulnerability of neurons derived from murine embryonic stem cells. Exp Neurol. 2005;191: 65-76. doi: 10.1016/j.expneurol.2004.08.025. pmid:15589513

Kent WJ. BLAT--the BLAST-like alignment tool. Genome Res. Cold Spring Harbor Lab; 2002;12: 656-664. doi: 10.1101/gr.229202. pmid:11932250

Speir ML, Zweig AS, Rosenbloom KR, Raney BJ, Paten B, Nejad P, et al. The UCSC Genome Browser database: 2016 update. Nucleic Acids Res. Oxford University Press; 2016;44: D717-25. doi: 10.1093/nar/gkv1275. pmid:26590259

Derrien T, André C, Galibert F, Hitte C. AutoGRAPH: an interactive web server for automating and visualizing comparative genome maps. Bioinformatics. Oxford University Press; 2007;23: 498-499. doi: 10.1093/bioinformatics/btl618. pmid:17145741

Airavaara M, Pletnikova O, Doyle ME, Zhang YE, Troncoso JC, Liu Q-R. Identification of novel GDNF isoforms and cis-antisense GDNFOS gene and their regulation in human middle temporal gyrus of Alzheimer disease. J Biol Chem. 2011;286: 45093-45102. doi: 10.1074/jbc.M111.310250. pmid:22081608

Djebali S, Wucher V, Foissac S, Hitte C, Corre E, Derrien T. Bioinformatics Pipeline for Transcriptome Sequencing Analysis. Methods Mol Biol. New York, NY: Springer New York; 2017;1468: 201-219. doi: 10.1007/978-1-4939-4035-6_14. pmid:27662878

Lequarré A-S, Andersson L, André C, Fredholm M, Hitte C, Leeb T, et al. LUPA: a European initiative taking advantage of the canine genome architecture for unravelling complex disorders in both human and dogs. Vet J. 2011;189: 155-159. doi: 10.1016/j.tvjl.2011.06.013. pmid:21752675

Hoeppner MP, Lundquist A, Pirun M, Meadows JRS, Zamani N, Johnson J, et al. An improved canine genome and a comprehensive catalogue of coding genes and non-coding transcripts. PLoS ONE. 2014;9: e91172. doi: 10.1371/journal.pone.0091172. pmid:24625832

Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. Oxford University Press; 2007;35: W345-9. doi: 10.1093/nar/gkm391. pmid:17631615

Wang L, Park HJ, Dasari S, Wang S, Kocher J-P, Li W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. Oxford University Press; 2013;41: e74-e74. doi: 10.1093/nar/gkt006. pmid:23335781

Li A, Zhang J, Zhou Z. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. BioMed Central Ltd; 2014;15: 311. doi: 10.1186/1471-2105-15-311. pmid:25239089

ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. Nature Research; 2012;489: 57-74. doi: 10.1038/nature11247. pmid:22955616

Uesaka M, Nishimura O, Go Y, Nakashima K, Agata K, Imamura T. Bidirectional promoters are the major source of gene activation-associated non-coding RNAs in mammals. BMC Genomics. 2014;15: 35. doi: 10.1186/1471-2164-15-35. pmid:24438357

Wu P, Zuo X, Deng H, Liu X, Liu L, Ji A. Roles of long noncoding RNAs in brain development, functional diversification and neurodegenerative diseases. Brain Res Bull. 2013;97: 69-80. doi: 10.1016/j.brainresbull.2013.06.001. pmid:23756188

Kim T-K, Hemberg M, Gray JM. Enhancer RNAs: a class of long noncoding RNAs synthesized at enhancers. Cold Spring Harb Perspect Biol. Cold Spring Harbor Lab; 2015;7: a018622. doi: 10.1101/cshperspect.a018622. pmid:25561718

Sigova AA, Mullen AC, Molinie B, Gupta S, Orlando DA, Guenther MG, et al. Divergent transcription of long noncoding RNA/mRNA gene pairs in embryonic stem cells. Proc Natl Acad Sci USA. 2013;110: 2876-2881. doi: 10.1073/pnas.1221904110. pmid:23382218

Lepoivre C, Belhocine M, Bergon A, Griffon A, Yammine M, Vanhille L, et al. Divergent transcription is associated with promoters of transcriptional regulators. BMC Genomics. BioMed Central; 2013;14: 914. doi: 10.1186/1471-2164-14-914. pmid:24365181

Koch CM, Andrews RM, Flicek P, Dillon SC, Karaöz U, Clelland GK, et al. The landscape of histone modifications across 1% of the human genome in five human cell lines. Genome Res. Cold Spring Harbor Lab; 2007;17: 691-707. doi: 10.1101/gr.5704207. pmid:17567990

Pekowska A, Benoukraf T, Zacarias-Cabeza J, Belhocine M, Koch F, Holota H, et al. H3K4 tri-methylation provides an epigenetic signature of active enhancers. EMBO J. EMBO Press; 2011;30: 4198-4210. doi: 10.1038/emboj.2011.295. pmid:21847099

Yamamoto M, Mitsuma N, Ito Y, Hattori N, Nagamatsu M, Li M, et al. Expression of glial cell line-derived neurotrophic factor and GDNFR-alpha mRNAs in human peripheral neuropathies. Brain Res. 1998;809: 175-181. pmid:9853108

Yamamoto M, Sobue G, Yamamoto K, Terao S, Mitsuma T. Expression of mRNAs for neurotrophic factors (NGF, BDNF, NT-3, and GDNF) and their receptors (p75NGFR, trkA, trkB, and trkC) in the adult human peripheral nervous system and nonneural tissues. Neurochem Res. 1996;21: 929-938. pmid:8895847

Matheson CR, Carnahan J, Urich JL, Bocangel D, Zhang TJ, Yan Q. Glial cell line-derived neurotrophic factor (GDNF) is a neurotrophic factor for sensory neurons: comparison with the effects of the neurotrophins. J Neurobiol. 1997;32: 22-32. pmid:8989660

Sapio MR, Goswami SC, Gross JR, Mannes AJ, Iadarola MJ. Transcriptomic analyses of genes and tissues in inherited sensory neuropathies. Exp Neurol. 2016;283: 375-395. doi: 10.1016/j.expneurol.2016.06.023. pmid:27343803

Melton C, Reuter JA, Spacek DV, Snyder M. Recurrent somatic mutations in regulatory regions of human cancer genomes. Nat Genet. Nature Publishing Group; 2015;47: 710-716. doi: 10.1038/ng.3332. pmid:26053494

Katainen R, Dave K, Pitkänen E, Palin K, Kivioja T, Välimäki N, et al. CTCF/cohesin-binding sites are frequently mutated in cancer. Nat Genet. Nature Publishing Group; 2015;47: 818-821. doi: 10.1038/ng.3335. pmid:26053496

Fakhouri WD, Rahimov F, Attanasio C, Kouwenhoven EN, Ferreira De Lima RL, Felix TM, et al. An etiologic regulatory mutation in IRF6 with loss- and gain-of-function effects. Hum Mol Genet. Oxford University Press; 2014;23: 2711-2720. doi: 10.1093/hmg/ddt664. pmid:24442519

Weedon MN, Cebola I, Patch A-M, Flanagan SE, De Franco E, Caswell R, et al. Recessive mutations in a distal PTF1A enhancer cause isolated pancreatic agenesis. Nat Genet. 2014;46: 61-64. doi: 10.1038/ng.2826. pmid:24212882

Boyd JG, Gordon T. Glial cell line-derived neurotrophic factor and brain-derived neurotrophic factor sustain the axonal regeneration of chronically axotomized motoneurons in vivo. Exp Neurol. 2003;183: 610-619. pmid:14552902

Boyd JG, Gordon T. Neurotrophic factors and their receptors in axonal regeneration and functional recovery after peripheral nerve injury. Mol Neurobiol. 2003;27: 277-324. doi: 10.1385/MN:27:3:277. pmid:12845152

Kapusta A, Feschotte C. Volatile evolution of long noncoding RNA repertoires: mechanisms and biological implications. Trends Genet. 2014;30: 439-452. doi: 10.1016/j.tig.2014.08.004. pmid:25218058

Chan JR, Cosgaya JM, Wu YJ, Shooter EM. Neurotrophins are key mediators of the myelination program in the peripheral nervous system. Proc Natl Acad Sci USA. National Acad Sciences; 2001;98: 14661-14668. doi: 10.1073/pnas.251543398. pmid:11717413

Guo W, Ji Y, Wang S, Sun Y, Lu B. Neuronal activity alters BDNF-TrkB signaling kinetics and downstream functions. J Cell Sci. The Company of Biologists Ltd; 2014;127: 2249-2260. doi: 10.1242/jcs.139964. pmid:24634513

Guven Y, Altunoglu U, Aktoren O, Uyguner ZO, Kayserili H, Kaewkahya M, et al. Twins with hereditary sensory and autonomic neuropathy type IV with preserved periodontal sensation. European Journal of Medical Genetics. 2014;57: 240-246. doi: 10.1016/j.ejmg.2014.02.014. pmid:24631696

Wang Y, Zhao X, Ju W, Flory M, Zhong J, Jiang S, et al. Genome-wide differential expression of synaptic long noncoding RNAs in autism spectrum disorder. Transl Psychiatry. Nature Publishing Group; 2015;5: e660. doi: 10.1038/tp.2015.144. pmid:26485544

Modarresi F, Faghihi MA, Lopez-Toledano MA, Fatemi RP, Magistri M, Brothers SP, et al. Inhibition of natural antisense transcripts in vivo results in gene-specific transcriptional upregulation. Nat Biotechnol. Nature Publishing Group; 2012;30: 453-459. doi: 10.1038/nbt.2158. pmid:22446693

Qiu J-J, Ren Z-R, Yan J-B. Identification and functional analysis of long non-coding RNAs in human and mouse early embryos based on single-cell transcriptome data. Oncotarget. Impact Journals; 2016;5

Henderson CE, Phillips HS, Pollock RA, Davies AM, Lemeulle C, Armanini M, et al. GDNF: a potent survival factor for motoneurons present in peripheral nerve and muscle. Science. 1994;266: 1062-1064. pmid:7973664

Trupp M, Rydén M, Jörnvall H, Funakoshi H, Timmusk T, Arenas E, et al. Peripheral expression and biological activities of GDNF, a new neurotrophic factor for avian and mammalian peripheral neurons. J Cell Biol. The Rockefeller University Press; 1995;130: 137-148. pmid:7790368

Pichel JG, Shen L, Sheng HZ, Granholm AC, Drago J, Grinberg A, et al. Defects in enteric innervation and kidney development in mice lacking GDNF. Nature. Nature Publishing Group; 1996;382: 73-76. doi: 10.1038/382073a0. pmid:8657307

Moore MW, Klein RD, Fariñas I, Sauer H, Armanini M, Phillips H, et al. Renal and neuronal abnormalities in mice lacking GDNF. Nature. Nature Publishing Group; 1996;382: 76-79. doi: 10.1038/382076a0. pmid:8657308

Pataskar A, Jung J, Smialowski P, Noack F, Calegari F, Straub T, et al. NeuroD1 reprograms chromatin and transcription factor landscapes to induce the neuronal program. EMBO J. EMBO Press; 2015;: e201591206

Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81: 559-575. doi: 10.1086/519795. pmid:17701901

Zhou X, Stephens M. Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nat Methods. Nature Research; 2014;11: 407-409. doi: 10.1038/nmeth.2848. pmid:24531419

Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. Nature Research; 2012;44: 821-824. doi: 10.1038/ng.2310. pmid:22706312

Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. Nature Publishing Group; 2012;9: 357-359. doi: 10.1038/nmeth.1923. pmid:22388286

McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. Cold Spring Harbor Lab; 2010;20: 1297-1303. doi: 10.1101/gr.107524.110. pmid:20644199

Nazarenko IA, Bhatnagar SK, Hohman RJ. A closed tube format for amplification and detection of DNA based on energy transfer. Nucleic Acids Res. Oxford University Press; 1997;25: 2516-2521. pmid:9171107

Myakishev MV, Khripin Y, Hu S, Hamer DH. High-throughput SNP genotyping by allele-specific PCR with universal energy-transfer-labeled primers. Genome Res. Cold Spring Harbor Laboratory Press; 2001;11: 163-169. pmid:11156625

Untergasser A, Nijveen H, Rao X, Bisseling T, Geurts R, Leunissen JAM. Primer3Plus, an enhanced web interface to Primer3. Nucleic Acids Res. Oxford University Press; 2007;35: W71-4. doi: 10.1093/nar/gkm306. pmid:17485472

Ewing B, Hillier L, Wendl MC, Green P. Base-calling of automated sequencer traces using phred. I. Accuracy assessment. Genome Res. 1998;8: 175-185. pmid:9521921

Ewing B, Green P. Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Res. 1998;8: 186-194. pmid:9521922

Gordon D. Viewing and editing assembled sequences using Consed. Curr Protoc Bioinformatics. Hoboken, NJ, USA: John Wiley & Sons, Inc; 2003;Chapter 11: Unit11.2-11.2.43

Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. Nature Publishing Group; 2012;7: 562-578. doi: 10.1038/nprot.2012.016. pmid:22383036

Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. Oxford University Press; 2015;31: 166-169. doi: 10.1093/bioinformatics/btu638. pmid:25260700

Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. Oxford University Press; 2010;26: 139-140. doi: 10.1093/bioinformatics/btp616. pmid:19910308

Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011;8: 469-477. doi: 10.1038/nmeth.1613. pmid:21623353

Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. Academic Press; 2001;25: 402-408. doi: 10.1006/meth.2001.1262. pmid:11846609

Sérandour AA, Avner S, Percevault F, Demay F, Bizot M, Lucchetti-Miganeh C, et al. Epigenetic switch involved in activation of pioneer factor FOXA1-dependent enhancers. Genome Res. Cold Spring Harbor Lab; 2011;21: 555-565. doi: 10.1101/gr.111534.110. pmid:21233399

Dignam JD, Lebovitz RM, Roeder RG. Accurate transcription initiation by RNA polymerase II in a soluble extract from isolated mammalian nuclei. Nucleic Acids Res. Oxford University Press; 1983;11: 1475-1489. pmid:6828386