**Deletion in the Bardet–Biedl Syndrome Gene** *TTC8* **Results in a Syndromic Retinal Degeneration in Dogs**

**Suvi Mäkeläinen 1, Minas Hellsand 2, Anna Darlene van der Heiden 1, Elina Andersson 3, Elina Thorsson 3, Bodil S. Holst 4, Jens Häggström 4, Ingrid Ljungvall 4, Cathryn Mellersh 5, Finn Hallböök 2, Göran Andersson 1, Björn Ekesten <sup>4</sup> and Tomas F. Bergström 1,\***


Received: 25 August 2020; Accepted: 15 September 2020; Published: 18 September 2020

**Abstract:** In golden retriever dogs, a 1 bp deletion in the canine *TTC8* gene has been shown to cause progressive retinal atrophy (PRA), the canine equivalent of retinitis pigmentosa. In humans, *TTC8* is also implicated in Bardet–Biedl syndrome (BBS). To investigate if the affected dogs only exhibit a non-syndromic PRA or develop a syndromic ciliopathy similar to human BBS, we recruited 10 affected dogs to the study. The progression of PRA for two of the dogs was followed for 2 years, and a rigorous clinical characterization allowed a careful comparison with primary and secondary characteristics of human BBS. In addition to PRA, the dogs showed a spectrum of clinical and morphological signs similar to primary and secondary characteristics of human BBS patients, such as obesity, renal anomalies, sperm defects, and anosmia. We used Oxford Nanopore long-read cDNA sequencing to characterize retinal full-length *TTC8* transcripts in affected and non-affected dogs, the results of which suggest that three isoforms are transcribed in the retina, and the 1 bp deletion is a loss-of-function mutation, resulting in a canine form of Bardet–Biedl syndrome with heterogeneous clinical signs.

**Keywords:** Bardet–Biedl syndrome (BBS); primary cilia; ciliopathy; BBS8; progressive retinal atrophy (PRA); retinitis pigmentosa

#### **1. Introduction**

Inherited retinal degenerations (IRDs) are a diverse group of retinopathies leading to visual impairment and blindness in humans and other species. Most of the IRDs, such as retinitis pigmentosa (RP) in humans (OMIM # 268000) and the canine equivalent, termed progressive retinal atrophy (PRA), are non-syndromic and only affect vision. Syndromic IRDs are less common and besides visual impairment, other organs are also affected. In golden retriever (GR) dogs with PRA, a 1 bp deletion in exon 7 of the Tetratricopeptide repeat domain 8 (*TTC8*) gene was identified in 2014 (CanFam3.1 Chr8:60,090,185delA, rs852355138, OMIA 001984-9615, here denoted as *TTC8*delA) [1]. This form of

PRA is generally referred to as GR-PRA2. The deletion was predicted to cause a frameshift of the open reading frame leading to a premature stop codon in exon 8, 15 codons downstream of the deletion. If translated, the truncated protein would lack most of the tetratricopeptide repeat motifs. In humans, mutations in the *TTC8* gene cause Bardet–Biedl syndrome (BBS; OMIM # 615985), a clinically and genetically heterogeneous autosomal recessive ciliopathy and the second most common human syndromic IRD after Usher syndrome [2,3].

BBS was first described by Georges Bardet and Artur Biedl in the early 1920s [4,5]. The symptoms are highly variable, even between patients from the same family, and are divided into primary and secondary characteristics [6,7]. Primary symptoms include retinal degeneration, obesity, polydactyly, kidney abnormalities, learning disabilities or cognitive impairment, hypogonadism in males, and genital abnormalities in females. Secondary features include speech delay, developmental delay, behavioral abnormalities, eye abnormalities, brachydactyly/syndactyly, ataxia/poor coordination/imbalance, short stature, mild hypertonia, diabetes mellitus, orodental abnormalities, cardiovascular anomalies, situs inversus, hepatic involvement, craniofacial dysmorphism, Hirschsprung disease, and anosmia [2,7,8]. For clinical diagnosis of BBS, it has been suggested that four of the primary characteristics or alternatively three primary and two secondary characteristics should be observed [8].

With 24 genes associated with BBS (OMIM # 209900) to date [9], the syndrome shows large non-allelic heterogeneity [10–13]. The most common cause of the disorder in humans are mutations in the *BBS1* and *BBS10* genes, each accounting for approximately 20% of human cases [14–16]. *TTC8* mutations account for approximately 2% of the cases, being amongst the less frequent causes for BBS [17]. In 2003, Ansley et al. discovered that the *TTC8* gene is associated with BBS and identified the syndrome as a basal body dysfunction of the ciliated cells [2]. The *TTC8* gene, also referred to as *BBS8*, is one of the eight *BBS* genes (*BBS1*, *BBS2*, *BBS4*, *BBS5*, *BBS7*, *TTC8*, *BBS9*, and *BBIP1*), encoding proteins that assemble into a stable octameric protein complex termed the BBSome [18,19]. The BBSome forms a membrane coat that sorts membrane receptors to the primary cilium and its dysfunction leads to the failure of cell-specific signal transduction affecting multiple cell-types and organs [20].

In addition to pleiotropic disorder, there are also reports of non-syndromic retinal degeneration caused by mutations in the BBS genes *BBS1*, *BBS2*, *ARL6*/*BBS3*, and *TTC8* [21–24], disrupting the normal function of photoreceptor cilia. At the time of the discovery of the canine *TTC8* deletion in 2014, there were indications that dogs homozygous for the deletion might exhibit clinical signs other than retinal degeneration [1]. However, only ophthalmic examinations of the affected dogs were performed and other BBS-associated characteristics could not be clinically investigated. In addition, it was not possible to analyze tissues from affected individuals. Here, we describe a detailed examination of two golden retrievers, homozygous for the *TTC8* variant (*TTC8*delA) that were followed from the time of the PRA-diagnosis until they were euthanized and a full necropsy was performed. This allowed for a rigorous clinical characterization of this canine form of a *TTC8*-mediated disease, and to investigate the effect of the canine *TTC8* mutation on the transcriptional and protein level in the canine retina.

#### **2. Materials and Methods**

#### *2.1. Animals and Samples*

A golden retriever sib-pair, a male (GR01) and a female (GR02), was followed from the time of their PRA-diagnosis until euthanasia, after which a necropsy was conducted and tissue samples were collected. The male dog was 6 years and 1 month old, and the female was 3 months older at the time of euthanization. Both dogs tested homozygous for the *TTC8*delA allele. Tissue samples were also collected from five unaffected dogs (RW01, BE02, LR02, GS01, and GSP01) euthanized for reasons unrelated to this study (Table S1). In addition, two unaffected dogs (BE01 and LR01; Table S1) were included in the clinical ophthalmic examination (see below). Interviews were conducted with owners of the sib-pair, as well as the owners of eight additional affected golden retrievers, all tested homozygous for the *TTC8*delA allele (Table S2). All samples were obtained with informed consent from the dog owners. Ethical approval was granted by the regional animal ethics committee (Uppsala ethics committee on animal experiments/Uppsala djurförsöksetiska nämnd; Dnr C12/15, Dnr 5.8.18-15533/2018, and Dnr 5.8.18-04682/2020).

#### *2.2. Ophthalmic Examination, Confocal Scanning Ophthalmoscopy (cSLO), Optical Coherence Tomography (OCT) and Electroretinography (ERG)*

Ophthalmic examination of GR01 and GR02, as well as an unaffected beagle (BE01) and an unaffected Labrador retriever (LR01) included reflex testing, testing of vision with falling cotton balls under dim and daylight conditions, and tonometry (Tonovet, Icare Finland Oy, Vantaa, Finland), as well as indirect ophthalmoscopy (Heine 500, Heine Optotechnik GmbH, Herrsching, Germany) and slit-lamp biomicroscopy (Kowa SL-17, Kowa Company Ltd., Tokyo, Japan) after dilation of pupils with tropicamide (Mydriacyl 0.5%, Alcon Nordic AS, København, Denmark).

cSLO- and OCT-imaging (Spectralis HRT + OCT, Heidelberg Engineering, Heidelberg, Germany) were performed after dilation of pupils with tropicamide and under light sedation with 5 μg/kg medetomidine (Sedator, Dechra Veterinary Products AB, Upplands Väsby, Sweden) and 50 μg/kg butorphanol (Dolorex, Intervet AB, Stockholm, Sweden). The cornea was kept moist using artificial tears throughout the procedure. Total retinal thickness in GR02 was measured as previously described [25], and compared to data from an unaffected 5-year-old female Labrador retriever (LR01) and a 7-year-old female beagle (BE01).

We recorded a bilateral, full-field ERG in GR02 under general anesthesia and compared to data from the unaffected LR01. Sedation with intramuscular acepromazine 0.03 mg/kg (Plegicil vet., Pharmaxim Sweden AB, Helsingborg, Sweden) was followed by induction with propofol 10 mg/kg, intravenously (Propovet, Orion Pharma Animal Health AB, Danderyd, Sweden). After intubation, inhalation anesthesia was maintained with isoflurane (Isoflo vet., Orion Pharma Animal Health AB). Corneal electrodes (ERG-JET, Cephalon A/S, Aalborg, Denmark) were used with isotonic eye drops (Comfort Shield, i.com medical GmbH, Munich, Germany) as a coupling agent. Gold-plated, cutaneous electrodes served as ground and reference electrodes (Grass, Natus Neurology Inc. Middleton, WI, USA) at the vertex and approximately 3 cm caudal to the lateral canthi, respectively. Light stimulation, calibration of lights, and processing of signals were performed as previously described [26] and the ECVO protocol for canine clinical ERGs was followed [27].

#### *2.3. Clinical Andrological Examination and Semen Analysis*

Semen from the affected male (GR01, age 6 years and 1 month), was collected using digital manipulation with a bitch in estrus present. The sperm concentration was measured using a Bürker chamber. Sperm morphology was evaluated using standard procedures in wet preparations of semen fixed in buffered formalin and in air-dried smears stained with carbolfuchsin–eosin. Due to a very low sperm concentration and volume, only 50 spermatozoa were examined in the wet preparation, under a phase-contrast microscope at 1000× magnification. All abnormalities on any given spermatozoon were counted and the overall frequencies were classified according to a system developed by Bane (1961) [28]. For a more detailed examination of the sperm heads, 100 spermatozoa were evaluated in a smear under a light microscope at 1000× magnification. Presence of spermatogenic cells was recorded in smears stained in accordance with standard well-established methodology. The head morphologies were classified as pear-shaped, narrow at the base, abnormal contour, undeveloped, loose and abnormal, narrow, and variable size.

#### *2.4. Echocardiography and Electrocardiographic (ECG) Examinations*

For the echocardiographic examinations of the sib-pair (GR01 and GR02, at an age of 6 years and 1 month), the dogs were placed in right and then left lateral recumbency on an ultrasound examination table. The echocardiographic evaluation was conducted by use of an ultrasonographic unit (EPIQ 7G, Philips Ultrasound, Bothell, WA, USA) equipped with a 5-1 matrix transducer and ECG monitoring. The heart was examined and subjectively assessed in standard right- and left-sided views [29]. Blood flow over heart valves was interrogated using color mode Doppler echocardiography and measured using spectral Doppler echocardiography. Left ventricular dimensions were measured using M-mode echocardiography in the right parasternal short axis view at the level of the papillary muscles. The left atrial diameter was measured in the right parasternal short axis view at the level of the aortic valve. Left ventricular dimensions were compared to published weight-based normal reference ranges [30] and left atrial diameter was indexed to the aortic diameter as previously described [31].

#### *2.5. Histopathological Examinations*

Tissue samples from the affected siblings GR01 and GR02 (Table S3), as well as a retinal sample from an unaffected 3-year-old Rottweiler (RW01), were fixed in 10% neutral buffered formalin for >24 h, paraffin wax embedded, sectioned at 4 μm and stained with hematoxylin and eosin (HE). In addition, kidney sections were also stained with periodic acid-Schiff (PAS). All sections were evaluated using light microscopy.

#### *2.6. Total RNA Extraction*

The retinal tissue sample from the affected female GR02 was collected directly in TRIzol Reagent (Invitrogen™, Carlsbad, CA, USA) for immediate RNA extraction. Retinal samples from three unaffected dogs—a beagle (BE02), a Labrador retriever (LR02), and a German shepherd (GS01)—were preserved in RNAlater (SigmaAldrich, Saint Louis, MO, USA) directly after euthanasia. The samples were then washed with 1 × PBS and between 50 to 100 mg of tissue was used for extraction. The samples from GR02, BE02, and GS01 were then homogenized in TRIzol Reagent with Precellys homogenizer (Bertin Instruments, Montigny-le-Bretonneux, France) and total RNA was extracted following the manufacturer's instructions (Pub. No. MAN0001271, Rev. B.0.). Total RNA from LR02 was extracted with RNeasy mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. RNA concentration was measured using Qubit RNA BR Assay kit (Invitrogen™, Waltham, MA, USA). RNA integrity and quality were inspected with Agilent TapeStation 4150 with an Agilent RNA ScreenTape (Agilent Technologies, Santa Clara, CA, USA). PolyA-selection was carried out utilizing Dynabeads mRNA Purification Kit (Invitrogen™, Waltham, MA, USA), applying the manufacturer-provided protocol.

#### *2.7. Long-Read cDNA Sequencing with Oxford Nanopore Sequencing Technology (ONT)*

Oxford Nanopore long-read cDNA sequencing libraries for the polyA-selected RNA of the canine retina from the affected GR02 as well as unaffected BE02 and GS01 were prepared using the Direct cDNA Sequencing kit (SQK-DCS109, Oxford Nanopore Technologies Ltd., Oxford, UK) following the manufacturer's instructions (Protocol version DCS\_9090\_v109\_revC\_04Feb2019; updated 02 May 2019). The input material for each library was 100 ng of polyA-selected RNA. The prepared libraries were sequenced on a MinION sequencer using R9.4.1 flow-cells and MinKNOW (v.19.06.7) software. The resulting raw fast5 files were basecalled using Guppy version 3.1.5+781ed57 with the –trim\_strategy DNA flag. The quality control for each run was performed using an R markdown script developed by ONT (Nanopore\_SumStatQC\_Tutorial.Rmd). To characterize the retinal *TTC8* transcripts, quality passed reads (quality score > 7) were mapped to the canine reference genome sequence CanFam3.1 using MiniMap2 (v.2.16) with -ax splice parameter (Li, 2018) [32]. A reference-guided assembly of the transcriptome was produced using StringTie2 (v.2.1.1) with default settings for long-reads (-L) [33], with Ensembl anotation build 100 (CanFam3.1). The assembled transcripts were compared to the raw sequnce data using The Integrative Genomics Viewer (IGV) [34,35]. Transcripts with less then 20× coverage were discared. To quantify the expression levels, quality passed reads (quality score > 7) were mapped to the canine reference transcriptome consisting of protein

coding gene sequences and non-coding RNA sequences from Ensembl build 100 (CanFam3.1) using Minimap2 (v2.16) with the settings -ax map-ont -N 100 [32]. The transcript counts were produced using Salmon (v.1.1.0) in aligned based mode and –libType A for automatic detection of library type [36]. The results for each transcript expressed as transcripts per million (TPM) were summarized to produce gene level TPM values. We then analyzed the expression levels of genes considered to be expressed in a cell-specific manner (retinal marker genes) [37–57]. The complete list of marker genes and their respective references can be found in Table S4.

#### *2.8. Quantitative RT-qPCR*

cDNA was synthesized using RNA prepared from retinal samples of the affected female (GR02) and three unaffected dogs (BE02, LR02, and GS01) using RT<sup>2</sup> First Strand kit (Qiagen, Hilden, Germany) with random hexamers provided in the kit. cDNA concentration was measured using Qubit ssDNA Assay kit (Invitrogen™, Waltham, MA, USA). To amplify the target and the reference genes, custom primers were designed using the software Primer3 [58,59] (Table S5). The cDNA fragments encoding regions of interest were amplified using RT<sup>2</sup> SYBR Green ROX qPCR Mastermix (Qiagen, Hilden, Germany) on a StepOnePlus Real-Time PCR system (Applied Biosystems™, Waltham, MA, USA), according to the manufacturer's instructions. Target gene expression was normalized to expression of *GAPDH* as well as *ACTB* reference genes, and shown relative to a control sample (BE02) (CT method). The results were confirmed in two independent experiments.

#### *2.9. Fluorescence Histochemistry*

Retina from the affected female (GR02) and an unaffected 10-year-old male German spaniel (GSP01) were analyzed by means of fluorescence histochemistry. The lens and anterior segment were removed and the vitreous punctured before the eye cups were fixed in 4% PFA in 1 × PBS for 15 min on ice, washed in 1 × PBS for 10 min on ice, and cryoprotected in 30% sucrose in 1 × PBS at 4 ◦C until saturated. The eye cups were embedded in Neg-50™ frozen section medium (Thermo Scientific, Waltham, MA, USA) and 10 μm sections were collected on Superfrost Plus slides (J1800AMNZ, Menzel-Gläser, Thermo Fisher Scientific, Waltham, MA, USA). Sections were re-hydrated in 1 × PBS for 10 min, incubated in blocking solution (1% donkey serum, 0.02% thimerosal, and 0.1% Triton X-100 in 1 × PBS) for 30 min at room temperature, and incubated in blocking solution containing Alexa Fluor™ 488-conjugated PNA (1:400, L21409, Invitrogen™, Waltham, MA, USA) and a rhodopsin primary antibody (1:5000, NBP2-25160, Novus Biologicals, Abingdon, UK) overnight at 4 ◦C. The sections were then washed in 1 × PBS, 3 × 5 min, incubated in Alexa 568 secondary antibody (1:2000, A10037, Invitrogen™, Waltham, MA, USA) for at least 2 h at room temperature, washed in 1 × PBS, 3 × 5 min, and the slides mounted using ProLong™ Gold Antifade Mountant with DAPI (P36931, Molecular Probes, Waltham, MA, USA). Images were captured using a Zeiss Imager.Z2 microscope equipped with an Axiocam 512 mono camera (Carl Zeiss Microscopy GmbH, Jena, Germany).

#### **3. Results**

#### *3.1. Desciption of Primary and Secondary Bardet-Biedl Syndrome Characteristics*

To investigate if the dogs exhibit other BBS-related problems, we interviewed the owners of 10 affected dogs. Two of the dogs were also part of the 2014 study [1], and we based the interviews on the questionnaire developed by Downs et al. (2014) [1]. Among the 10 affected dogs, 4 were reported to exhibit a minimum of four characteristics matching human primary BBS-signs (Table S2). The remaining dogs displayed three (two dogs) and two (three dogs) primary signs, respectively, and for one dog, PRA was the only reported primary characteristic. All of the dogs had been diagnosed with PRA, and the average age of diagnosis was 4 years and 8 months. None of the dogs were reported to have polydactyly or any other malformations of the paws, but obesity or polyphagia, renal problems, cognitive dysfunction, irregular estrous cycles, and decreased libido in males were among the reported

primary signs. As for secondary BBS signs, many dogs were described to have a short stature relative to the standard for the breed. Partial or complete loss of the sense of smell (anosmia), worsening with age, was noticed by most of the owners (Table S2). Half of the individuals were also diagnosed with cataracts. At the time of this study, seven of these dogs were no longer alive. The oldest dog was reported to have died in his sleep at 10 years of age, while the other dogs had been euthanized at an average age of 7.5 years. Among the dogs in the questionnaires was the BBS-affected sib-pair (GR01 and GR02), which were further investigated in this study.

#### *3.2. Ophthalmic Examinations of the A*ff*ected Sib-Pair*

The male dog GR01 was initially diagnosed with PRA at the age of 4 years and 3 months, after which he and his female littermate GR02 were genotyped and found to be homozygous for the *TTC8*delA allele and subsequently recruited for this study. The two dogs were then examined ophthalmologically on three occasions: at the age of 4.5 and 5.5 years, and finally at the age of 6 years and 1 month. At the first examination, both dogs had mildly dilated pupils in room light. Pupillary light reflexes and menace responses were considered normal, but the dazzle reflex was mildly impaired in both eyes (oculi uterque, OU). The results of the cotton ball test in room light were unremarkable for both dogs, whereas the female detected approximately 50% and the male only between 10–20% of the cotton balls under dim light conditions. On indirect ophthalmoscopy, generalized abnormal tapetal reflection going from hypo- via normo- to hyperreflection (in the male) was seen when the indirect ophthalmoscopy lens was tilted back and forth. Furthermore, retinal vascular attenuation and pigment clumping in the non-tapetal fundus were observed. Compared to the female dog, the male showed signs of more advanced retinal degeneration. Findings were symmetrical between eyes. In addition to the abnormal retinal appearance, equatorial and anterior cortical cataracts were seen in the eyes of both individuals. The male had multiple, pigmented iridociliary cysts mainly emerging through the nasal portion of the pupils (OU), but free-floating and collapsed cysts were also seen in the anterior chambers. Intraocular pressures (IOP) were normal (OU) in both dogs.

To further investigate the retinal phenotype, we performed cSLO and OCT on the affected female (GR02). Reduced fundus autofluorescence (FAF) in the entire tapetal fundus was seen on cSLO (Figure 1a), whereas the non-tapetal fundus appeared slightly hyperreflective indicating accumulation of lipofuscin in the degenerating retina. Retinal vascular attenuation in the affected fundus was evident, as well as islets with bright tapetal hyperreflection (Figure 1b), compared to the unaffected dog (Figure 1c,d).

**Figure 1.** Fundus autofluorescence (FAF) and infrared confocal scanning laser ophthalmoscopy (IR-cSLO). cSLO FAF- and IR-cSLOs of the left eye of GR02 (panels (**a**,**b**), respectively) and age-matched, unaffected Labrador retriever LR01 (panels (**c**,**d**), respectively). Both FAF-images are inverted and hence, the tapetal reflection in GR02 (**a**) is remarkably faint. Detection of FAF in the tapetal fundus is difficult because of the normal, bright reflection from the tapetum lucidum. Black arrows indicate areas with bright reflection in the tapetal fundus on FAF (**a**), but hyporeflection on IR (**b**), which is suggestive of lipofuscin accumulation. The slight hyperreflection (the grey shade below the optic nerve head in the inverted image) (**a**) in the non-tapetal area autofluorescence indicates storage of lipofuscin. The IR-image (**b**) shows the large, myelinated optic nerve head often seen in the golden retriever breed, but very attenuated retinal vessels (red arrows) and islets of hyperreflectivity. Cast shadows in the cSLOs from GR02 were caused by the cataracts present in this dog (one shadowed area indicated by the blue oval).

OCT showed considerable retinal thinning (Figure 2a) compared to two unaffected dogs (LR01 and BE01). Mainly the outer retinal layers were reduced in thickness and outer retinal landmarks, such as the external limiting membrane (ELM), ellipsoid zone (EZ), and interdigitation zone (IZ), could not easily be identified in the affected dog (Figure 2b). The inner nuclear layer (INL) was thin and often fragmented, making the segmentation of the outer plexiform layer (OPL) and inner plexiform layer (IPL) difficult. In the ventral (inferior) non-tapetal fundus, retinal thickness was also reduced, but more irregularly, giving a patchier appearance. Again, segmentation of the outer retinal layer on the OCT image was difficult and a distinct outer nuclear layer ONL occasionally missing. Both fragmentation and thickening of the retinal pigment epithelium (RPE)was observed (Figure 2c). The cataracts in the male (GR01) hindered the cSLO and OCT examinations.

**Figure 2.** Retinal thickness and optical coherence tomography (OCT). (**a**) Total retinal thickness temporally from the optic nerve head (ONH) towards the periphery of the retina is reduced in the affected dog (GR02) compared to the thickness observed in an unaffected Labrador retriever (LR01) and a Beagle (BE01). (**b**) OCTs from the nasal part of the fundus (red rectangles) from the affected dog (GR02, top) and the unaffected dog (LR01, bottom) showing both the marked thinning of the outer nuclear layer (ONL) (asterisk) and less distinct segmentation of retinal layers in the affected dog. (**c**) Appearance of the retina in the non-tapetal, ventral fundus of GR02. The vertical red bar in the cSLO-image shows the area outlined by the red rectangle on the OCT photograph, which is then magnified below (the area inside the large red rectangle). Cast shadows on the cSLO-image were caused by cataracts. The thinned retina is more variable in thickness than in the tapetal fundus dorsal to the ONH. The retinal pigment epithelium (RPE) is occasionally fragmented (red arrows) or appears thickened (blue arrows). The segmentation of the outer retina is difficult and dark areas probably representing clusters of nuclei in the ONL are seen occasionally (white arrows).

Flash-electroretinography (FERG) was performed in the female, but not in the male, because of the ophthalmoscopic signs of more advanced retinal degeneration. In the female, rod responses were non-detectable throughout the 20 min dark adaptation. Neither was the a-wave of the mixed, dark-adapted rod-cone response discernible and the b-wave had profoundly subnormal amplitude and was dominated by the early, mainly cone-driven part. The light-adapted cone-driven transient b-wave was relatively better preserved, but biphasic with a late, second peak, whereas the a-wave was

essentially missing. The cone flicker response was more abnormal, as the cone-driven responses were unable to follow the 30 Hz stimuli (Figure 3).

**Figure 3.** Flash-electroretinographs (FERGs). Both dark-adapted (DA) and light-adapted (LA) responses of the FERG had subnormal amplitudes in the affected golden retriever (GR02; red tracings) compared to the unaffected Labrador retriever (LR01) (blue tracings). Rod responses were essentially non-detectable and even the dark-adapted responses to the bright rod-cone stimulus seem to be mainly cone-driven. Light-adapted cone transients lacked a-wave and had a biphasic waveform. When a 30 Hz stimulus was employed, the cone responses of GR02 came out of sync and were not time-locked to the individual stimuli. Calibrations show the time base and the different amplitude scales for GR02 and LR01.

The siblings were re-examined 11 months later at the age of 5.5 years. Pupils were now moderately to widely dilated in both dogs, pupillary light reflexes sluggish and incomplete, and dazzle reflexes impaired. The male dog showed poor menace responses and did not respond to falling cotton balls regardless of lighting, but did occasionally follow hand movements at 30 cm in bright light. Vision was less impaired in the female; menace responses were retained and approximately 50% of the falling cotton balls were perceived in room light, whereas no cotton balls were detected in dim light. Ophthalmoscopy showed generalized tapetal hyperreflection in both dogs with a brighter hyperreflection in the male than in the female. Retinal vascular attenuation was now moderate to advanced and the optic nerve heads appeared pale. In the non-tapetal area, pigment clumping was observed, as on the previous examination. Cataractous changes had progressed in both dogs and were seen in both the anterior and posterior cortices of both eyes, as well as in the equatorial region. Iridociliary cysts were still present in the male and his eyes were normotensive.

Both dogs were examined a third time at the age of 6 years and 1 month. At this time, the male was severely visually impaired both under daylight and dim light conditions, and had developed bilateral mature-hypermature total cataracts. In addition, the dog showed both intact and ruptured iridociliary cysts bilaterally (Figure 4a) and had normal IOP. However, the ocular fundi could not be examined in vivo because of the cataracts. Postmortem gross examination of the fundi revealed bright hyperreflection in the tapetal area (Figure 4b). The network of retinal blood vessels was severely attenuated (even for a postmortem specimen) and difficult to follow towards the periphery (Figure 4b). The pupils of the female dog (GR02) were dilated in room light conditions, the pupillary light reflexes (PLRs) were sluggish and incomplete, dazzle reflexes were poor, and menace responses were present in room light, but not in dim light conditions. The female was unable to detect falling cotton balls regardless of lighting. The cataractous changes had progressed (OU), particularly in the posterior than in the anterior cortex, and the tapetal hyperreflectivity and retinal vascular attenuation could now only be seen like through a haze.

**Figure 4.** Iridociliary cysts and fundus imaging of the affected male (GR01). (**a**) Iridociliary cysts both at the posterior side of the iris and free-floating in the anterior chamber were seen in both eyes of the affected male (GR01). Pigment and remnants of ruptured cysts were seen on the anterior capsules of the cataractous lenses. (**b**) The eyecup of the affected male (GR01), *postmortem.* Bright tapetal hyperreflection is seen close to the optic nerve head (ONH), whereas the rest of the tapetum is faintly colored. The retinal vessels, including the larger venules, are almost impossible to follow.

#### *3.3. Andrological Examination*

The testes of the affected male dog GR01 were both palpable in the scrotum and found to be smaller than normal for the age and breed. An ejaculate of a very small volume (<0.5 mL), and light grey in color, was obtained. The sperm concentration was 15 <sup>×</sup> 106/mL and motility was <sup>&</sup>lt;5%. The total sperm count could not be calculated exactly, but was estimated to be <7 <sup>×</sup> 106/mL. There was a high proportion of spermatozoa with abnormalities, and less than 1% were considered morphologically normal. Midpiece defects were predominant (82% of spermatozoa), consisting of "Dag defect" [60], with strong folding, coiling, and fracture of the distal part of the midpiece, and "tail stump defect", where, in place of normal tails, short "stumps" were found. Other defects were coiled tails (12%), and double bent tails (6%), as well as proximal droplets (12%), acrosome defects (4%), and nuclear pouches (20%). Head defects were detected in 36% of the spermatozoa, the most common being pear shaped (11%), variable size (10%), and loose abnormal heads (7%). In addition, a large number of spermatogenic cells, inflammatory cells, and epithelial cells, as well as necrotic cells of varying sizes, were present (Figure 5a).

**Figure 5.** Light microscopy of sperm and seminiferous tubules from the affected male (GR01). (**a**) Papanicolaou staining of semen showing several abnormal spermatozoa. Different tail defects are evident, including coiled tails (C), tail stumps (\*) and loose heads (L). S: spermatogenic cell, N: neutrophil leukocyte. (**b**) Tubuli seminiferi in testis stained with hematoxylin and eosin, showing mainly early stages of spermatogenesis and few elongated spermatids (ES).

#### *3.4. Cardiac Examination*

Both dogs (GR01 and GR02) presented normal respiratory sinus arrhythmia at a normal heart rate (between 80–90 BPM). Echocardiography showed normal cardiac morphology, dimensions, and motion in both dogs. Minimal regurgitation over the mitral valve was detected in both dogs, but valve morphology appeared normal. None of the dogs manifested *situs inversus*.

#### *3.5. Necropsy Findings*

To further investigate clinical features associated with BBS, a thorough necropsy was conducted on both affected dogs. The body condition score (BCS) of the male (GR01) was estimated to be 9/9, consistent with obesity [61]. In gross examination of the male, an abnormal head shape with a broadened muzzle and flat forehead was perceived, as was a low withers height. Furthermore, there was gingival hyperplasia and moderately increased interdental spaces (diastema) in both lower and upper jaw. The testicles were also seemingly small for the dog's breed and age. Gross examination of the kidneys revealed suspected chronic infarcts bilaterally, with the most marked changes in the right kidney (Table S3). Mild myxomatous valvular degeneration was seen in the heart of the affected male. The body condition score (BCS) of the female (GR02) was assessed to be 7/9, consistent with heavy [61]. The interdental spaces were also moderately increased. The kidneys showed signs of mild chronic glomerulonephritis bilaterally (Table S3).

#### *3.6. Histopathology Findings*

Next, we examined the histopathology of tissue samples collected during the necropsy. Abnormal changes were observed in the retina, kidney, and testis. The majority of the retina from the male (GR01) exhibited severe retinal thinning with complete loss of normal architecture with most pronounced degenerative changes in the non-tapetal fundus. Compared to a normal retina (Figure S1a), the ONL of the affected male GR01 was severely affected in the tapetal fundus and completely missing in the non-tapetal fundus (Figure S1b,c).

Histopathological examination of the testicles revealed only few (or no) late spermatids in the tubuli seminiferi (Figure 5b). A reduced number of spermatocytes was also observed. Microscopically, sections from the lesions in the right kidney of the male (GR01) corresponded to segmental areas of fibrosis, infiltrated by a moderate number of lymphocytes and plasma cells, extending from cortex to medulla (Figure 6a,b). In the fibrotic tissue, there were some degenerated glomeruli and lack of tubular

structures. Furthermore, a disarray of occasional small glomeruli with peripheral nuclei and inapparent capillaries (fetal glomeruli), as well as atypical tubular structures outlined by pseudostratified cuboidal epithelium with large, basophilic, plump cells, were seen. No cellular atypia or mitotic activity were evident. Occasional dilated tubular structures were seen (Figure 6b). No chronic infarcts were confirmed histologically in the left kidney. On microscopic examination of the kidneys from the female (GR02), changes consistent with mild chronic glomerulonephritis was evident.

**Figure 6.** Renal histopathology of the affected male dog (GR01). Light microscopy of the right kidney stained with hematoxylin and eosin, showing a segmental area of fibrosis, (**a**) infiltrated by lymphocytes and plasma cells, extending from cortex, with the presence of a disarray of atypical tubular structures, and (**b**) with occasional small glomeruli with peripheral nuclei and inapparent capillaries (fetal glomeruli, marked with an arrow). Furthermore, multiple degenerated glomeruli are seen (arrowheads).

#### *3.7. Characterization of Canine TTC8 Transcripts in the Retina*

To comprehensively characterize the different *TTC8* transcripts expressed in the canine retina and to investigate the effect of the 1 bp deletion (*TTC8*delA) in an affected dog, we performed full-length cDNA sequencing of the neural retina from the female (GR02), as well as of two 12-year-old unaffected female dogs, a beagle (BE02) and a German shepherd (GS01), using Oxford Nanopore Technologies (ONT) for long-read sequencing. This produced 3.6 M reads and 5.86 Gb of quality passed DNA sequence for GR02, 7.8 M reads and 10.68 Gb for BE02, as well as 7.6 M reads and 7.06 Gb for GS01 (Table 1).


The cDNA sequence data was assembled into transcripts using reference-guided assembly. The predicted transcripts of the *TTC8* gene were then manually curated with the raw sequence data. We could not assemble any full-length *TTC8* transcripts from the sample of the affected dog (GR02) due to lack of read coverage over the *TTC8* locus. In contrast, sequence data from the unaffected dogs (BE02 and GS01) suggested transcription of three different retinal isoforms at the *TTC8* gene, here denoted tr1, tr2, and tr3 (Figure 7a). Fourteen exons were detected in tr1, 15 in tr2, and 16 in tr3, comprising 6.5%, 76.1%, and 17.4% of the *TTC8* transcripts, respectively. In tr1, exons 1b and the retinal specific exon 2a were missing. In tr2, exon 1b was absent, whereas tr3 featured both exons 1b and 2a. The comparison of the identified transcripts to existing annotations showed that the exon-intron boundaries and

sequence between the first and the last exon of tr1 corresponds to the curated NCBI RefSeq transcript NM\_001284469.1, tr2 corresponds to the predicted NCBI RefSeq transcript XM\_003639207.4, as well as Ensembl transcript ENSCAFT00000027700.5, and tr3 corresponds to the predicted NCBI RefSeq transcript XM\_014115661.2. The exon-intron boundaries of the identified transcripts show similarity to human *TTC8* transcripts (Figure 7b). In all three transcripts, a putative translation initiation site (TIS) was identified in exon 1 (Chr8:60,061,732) containing an AUG start codon with an adjacent optimal (GCCRCCAUGG) Kozak consensus sequence (Figure 7a) [62,63]. Using the TIS, tr1 and tr2 have an open reading frame extending to the last exon, encoding for 505 and 515 aa protein products, respectively. However, in tr3, the TIS would lead to a termination codon two exons downstream of exon 1, producing a 49 aa protein product and would lack the tetratricopeptide repeat (TPR) motifs. An alternative TIS is found at position Chr8:60,077,223 of exon 2 and, although the adjacent Kozak consensus is not optimal, it is classified as strong (NNNRNNAUGG) [63]. This alternative TIS is in-frame with the original coding sequence, and its open reading frame would produce a 455 aa protein product for all three transcripts, including the TPR motifs (Figure 7a). We did not observe any transcripts skipping exon 7, where the deletion Chr8:60,090,185delA is located, nor any exons downstream the deletion.

**Figure 7.** Retinal isoforms of canine the *TTC8* gene and human variants associated with Bardet-Biedl syndrome (BBS). (**a**) On top, a schematic representation of the TTC8 protein including the approximate

location of the TPR motifs, and below, the three identified canine retinal transcripts. The black dashed lines indicate the two alternative TIS including the Kozak consensus sequence, and the translation termination site in exon 14. The red dashed lines indicate the 1 bp deletion of adenine in exon 7 and the premature stop codon in exon 8. (**b**) An overview of identified human *TTC8* coding variants (red dashed line) and splice-site variants (blue dashed line), see Table S6. The positions of the human mutations are based on the reported positions in the original publications, and the exon numbering may therefore differ from the positions in the presented transcript (ENST00000345383.9). tr1: transcript 1, tr2: transcript 2, tr3: transcript 3, TPR: tetratricopeptide repeat, TIS: translation initiation site.

Next, we compared *TTC8* gene expression levels of the unaffected BE02 and GS01 with the *TTC8* levels of the affected GR02 with 1 bp deletion in the exon 7 (Figure 8). We first estimated the transcript abundance for each annotated gene with Salmon (v.1.1.0) [36] using canine transcriptome from Ensembl build 100 and then summarized the level of transcript expression of each gene in the annotation. In BE02, the *TTC8* gene was expressed at a level of 99.8 transcripts per million (TPM) (856th most highly expressed gene), and in GS01 at a level of 60.6 TPM (859th most highly expressed gene). As suggested by the low read coverage of GR02 of *TTC8* transcripts, the expression levels were low in the affected dog, with 16.8 TPM (3625th most highly expressed gene). We then summarized the expression levels of known marker genes for different cell types in the retina (Table S4, Figure 8). The expression of genes transcribed in the rod photoreceptor cells (*PDE6A*, *PDE6B*, *CNGB1*, *CNGA1*, and *GNAT1*) [39,43,64–66] accounted for 25–27% of the total marker gene expression (TPM) in both the unaffected dogs (Figure 8a,b). For the affected GR02, the expression was estimated to 0.2%, indicating that the sample only included a small fraction of rod cells compared to the unaffected dogs (Figure 8c). The rhodopsin gene (*RHO*), encoding for a specialized G protein-coupled receptor known to be expressed exclusively in rod photoreceptors, was not included in the most recent Ensembl annotation release (build 100) used in the analysis, and therefore not included in the rod marker gene list. Similarly, the POU Class 4 Homeobox 1 (*POU4F1*), marker for retinal ganglion cells [67], was not among the genes annotated in Ensembl build 100. Both *RHO* and *POU4F1* were instead included in the quantification using reverse transcription quantitative real-time PCR (see below). The expression levels of cone photoreceptor cell markers (*ARR3*, *GUCA1C*, *PDE6C*, *PDE6H*, and *OPN1LW*) [37–43] were also lower in GR02, whereas the expression levels of macroglial cell (MG) marker genes (*RLBP1*, *SLC1A3*, *GLUL*, *CLU*, and *GFAP*) [43–46] were elevated (Figure 8a–c). The marker gene expression levels of the horizontal cells (HC) [47,48], retinal ganglion cells (RGC) [49–54], RPE cells [50,55], amacrine cells (AC) [47,56], or bipolar cells (BC) [57] did not differ drastically between the three dogs, although the proportional expression of BC and RPE markers was higher in the affected female. The complete list of marker genes and their expression values in each sample can be found in Table S4.

**Figure 8.** Expression of retinal marker genes. (**a**) The average expression of marker genes for each retinal cell type in the unaffected BE02 (**a**) and GS01 (**b**) dogs as well as the affected GR02 (**c**), based on cDNA sequencing. (**d**) Relative mRNA expression levels by quantitative RT-qPCR in two different regions

(exons 7–8 and exons 13–14) of the *TTC8* gene, as well as the retinal marker genes *GFAP* (macroglial cells), *OPN1LW* and *OPN1SW* (cone photoreceptors), *RHO* (rod photoreceptors), *RPE65* (RPE cells), and *POU4F1* (retinal ganglion cells) expression in three unaffected dogs (BE02, LR02, GS01) and the affected GR02, normalized to *GAPDH* and *ACTB* gene expression. Rod: rod photoreceptors, cone: cone photoreceptors, RGC: retinal ganglion cells, AC: amacrine cells, HC: horizontal cells, BC: bipolar cells, MG: macroglial cells, RPE: retinal pigment epithelium, BE02: unaffected beagle, GS01: unaffected German shepherd, GR02: affected golden retriever, LR02: unaffected Labrador retriever.

To verify the relative differences in the expression levels of *TTC8* and retinal marker genes in the affected GR02 and unaffected BE02 and GS01, as well as in an unaffected Labrador retriever (LR02, not used in the cDNA sequencing), we performed reverse transcription quantitative real-time PCR (RT-qPCR). We amplified two separate regions of the *TTC8* gene over the exons 7 to 8 and 13 to 14 (Figure 8d). In addition, we amplified canine long-wave (*OPN1LW)*, and short-wave (*OPN1SW*) opsins expressed in cone photoreceptors, as well as, retinoid isomerohydrolase RPE65 (*RPE65*) and glial fibrillary acidic protein (*GFAP*) genes (markers for RPE cells and macroglia, respectively) to compare their cDNA sequencing levels to RT-qPCR. We also used *RHO* to evaluate relative rod photoreceptor levels and POU Class 4 Homeobox 1 (*POU4F1*) to evaluate retinal ganglion cells in these retinal samples. With the exception of *RPE65* expression, the RT-qPCR results reflected the cDNA sequencing quantification showing that the expression of both *TTC8* amplicons, as well as *RHO* expression, were lower in the GR02 compared to the three unaffected dogs and *GFAP* expression was higher than in the unaffected dogs (Figure 8d).

#### *3.8. Fluorescence Histochemical Analysis of the TTC8delA Retina*

RT-qPCR analysis revealed low expression of markers of rod and, to a lesser extent, cone photoreceptors in the affected female (GR02; Figure 8), and the histopathology analysis for the affected male (GR01; Figure S1) and OCT analysis for the affected female (GR02; Figure 2) showed that both the ONL and photoreceptor layer were thinner compared to control retinas. To investigate the presence of photoreceptors in the *TTC8*delA retina, we performed fluorescence histochemistry using an antibody directed against rhodopsin and fluorophore-conjugated peanut agglutinin (PNA), which binds selectively to cone photoreceptors in the retina [68]. Retinal sections from an unaffected, 10-year-old German spaniel (GSP01) and the affected, 6-year-old female golden retriever (GR02) were analyzed (Figure 9). We found that the ONL and the photoreceptor layer were thinner in the GR02 retina compared to control. In the unaffected retina, both rhodopsin and PNA staining were found in the photoreceptor layer (Figure 9a), whilst in the affected GR02 retina, rhodopsin staining was absent and PNA labeled what appeared to be truncated cone outer segments in the photoreceptor layer (Figure 9b). The fluorescence data corroborate the histopathology and OCT findings.

**Figure 9.** Fluorescence histochemistry of dog retinas for rhodopsin and cone photoreceptors. Fluorescence micrograph showing rhodopsin expression (red), Alexa™ 488-conjugated PNA (green), and DAPI (blue) in (**a**) a male, unaffected, 10-year-old German spaniel (GSP01) and (**b**) a female, affected, 6-year-old golden retriever (GR02). Note the thin ONL, lack of rhodopsin staining, and the truncated cone outer segments (exemplified with arrowheads) in the retina of the affected compared that of the unaffected dog. PNA: peanut agglutinin, ONL: outer nuclear layer, INL: inner nuclear layer, GCL: ganglion cell layer. Scale bar: 50 μm and applies to both images.

#### **4. Discussion**

When the deletion at the *TTC8* locus (*TTC8*delA) in golden retrievers was discovered in 2014, there were indications that the mutation did not only lead to progressive retinal degeneration in the affected dogs, but could also cause additional clinical signs, suggesting a syndromic type of disease [1]. However, at the time, no cases were available for a thorough clinical examination apart from ophthalmoscopy, and no tissue samples could be collected. With the increasing number of genetically affected dogs detected by DNA-testing, we were able to continue the questionnaires and interviews initiated by Downs et al. The results strengthened the view that the disease in golden retrievers (generally referred to as GR\_PRA2) may indeed be similar to Bardet–Biedl syndrome (BBS) in humans. To diagnose BBS patients correctly is often challenging due to the heterogeneity of the disease with symptoms that vary even between individuals within families carrying identical mutations. In addition, the symptoms overlap with other diseases, such as Laurence–Moon syndrome, which may in fact, to some extent be considered a variation of the same condition [69]. Polydactyly or similar congenital digit abnormalities can be cues in the search for a correct diagnosis. In the absence

of digit malformations, the relatively early onset of retinal degeneration [6,7], can guide clinicians towards a BBS diagnosis. Both polydactyly and retinal degeneration are primary characteristics of BBS, as well as obesity, kidney abnormalities, cognitive impairment, hypogonadism in males, and genital abnormalities in females [7]. To date, 15 mutations in the human *TTC8* gene have been reported to cause BBS [2,17,70–76] (Table S6; Figure 7b). These mutations include six intronic splice-site variants, as well as three exonic splice-site variants (one missense, one nonsense, and one silent mutation). In addition, the identified mutations include three frameshift insertions/deletions, one complex variant, one nonsense variant, and one non-frameshift deletion in the coding sequence. Moreover, one *TTC8* splice-site variant has been associated with non-syndromic RP in humans [24]. Interestingly, one of the exonic splice-site variants (NM\_144596.3: c.1347G>C; p.Gln449His) has been reported to cause both non-syndromic RP [77] and BBS [70], suggesting that the genetic background of each patient may play a significant role in the development of BBS symptoms in humans. Similar to humans, the genetic background is also likely to result in the heterogeneous clinical signs in dogs (Table S2).

All the dogs in this study were diagnosed with PRA initially causing visual problems in dim light, later also under daylight conditions and eventually causing blindness. This is consistent with BBS, where rod-cone dystrophy is the most frequently observed clinical sign, diagnosed in about 93% of the patients with approximately three out of four becoming legally blind by the second to third decade of life [6,7]. FERGs from the middle-aged female dog (GR02) showed bilateral loss of rod function, as well as profoundly subnormal and delayed cone-amplitudes, supporting a diffuse rod-cone degeneration (Figure 3). This is similar to human patients with *TTC8* mutations, where loss of rod and cone function has been reported [24,72,77]. Early in the course of the disease, human BBS patients may have cone flicker responses with marked delays and near normal amplitudes that deteriorate over time [78]. *Ttc8*-knockout mouse models show both reduced rod and cone ERGs at an early age [79]. Interestingly, the mouse models also demonstrate that cone structure, function, and viability depend on the normal expression of *Ttc8.* Loss of the TTC8 protein results in shortened and disorganized photoreceptor outer segments already before retinal maturation, as early as post-natal day 10 in *Ttc8*-knockout mice [79]. The lack of rhodopsin staining, and the appearance of truncated cone outer segments observed in the retina of the affected female using fluorescence histochemistry (Figure 9), suggests that the deletion at the *TTC8* gene (*TTC8*delA) eventually results in a photoreceptor phenotype similar to the *Ttc8*-knockout mice. Taken together, the progressive reduction of cone photoreceptors and concomitant decline of cone ERGs may not only be secondary to degeneration of rods, as the TTC8 protein appears to be essential for normal cone function.

We were able to investigate the retina of the affected female using OCT, and considerable outer retinal thinning was observed (Figure 2), corroborating the results from ERG and postmortem findings. The abnormal appearance and thinning of the INL was judged as secondary to the photoreceptor degeneration. OCT-scans from the ventral (inferior) retina showed a more patchy degeneration with irregular RPE lining, also reported in some human patients with BBS [80], and clusters of nuclei in the ONL. Accumulation of lipofuscin, most clearly seen as increased autofluorescence in the non-tapetal fundus (Figure 1) has also been reported in human patients with BBS. However, the precise distribution of lipofuscin in the canine fundus is difficult to determine on cSLO because of the tapetal reflection in this species.

In addition to retinal degeneration, both the clinically investigated siblings, as well as three other affected dogs, were diagnosed with cataracts, a common secondary characteristic in human BBS [78]. In contrast to the posterior polar cataracts typically seen in human patients, both the affected male and the female developed cataracts that rather rapidly spread in the posterior, equatorial, and anterior cortices. The cataracts contributed to the dogs' visual impairment. The prevalence of posterior polar cataracts in human BBS has been reported to be similar in patients with different genotypes, although none of the patients studied was homozygous for pathogenic mutations in the *TTC8* gene [81]. We cannot exclude that the cataracts of the affected dogs were secondary to the retinal degeneration, because secondary cataracts are frequently seen in late-stage hereditary photoreceptor degenerations

(PRA) in dogs [82]. Toxic dialdehydes from decaying photoreceptors have been proposed to induce cataract formation [83] and the location of the cataractous changes in GR01 and GR02 in the more metabolically active equatorial and cortical regions of the lenses may also suggest a toxic effect on the lens. However, the cataracts in the golden retriever sib-pair were observed already at the initial examination when the degeneration of the retinae was judged as moderately advanced.

The affected male was also diagnosed with iridociliary cysts (Figure 4a). Iridociliary cysts have previously been reported in golden retrievers with glaucoma [84] and retinal detachments, but have not been associated with progressive retinal degenerations. Our observation of iridociliary cysts in the affected male is most likely a coincidental finding, but could possibly be a result of the deletion at the *TTC8* gene and a consequence of abnormal ciliary trafficking.

When Ansley et al. identified *TTC8* as a novel BBS gene in 2003 [2], situs inversus was observed in one of the patients with a *TTC8* splice-site mutation. This led to the important conclusion that BBS is a disease caused by ciliary dysfunction. The TTC8 protein is part of the BBSome, an octameric protein complex which functions in the exit of activated G protein–coupled receptors (GPCRs) from cilia [19,85,86]. Thus, the underlying cause of the syndrome is malfunction of primary cilia, an organelle emanating from the cell surface of most mammalian cell types during growth arrest [87]. Primary cilia function in cell signaling during development and in homeostasis, thus explaining the multitude of organs affected in each BBS patient and the repertoire of different symptoms experienced. Nodal cilia, a type of motile primary cilia transiently present during embryonal development, are crucial for breaking the left-right symmetry during the embryogenesis, and defects in this process may cause different forms of heterotaxy such as partial or total situs inversus, as well as congenital heart defects [88,89]. Manifestation of situs inversus and other forms of heterotaxy is rare in the general population (1/10,000). However, in a study of 368 participants in the Clinical Registry Investigating Bardet–Biedl Syndrome (CRIBBS), six patients (1.6%) were found to have disorders of asymmetry, suggesting a 170-fold increase compared to the general population [90]. Cardiovascular diseases and congenital heart defects have also been reported for human BBS patients [7,91]. Of the 19 patients with *TTC8* mutations (Table S6), 3 were diagnosed with hypertension or haemophilia [2,72,76]. We therefore investigated the cardiac status of the two affected golden retrievers. Both dogs showed normal cardiac morphology, dimensions and motion, and a normal left-right body axis. In the postmortem gross examination, a mild myxomatous valvular degeneration was seen in the heart tissue of the affected male. While the finding is common in dogs [92], it has recently been shown that defects of the primary cilia in the extra-cellular matrix can lead to myxomatous mitral valve disease in humans and mice [93].

Hypogonadism in males and genital abnormalities in females are considered primary characteristics of human BBS [6,7]. Hypogonadism was also common for the male patients with *TTC8* mutations (Table S6). The affected the males had a low total sperm count and a high number of abnormal sperm cells, mainly lacking flagellum or having other tail-defects, indicating infertility (Figure 5a). He did not display any interest in the female dog in heat at the time of sampling. The histopathological examination of the testes revealed only a few elongated spermatids (Figure 5b). Severe sperm defects, including tail abnormalities as well as a low total sperm count, have previously been reported in Hungarian Puli dogs with a loss-of-function mutation in the Bardet–Biedl syndrome 4 (*BBS4*) gene [94]. The importance of BBS-genes for the normal formation of flagella in spermatogenesis has further been shown in mouse models. Genetically modified mice for any of the genes *Bbs1* [95], *Bbs2* [96], *Bbs3* (*Arl6*) [97], *Bbs4* [98], *Bbs6* (*Mkks)* [99], or *Bbs7* [100], result in failure to form normal flagella. In a *Bbs8*-knockout mouse model, the sperm defects were not described, but the males were found to be infertile [79].

Interestingly, many of these mouse models, including the *Ttc8* knockout mouse [101], also showed partial or complete anosmia. This is a secondary characteristic experienced by many BBS patients in general [102], but it was not reported for the 19 human patients with *TTC8* mutations. This is in contrast with the present study, where most of the owners of the affected dogs (8/10) reported their dogs having poor or gradually worsening sense of smell (Table S2). This appear to be a characteristic that is shared between dogs and mice, but not with humans. Moreover, as was the case for the dogs, none of the mouse models exhibited polydactyly, and the phenotype of dogs may therefore resemble the mouse phenotype more closely than that of human patients with *TTC8* mutations.

Obesity was reported for 16 of the 19 human patients. Although we did not formally investigate the body condition score for more than two dogs, the majority of them appeared overweight or heavy according to the owners (Table S2). It should, however, be noted that we have not compared this to the general golden retriever population. Similar to the two golden retriever siblings in this study, dental anomalies were reported for one human patient with a *TTC8* splice-site mutation [76] and for 27% of the BBS patients in general [7].

The histopathological changes in the right kidney of the affected male were most similar to a renal dysplasia-like lesion (Figure 6). The female sibling had mild chronic bilateral glomerulonephritis, but such mild chronic inflammatory changes are an unspecific finding, which probably did not affect the animal clinically. In total, 4 of the 10 dogs included in this study had renal problems (Table S2), suggesting that, as in humans, mutations in the *TTC8* gene may result in variable renal phenotypes. Among the 10 golden retrievers included in this study (Table S2), 1 dog was reported to be euthanized due to kidney failure. Other reasons for euthanasia included gastrointestinal problems (2 dogs), as well as problems related to sensory deprivation (visual and olfactory impairments: 2 dogs) and neoplastic disease (1 dog). Interviews indicated a rather low average lifespan of the affected golden retrievers: 7 years and 8 months for the nine dogs in the questionnaires (one of the ten dogs is still alive as of preparing this manuscript, 5 years of age). Data from 1995–2002 show that almost 90% of golden retrievers survive to 8 years of age and more than 80% to the age of 10 [103]. The 1 bp deletion in the canine *TTC8* gene had detrimental effects on the dogs' health and longevity. In human BBS patients, kidney failure is a frequent cause of death and the average survival of the patients overall is substantially reduced [104]. Taken together, the deletion in the canine *TTC8* gene is associated with additional clinical features apart from PRA in the affected dogs, such as obesity, renal and genital anomalies, anosmia, short stature, and dental anomalies that are similar to human BBS.

The major isoforms of the TTC8 protein are predicted to consist of approximately 500 amino acid residues, and the exon-intron structure is well conserved between human, dog, and mouse. Two canine retinal transcripts were previously predicted to exist based on Sanger sequencing of cDNA by Downs et al (2014) [1], one with and one without the retinal specific exon termed 2a. The transcript with the retinal specific exon was originally reported in humans [24]. It was shown that an in-frame splice-site mutation leading to the skipping of exon 2a is sufficient to cause a non-syndromic retinal degeneration in humans. In mice, the retinal specific exon 2a was found to be exclusively expressed in the ocular tissue, having the highest expression in the outer segments (OS) of both rod and cone photoreceptor cells [24,105]. In addition to retinal photoreceptor cells, exon 2a has also been found expressed in pinealocytes of rat pineal gland, and it has been suggested that these two cell types have evolved from a common precursor photodetector cell [106]. In our study, we used Oxford Nanopore Sequencing Technology for long-read sequencing, capable of reading through full-length transcripts and therefore enhancing the detection of different splice sites. Notably, the data only gave support for three different transcripts expressed in the retinal tissue (Figure 7), while the Ensembl annotation (release 100) predicts the existence of seven canine transcripts, of which only one (ENSCAFT00000027700.5 corresponding to tr2) was identified in this study. The six additional Ensembl canine transcripts are either skipping exon(s): E5 (ENSCAFT00000043691.2), E6 (ENSCAFT00000050179.3), E5-6 (ENSCAFT00000093101.1), E2a-10 (ENSCAFT00000070341.1), E1-2a (ENSCAFT00000086744.1), or E1-5 (ENSCAFT00000087482.1). We found no evidence for the presence of these other transcripts in our data from the canine retina.

As expected, the tr2 transcript, which includes the retinal specific exon 2a, was found to be the most highly expressed of the three identified transcripts. We also identified a short transcript (tr1) without the retinal specific exon 2a, with low expression in the dog retina. In mice, despite having a lower expression than the retinal specific transcript, the short transcript lacking the exon 2a is expressed throughout the retinal cell layers and shows the highest expression in the RPE [24]. The retinal

samples investigated in this study did not exclusively include photoreceptor cells, but also consisted of other neural cells and RPE cells (Figure 8), which likely express the short transcript tr1. In addition to these two previously identified canine retinal transcripts, we also found evidence for a third transcript (tr3) including both exon 2a as well as a short exon with 22 nucleotides, here denoted 1b. A transcript similar to tr3 is also reported in GenBank (XM\_014115661.2), likely annotated by an ab initio prediction detecting an open reading frame of sufficient length, but likely without supporting alignment (The NCBI Eukaryotic Genome Annotation Pipeline; A translation initiation site (TIS) with an optimal Kozak consensus sequence was identified in exon 1, but in tr3 this would extend an ORF of 147 nucleotides over the exons 1b and 2a, and reach a stop codon in exon 2 (Figure 7). This short transcript would likely be degraded by nonsense-mediated decay (NMD) [107]. If an alternative TIS was used in exon 2, tr3 would be in-frame with the other transcripts and produce a protein product of 455 amino acid residues featuring the complete TPR motif, which is involved in protein-protein interaction [108]. It is yet to be defined if this transcript is translated into a protein product, and if its expression is tissue specific.

The expression levels for all three *TTC8* transcripts found in this study were markedly reduced in the affected female dog. This was seen in the quantification of the cDNA sequencing data (Figure 8a–c), but also by the lack of read coverage of transcripts mapping to the *TTC8* locus in the dog reference genome sequence when manually inspecting the data in IGV. The highest read depth in the locus was 13× (TPM = 16.8), while the read coverage for the two unaffected dogs (250–500×, TPM = 60.6 and 99.8, respectively) was clearly higher than in the affected female. The reads detected in the affected dog appeared to be only partly spliced and none of the reads reached full-length over all exons. The low levels of *TTC8* transcripts in the affected dog was further supported by RT-qPCR (Figure 8d). The data suggested that the deletion in exon 7 (*TTC8*delA) and the subsequent premature stop codon in exon 8 lead to degradation of the *TTC8* transcripts in the affected individual. To estimate the cell type constitution of the samples, we also studied the expression levels of retinal marker genes in the affected female compared to the two unaffected female dogs. As suggested by the histological sectioning of the affected male retina (Figure S1) and the fluorescence histochemistry of the affected female (Figure 9), the sequencing results showed a drastically lower expression of rod photoreceptor cell-specific genes, and cone photoreceptor cell marker expression was only one third compared to the levels detected in the two other dogs (Figure 8a–c, Table S4). The same pattern was observed in the RT-qPCR, where we quantified the relative expression of *RHO* and *OPN1LW* genes (Figure 8d). The retinal degeneration of both of these dogs (GR01 and GR02) was therefore likely advanced. BBSome deficiency has been shown to cause defects in the transport of phototransduction proteins between the inner and outer segments of the photoreceptors and, ultimately, these defects lead to photoreceptor cell death [96,109]. In addition to aberrant photoreceptor marker gene expression, our data also suggested that the expression of macroglial genes (Müller glia and astrocytes) was higher in the affected dog (Figure 8), most notably clusterin (*CLU*) and glial fibrillary acidic protein (*GFAP*) expression (Table S4), both of which are known to be upregulated under retinal stress and retinal degeneration [110–113].

#### **5. Conclusions**

The *TTC8* gene encodes for one of the proteins forming the BBSome, and has in humans been implicated in Bardet–Biedl syndrome (BBS). Long-read cDNA sequencing of non-affected dogs suggested the expression of three retinal *TTC8* transcripts and that the 1 bp deletion is a loss-of-function mutation. Golden retriever dogs homozygous for the deletion develop an autosomal recessive form of RP-like retinal degeneration (PRA), but it has hitherto been unclear if the affected dogs develop a non-syndromic PRA or a syndromic ciliopathy similar to human BBS. In addition to PRA, we have shown that the loss-of-function mutation indeed causes additional clinical features, such as obesity, renal and genital anomalies, anosmia, short stature, and dental anomalies. We therefore conclude that the deletion can result in a canine form of BBS. As in humans, BBS in dogs appear to be a heterogeneous disorder with variable severity of clinical and morphological signs. A canine model for BBS may be of

importance for novel therapeutic management of human patients. Canine models have successfully been used to establish protocols for gene therapy of other inherited retinal diseases, and a *TTC8*-dog model could potentially be developed to restore vision and improve the quality of life for BBS patients.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/9/1090/s1, Figure S1: Retinal changes associated with *TTC8*delA, Table S1: List of dogs included for OCT and/or collection of tissue samples, Table S2: Primary and secondary characteristics, Table S3: Macroscopic and microscopic examination, Table S4: Marker gene expression, Table S5: Primers, Table S6: Human *TTC8* mutations.

**Author Contributions:** Conceptualization, S.M. and T.F.B.; methodology, S.M., M.H., A.D.v.d.H., E.A., E.T., B.S.H., J.H., I.L., C.M., F.H., G.A., B.E. and T.F.B.; software, S.M.; validation, S.M. and A.D.v.d.H.; formal analysis, S.M., M.H., A.D.v.d.H., E.A., E.T., B.S.H., B.E. and T.F.B.; investigation, S.M., M.H., A.D.v.d.H., E.A., E.T., B.S.H., J.H., I.L., C.M., F.H., G.A., B.E. and T.F.B.; resources, S.M., M.H., E.A., E.T., B.S.H., J.H., I.L., F.H., G.A., B.E. and T.F.B.; data curation, S.M., M.H., E.A., E.T., B.S.H., J.H., I.L., F.H., G.A., B.E. and T.F.B.; writing—original draft preparation, S.M. and T.F.B.; writing—review and editing, S.M., M.H., A.D., V.D.v.d.H., E.A., E.T., B.S.H., J.H., I.L., C.M., F.H., G.A., B.E. and T.F.B.; visualization, S.M., M.H., A.D.v.d.H., E.A., E.T., B.S.H., F.H., B.E. and T.F.B.; supervision, C.M., F.H., G.A., B.E. and T.B.; project administration, G.A. and T.F.B.; funding acquisition, G.A. and T.F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by "The Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning" (FORMAS; http://formas.se), grant number: 221-2014-1005) and by the "Agria and SKK Research Fund" (https://www.skk.se/sv/Agria-SKK-Forskningsfond/, grant numbers: P2015-0012, N2016-0015, N2017-0016, N2018-0013.

**Acknowledgments:** The authors would like to acknowledge the support of the National Genomics Infrastructure (NGI)/Uppsala Genome Center and UPPMAX for providing computational infrastructure and the Swedish Bioinformatics Infrastructure Sweden at SciLifeLab for bioinformatics advice. We would also like to acknowledge Sofia Ryberg at the Department of Clinical Sciences, SLU for assisting in the ophthalmic examinations and support of the dedicated dog owners who allowed their dogs to take part in this study.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **X-Linked Duchenne-Type Muscular Dystrophy in Jack Russell Terrier Associated with a Partial Deletion of the Canine** *DMD* **Gene**

**Barbara Brunetti 1,\*, Luisa V. Muscatello 1, Anna Letko 2, Valentina Papa 3, Giovanna Cenacchi 3, Marco Grillini 4, Leonardo Murgiano 2,5, Vidhya Jagannathan <sup>2</sup> and Cord Drögemüller <sup>2</sup>**


Received: 3 September 2020; Accepted: 5 October 2020; Published: 8 October 2020

**Abstract:** A 9-month old male Jack Russell Terrier started showing paraparesis of the hindlimbs after a walk. Hospitalized, the dog went into cardiac arrest, and later died. Necroscopic examination revealed a severe thickness of the diaphragm, esophagus, and base of the tongue, leading to the diagnosis of muscular dystrophy. The histology confirmed the marked size variation, regeneration, and fibrosis replacement of the skeletal muscle fibers. Immunohistochemistry demonstrated the absence of dystrophin confirming the diagnosis. Transmission electron microscopy showed disarrangement of skeletal muscle fibers. Finally, whole-genome sequencing identified a ~368kb deletion spanning 19 exons of the canine dystrophin (*DMD*) gene. This pathogenic loss-of-function variant most likely explains the observed disease phenotype. The X-chromosomal variant was absent in seven controls of the same breed. Most likely, this partial deletion of the *DMD* gene was either transmitted on the maternal path within the family of the affected dog or arose de novo. This study revealed a spontaneous partial deletion in *DMD* gene in a Jack Russell Terrier showing a Duchenne-type muscular dystrophy due to non-functional dystrophin.

**Keywords:** canine; dystrophinopathy; Duchenne; immunohistochemistry; precision medicine

#### **1. Introduction**

Duchenne and Becker muscular dystrophies are X-linked recessive disorders, therefore these forms occur predominantly in males. They are caused by genetic variants in the dystrophin gene (*DMD*), the largest gene in the human genome, and as such are called dystrophinopathies [1]. The *DMD* gene spanning 2.4 Mb on the X chromosome encodes 79 exons for a 427-kD protein, called dystrophin, a rod-shaped protein located on the inner face of the plasma membrane of all types of myofibers and anchors cytoskeletal F-actin to the extracellular matrix protein laminin [2,3]. Dystrophin protein has four main functional domains. The N-terminus is the cysteine-rich actin-binding domain, while the carboxy-terminal domain interacts with β-dystroglycan as well as with dystrobrevin and syntrophin. The central rod domain comprises the majority of the mass of the dystrophin molecule, forming a

flexible, rod-shaped structure [4]. In humans, a causative genetic variant in *DMD* was found in 96% of Duchenne muscular dystrophy (DMD) cases and 82% of Becker muscular dystrophy (BMD) cases. Around one third of pathogenic variants in *DMD* are spontaneously occurring de novo mutations in the affected male patients [1]. The most common genetic variants within *DMD* are large deletions (approximately 70%) or duplications (10–14%) often encompassing numerous exons of the gene [1].

Mutations may occur throughout the 79 exons of the *DMD* gene but concentrate in major (exons 45–53) and minor (exons 2–20) hotspot areas. According to Leiden's database, ~40% of *DMD* gene variants are deletions of a mean length of 6.5 exons, with exon 47 being most commonly affected [5]. These deletions tend to predominate in one of two hotspots, namely, the central rod domain around exons 44–53 (~80%) and, to a lesser extent (~20%), at the 5 -end of the gene [6,7]. Duplications occur most frequently in the region of exon 20 [5]. Dystrophin is important for the maintenance of the structural integrity of muscle fibers during contraction. Both Duchenne- and Becker-types of muscular dystrophies have similar signs and symptoms and are due to different mutations in the same gene. In human Duchenne muscular dystrophy, variants affecting the *DMD* gene result in a severely truncated, non-functional dystrophin, while in Becker muscular dystrophy, the mutations result in a truncated semi-functional dystrophin protein [8]. The two conditions differ in their severity, age of onset, and rate of progression. The signs of BMD are usually milder and more variable [9].

The diagnostic evaluation of muscular dystrophy often includes a muscle biopsy to demonstrate fibrosis, muscle fiber size variation, internalization of muscle nuclei, and more importantly, the abnormal expression of dystrophin by immunostaining or Western blot analysis. Fibro-fatty replacement, inflammation, and degenerative fibers are commonly described histological features [10]. Similar muscular dystrophies, homologous to human DMD and BMD, also occur in dogs (OMIA 001081-9615), cats (OMIA 001081-9685), pigs (OMIA 001081-9823), and mice [9]. Thus far, 12 different causative *DMD* variants for DMD have been described in several dog breeds (OMIA 001081-9615). For instance, the Golden Retriever muscular dystrophy model is widely used as a useful animal model, because these severely affected dogs show a phenotype strongly resembling human DMD on a clinical and histopathological level [10].

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

The dog in this study was privately owned, and samples were collected during the necropsy requested from the owner because the dog died unexpectedly.

#### *2.2. Clinical History and Necropsy Request*

A 9-month old male Jack Russell Terrier started showing paraparesis of the hindlimbs after a walk, but he was still able to stand and walk; during the next day, dyspnea developed and the dog was immediately taken to the veterinarian. At arrival in the practice, he went into cardiac arrest, and although resuscitation was tried, he finally died. The owner reported that the dog was originally imported from Poland, and he confirmed that the dog had always been healthy and regularly vaccinated. As a result of the sudden and unexpected death of the subject, the owner requested a necropsy to identify the causes of death.

#### *2.3. Histopathology*

During necropsy, tissue samples for histological evaluation were fixed in 10% buffered formalin. Particular attention was paid to sampling the various muscles both affected by the pathology and those that were apparently normal. Diaphragm, esophagus, tongue, thigh (quadriceps femoral), and heart muscles were sampled. The quadriceps femoris muscle macroscopically looked normal.

Tissue was fixed in 10% formalin for 24 h at room temperature, then it was embedded in paraffin and cut at a thickness of 4 microns. The sections were stained with hematoxylin and eosin (H&E) and Masson's trichrome stain.

#### *2.4. Immunohistochemistry*

Immunohistochemistry was performed using three antibodies against-dystrophin: C-terminus simultaneously specific for anti-human, -dog, and -mouse (corresponding to the 3 end of the dystrophin gene); anti-human rod domain and anti-human N-terminus (corresponding to the 5 end of the dystrophin gene). See specific technical data in Table S1 of supplementary material.

#### *2.5. Transmission Electron Microscopy*

Samples were fixed in 1% OsO4 in cacodylate buffer, dehydrated in graded ethanol, and embedded in Araldite. Thin sections, stained with uranyl acetate and lead citrate, were examined using a Philips TEM CM100 Transmission Electron Microscope.

#### *2.6. Whole-Genome Sequencing (WGS)*

In order to investigate the underlying molecular basis, whole-genome sequencing (WGS), using genomic DNA isolated from the blood sample of the affected dog, was performed as described before [11]. Data corresponding to roughly 29× coverage of the genome were collected on an Illumina NovaSeq6000 instrument (2 × 150 bp). Read mapping, re-alignment, and variant calling were carried out as previously described [11] with respect to the CanFam3.1 genome reference assembly and the NCBI annotation release 105. The *DMD* gene representing the functional candidate was visually inspected for structural variants using Integrative Genomics Viewer [12]. WGS of the affected dog is available at the European Nucleotide Archive [13] sample accession SAMEA6249497. Numbering within the canine *DMD* gene as reported in the paper refers to the mRNA accession no. NM\_001003343.1 and the protein accession no. NP\_001003343.1.

#### **3. Results**

#### *3.1. Necropsy Examination*

During necropsy, on external examination, the muscle conformation was normal without signs of muscle hypotrophy or hypertrophy and in a good state of nutrition. Interestingly, the diaphragm muscles were diffusely severely thickened, contracted with severe fibrosis, and the thickness of the muscle was about 1 cm (Figure 1A). In addition, the muscle of the base of the tongue and the distal part of esophagus were severely thickened (Figure 1B). The heart was increased in volume (cardiomegaly). As seen in the section, the right ventricular lumen was markedly dilated (dilated cardiomyopathy). The left and right atrioventricular valves were moderately swollen and edematous (valvular endocardiosis). Some sero-hemorrhagic fluid (around 10 mL) was found in the thoracic cavity with a severe acute pulmonary edema; the lungs showed multifocal petechiae and hemorrhages; whitish foam was found in the trachea and large bronchi. The sero-hemorrhagic fluid and the lung's petechiae and hemorrhages could be due to resuscitation procedures.

**Figure 1.** (**A**) Macroscopic image of diaphragm: it was severely thickened until 1 cm in section as well as the distal part of esophagus (4–5 mm thick) (**B**).

#### *3.2. Histopathology*

All the muscles examined (diaphragm, esophagus, tongue, heart, thigh muscle) showed marked alterations. In the longitudinal sections, hematoxylin-eosin (H&E) rows of central nuclei in the fibers (regeneration) and fiber splitting (Figure 2A) were evident, and a moderate fibrosis was appreciated and confirmed with Masson's trichrome stain. Cross-sections revealed marked variations in the fibers' diameter with atrophy in some and hypertrophy in others (Figure 2B). Multifocal myofibers were fragmented and hypereosinophilic with a lack of a transversal band (coagulative necrosis) and internal nuclei and were multifocally infiltrated by lymphocytes and macrophages. There was also mild to moderate adipose tissue replacement, and fiber mineralization was rare. In the myocardium, there was slight fibrosis, whereas the lymphocytic and macrophage interstitial infiltrate was multifocal and moderate. A large focal area of coagulative necrosis was present.

**Figure 2.** Longitudinal section of diaphragm stained with hematoxylin-eosin (H&E) (**A**), showing the rows of nuclei in the fibers (regeneration) and fiber splitting. (**B**). Cross-section of diaphragm stained with H&E showing marked atrophy in some muscle fibers and hypertrophy in others, and fiber splitting. A diffuse and mild fibrosis was also present. Bar = 100 micron.

#### *3.3. Immunohistochemistry*

All three antibodies tested for dystrophin gave a strong positivity and membrane positivity on the diaphragmatic muscle of a dog of the same breed used as a positive control (Figure 3A–C). On the contrary, all examined sections of the diaphragm, esophagus, tongue, and thigh muscle of the affected dog did not show anti-C-terminal (Figure 3D) and N-terminal (Figure 3F) antibody positivity, with the exception of some multifocal revertant fibers. There was weak and multifocal positivity to the anti-Rod antibody (Figure 3E).

**Figure 3.** Control dog. Immunostains for the C-terminus (**A**), rod domain (**B**), and N-terminus (**C**) of the dystrophin in the diaphragm of a normal dog showing diffuse and strong sarcolemmal staining (positive control, transverse sections, objective 20×). In the affected dog, there was a total loss of positivity using both the anti-C-terminal (**D**) and N-terminal antibodies (**F**), and a very mild and multifocal signal with the anti-Rod antibody (**E**) (transverse sections, objective 20×).

#### *3.4. Transmission Electron Microscopy*

Transmission electron microscopy (TEM) analysis performed on the diaphragm and tongue muscle showed mild myopathic changes as Z-band streaming (Figure 4A) and myofibrillar disarray (Figure 4B) with mitochondria hyperplasia (Figure 4C).

**Figure 4.** Electron microscopy of diaphragm muscle (**A**,**B**) and tongue (**C**): mild myopathic changes were seen as Z-line streaming (z), myofibrillar disarray (d), and mitochondrial hyperplasia (m): mitochondria seemed to be increased in size compared to the control (insert). 10,500×.

#### *3.5. Whole-Genome Sequencing*

In the region of the *DMD* gene, a large hemizygous structural variant was detected in the genome of the affected dog (Figure 5). The variant encompasses a 367,633 bp region upstream of exon 3, up to part of intron 21. The deletion truncates the coding sequence of 19 of the 79 *DMD* exons annotated in the canine reference genome (exon 3 to exon 21). The formal variant designation is CFAX

g.[27,615,280\_27,982,912del]. No further single nucleotide variants affecting the coding regions of this gene were detected. Furthermore, the variant was absent from seven (five males, two females) publicly available WGS of Jack Russell Terriers as well as 576 dogs of 125 breeds and eight wolves [11].

**Figure 5.** Large hemizygous partial deletion (CFAX g.[27,615,280\_27,982,912del]) of the canine *DMD* gene in the Duchenne muscular dystrophy (DMD)-affected dog. The deletion spans exons 3–21 (highlighted by red box) and part of the flanking introns. The lower part shows an integrative genomics viewer screenshot of the region of interest from the whole-genome sequence data. Note the absence of mapped reads (shown in grey) within the region of 19 coding exons (red arrow). The light grey lines indicate read pairs, mapping across the breakpoints.

#### **4. Discussion**

This report details a dystrophin-deficient 9-month-old Jack Russell Terrier with an abrupt onset of clinical signs that resulted in death the next day. The Polish owner reported that the affected dog had never shown any health problems before, and apart from the known Polish geographical origin, there was no pedigree information available. The juvenile age of death of this subject was in line with other cases of canine DMD reported previously, with a range from a few weeks to a few months of life [5,14,15]. The great majority of Golden Retriever Muscular Dystrophy (GRMD) dogs do not survive beyond age two [16]. Only some rare cases of similarly affected dogs of up to 5.6 years of age were reported before [17]. The clinical signs probably were not observed by the owner, but in any case, considering the conditions of the diaphragm, esophagus, tongue, and heart, we may suppose that it took place. It is obviously speculative, but we can assume the occurrence of heart failure, respiratory failure, and dysphagia.

It is reported that dystrophin deficiency results in progressive gross muscle atrophy in most breeds of dogs, but it causes a marked muscular hypertrophy in the mice, cats, and in dogs such as the Rat Terrier. Muscle fiber hypertrophy occurs, especially in early stages of the disorders, but extensive fiber necrosis in most cases leads to a general muscle atrophy. At this time, there is no explanation for this phenomenon, although it would appear that muscle hypertrophy is more apparent in animals of small stature [9].

Immunostaining in dystrophinopathies routinely involves the use of three antibodies targeting the different region of dystrophin, namely, the C- and N-terminus, as well as the central rod domain. In these cases, we tested all the three antibodies on formalin-fixed paraffin-embedded (FFPE) tissue of the case and one control dog. The results obtained in the affected dog were very similar to those observed in DMD-affected human patients [18], with total loss of positivity using both the anti-C-terminal and N-terminal antibodies, and a very mild and multifocal signal with the anti-Rod antibody. These results based on the use of these antibodies working on FFPE tissue are important in

veterinary medicine, as it is not always possible to obtain fresh muscle biopsies and to then freeze them to send to suitable laboratories. Furthermore, the use of these antibodies allows the diagnosis of muscular dystrophy to be verified on the sample available to the pathologist. It is also possible to work retrospectively on archived samples. The key tests performed on the biopsy sample for muscular dystrophies were immunohistochemistry and immunoblotting for dystrophin, and therefore, the presented immunohistochemistry findings in the studied Jack Russell Terrier are highly suggestive of the presence of a Duchenne-type disease phenotype.

The first identified case of canine muscular dystrophy was in a Golden Retriever in 1958. Since then, various muscular dystrophy phenotypes have been reported in at least 15 other breeds [19]. Therefore, the obtained genetic results confirmed the previous clinical and pathological suspicion of a Duchenne-type muscular dystrophy in the affected Jack Russell Terrier. Different but similar gross deletions have been reported before in the Duchenne type of muscular dystrophy in dogs encompassing the entire *DMD* gene (OMIA 001081-9615) [2,15]; or, similarly to the findings in the herein presented dog, affecting a significant number of the coding exons, as in an affected Tibetan Terrier showing a deletion encompassing *DMD* exons 8 to 29 [20]. Other previously reported *DMD*-associated loss-of-function variants causing canine Duchenne-like muscular dystrophies include structural variants such as smaller sized deletions or gross insertions, as well as variants affecting single or smaller numbers of nucleotides such as splicing variants, inversions, or nonsense (stop-gain) variants (OMIA 001081-9615).

#### **5. Conclusions**

The herein reported partial deletion of approximately a fifth of the canine *DMD* gene leads to a true null allele of the *DMD* gene, with the lack of protein expression being experimentally confirmed. The identified deletion variant might have occurred de novo during meiosis of the maternal germ cells and subsequently passed to the offspring as a consequence of low-level mosaicism in the mother, or it may have been transmitted as an X-linked recessive variant on the maternal path within the family of the affected dog. No evidence for further cases was obtained; we can, therefore, speculate that the mutation might also have occurred post-zygotically during early embryonic development of the affected dog explaining this single case. Taken together, this is the first report of a *DMD*-associated Duchenne-like muscular dystrophy in the Jack Russell Terrier breed. This study provides an example of a pathogenic disease-causing variant underlying a sporadic syndrome observed in dogs. Precision medicine using WGS has been proven to be suitable in recent years to help diagnose rare diseases in routine diagnoses in veterinary medicine: to both confirm genetic etiology and to obtain an insight into the molecular mechanisms involved. Nonetheless, the Golden Retriever muscular dystrophy continues to be the best-documented of the canine dystrophinopathies [21]. Thus far, there have been reports of 12 different canine breed-specific DMD-associated variants in the *DMD* gene (OMIA 001081-9615). The herein described phenotype of the affected Jack Russell terrier was most likely caused by the identified large genomic ~368kb deletion, spanning 19 coding exons of the *DMD* gene. This partial deletion truncates 18% of the coding sequence from exon 3 to 21 from the 79 exons annotated in the canine reference genome.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/10/1175/s1, Table S1. The table shows the technical specifications of the 3 antibodies used to highlight the various parts of dystrophin.

**Author Contributions:** Conceptualization, B.B. and L.V.M.; methodology, B.B, L.V.M., A.L., V.P, G.C., and M.G. ; software, V.J. and L.M.; validation, V.P., G.C., and M.G.; data curation, B.B.; writing—original draft preparation, B.B.; writing—review and editing, B.B., L.V.M., A.L., and C.D.; visualization, B.B., L.V.M., A.L., and V.P.; supervision, C.D.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The Next Generation Sequencing Platform of the University of Bern is acknowledged for performing the whole-genome re-sequencing experiments and the Interfaculty Bioinformatics Unit of the University of Bern for providing high performance computing infrastructure.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Case Report*

### *SLC19A3* **Loss-of-Function Variant in Yorkshire Terriers with Leigh-Like Subacute Necrotizing Encephalopathy**

**Michaela Drögemüller 1,**†**, Anna Letko 1,**†**, Kaspar Matiasek 2,**†**, Vidhya Jagannathan 1, Daniele Corlazzoli 3, Marco Rosati 2, Konrad Jurina 4, Susanne Medl 5, Thomas Gödde 6, Stefan Rupp 7, Andrea Fischer 8, Alejandro Luján Feliu-Pascual <sup>9</sup> and Cord Drögemüller 1,\***


Received: 17 September 2020; Accepted: 15 October 2020; Published: 16 October 2020

**Abstract:** Sporadic occurrence of juvenile-onset necrotizing encephalopathy (SNE) has been previously reported in Yorkshire terriers. However, so far, no causative genetic variant has been found for this breed-specific form of suspected mitochondrial encephalomyopathy. Affected dogs showed gait abnormalities, central visual defects, and/or seizures. Histopathological analysis revealed the presence of major characteristics of human Leigh syndrome and SNE in Alaskan huskies. The aim of this study was to characterize the genetic etiology of SNE-affected purebred Yorkshire terriers. After SNP genotyping and subsequent homozygosity mapping, we identified a single loss-of-function variant by whole-genome sequencing in the canine *SLC19A3* gene situated in a 1.7 Mb region of homozygosity on chromosome 25. All ten cases were homozygous carriers of a mutant allele, an indel variant in exon 2, that is predicted to lead to a frameshift and to truncate about 86% of the wild type coding sequence. This study reports a most likely pathogenic variant in *SLC19A3* causing a form of SNE in Yorkshire terriers and enables selection against this fatal neurodegenerative recessive disorder. This is the second report of a pathogenic alteration of the *SLC19A3* gene in dogs with SNE.

**Keywords:** *Canis familiaris*; whole-genome sequencing; rare disease; precision medicine; neurometabolic disorder

#### **1. Introduction**

Subacute necrotizing encephalomyelopathy (SNE), also termed Leigh syndrome (LS; OMIM 256000) represents a devastating neurodegenerative disorder in people, characterized by a wide variety of clinical signs, ranging from severe neurologic problems to a near absence of abnormalities with the central nervous system most frequently affected [1]. Originally, Archibald Denis Leigh, a British neuropsychiatrist described the condition in 1951 [2]. SNE is characterized by focal and bilaterally symmetrical, necrotic lesions involving the thalamus, brainstem, and posterior columns of the spinal cord [3]. In SNE, various mutations in mitochondrial respiratory chain complexes lead to the disruption of ATP synthesis resulting in the characteristic pathology of SNE [4]. Mitochondrial encephalomyelopathies, such as SNE or LS, represent rare inherited neurometabolic disorders showing considerable genetic heterogeneity and associated pathogenic variants affecting over 85 different genes of the mitochondrial or nuclear genome [3]. Therefore, they represent mitochondrial disorders with the largest genetic heterogeneity [1].

As human SNE is rare and heterogeneous, studying domestic animal species showing resembling conditions might add to the understanding of such a complex group of disorders. Rare forms of SNE were described e.g., in cattle [5–7] and dogs (OMIA 001097-9615). The first report of this disorder was described in Alaskan huskies [8,9] and subsequently, a similar form of SNE was reported in Yorkshire terriers [10] and American Staffordshire bull terriers [11]. Neuropathologically, SNE in Yorkshire terriers is nearly identical to the Alaskan husky form and very similar to human Leigh syndrome [10]. An initial genetic investigation of SNE-affected Yorkshire terriers revealed no indication for disease-causing variants in the mitochondrial genome [10], whereas more recently in Alaskan huskies the pathogenesis of recessively inherited SNE was unraveled [12,13]. This breed-specific fatal brain disorder in Alaskan huskies is associated with a deleterious loss-of-function variant in *SLC19A3* encoding for a thiamine transporter 2 (THTR2) with a predominately central nervous system (CNS) distribution [12,13]. The *SLC19A3* gene product controls the uptake of thiamine in the CNS via expression of the thiamine transporter protein THTR2. Pathogenic variants are associated with thiamine metabolism dysfunction syndrome-2 in people (THMD2; OMIM 607483), also known as biotin-responsive basal ganglia disease (BBGD) or thiamine-responsive encephalopathy [14]. This *SLC19A3*–related condition is an autosomal recessive disorder with childhood-onset that presents as a subacute encephalopathy and progresses to severe cogwheel rigidity, dystonia, quadriparesis, and eventually death if left untreated (OMIM 606152). The *SLC19A3*–related SNE of Alaskan huskies was proposed as a possible large animal model that may allow prospective investigations into the mechanisms of *SLC19A3*-related syndromes and the potential role of thiamine and/or biotin as a therapeutic strategy [12,13].

To elucidate the disease mechanism underlying monogenic autosomal recessive inherited SNE in Yorkshire terriers, we applied homozygosity mapping and whole-genome sequencing revealing a most likely pathogenic variant in the canine *SLC19A3* gene.

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

All animal experiments were performed according to the local regulations. All animals in this study were examined with the consent of their owners. Sample collection was approved by the Cantonal Committee for Animal Experiments (Canton of Bern; permit 71/19).

#### *2.2. Animals*

In total, 172 blood samples of Yorkshire terriers were collected. Ten dogs were diagnosed with Leigh-like subacute necrotizing encephalopathy (SNE) according to Baiker et al. [10]. These affected dogs were unrelated, apart from two full siblings for which their sire and dam (obligate carriers) as well as a single normal littermate were also available. The remaining 159 dogs represented unrelated purebred controls. Genomic DNA was isolated from EDTA blood samples using the Maxwell RSC Whole Blood DNA kit (Promega, Dübendorf, Switzerland).

#### *2.3. Single Nucleotide Polymorphism Array Genotyping*

Four selected SNE-affected Yorkshire terriers were genotyped on Illumina CanineHD BeadChip array (Illumina, San Diego, CA, USA). PLINK v1.9 [15] was used to perform the quality control filtering steps of the obtained genotyping data and the subsequent homozygosity analysis. Single nucleotide polymorphisms (SNP) with a call rate <90% were removed leaving 167,185 markers. All individuals had call rates >90%. Homozygosity analysis was carried out with PLINK v1.9 [15] to determine intervals of extended homozygous regions with alleles shared by all four affected dogs.

#### *2.4. Whole-Genome Sequencing*

Whole-genome sequence (WGS) data of a single affected dog was obtained at 19.7× coverage in order to identify the causative variant for SNE. The sequence data analysis and calling of single nucleotide variants and small indels (SNVs) including the prediction of functional effects were described before [16]. The dog reference genome assembly CanFam3.1 and NCBI annotation release 105 was used. Additionally, a publicly available control genomes cohort of 720 dogs from 130 various breeds, and 9 wolves [16] was used to filter variants private in the sequenced SNE-affected dog; this also included 60 unrelated Yorkshire terriers (Supplementary Table S1). The Integrative genomics viewer (IGV) software [17] was used for visual inspection and screening for structural variants in the associated regions.

#### *2.5. Sanger Sequencing and Targeted Genotyping*

Polymerase chain reaction (PCR) and Sanger sequencing were used to validate and characterize the *SLC19A3* indel variant (XM\_022409850.1:c.205\_210delins35) identified from whole-genome sequencing. PCR primers were designed using primer 3 [18]. PCR products from genomic DNA were amplified using AmpliTaqGold360 MasterMix (Thermo Fisher Scientific, Waltham, MA, USA) and the purified PCR amplicons were directly sequenced on an ABI3730 capillary sequencer (Thermo Fisher Scientific) using the following primers: GGCAGTCACCATCCCATAGA (forward) and GATATTGGGCAAGCCACCTA (reverse) generating 309 bp products. The sequence data were analyzed with Sequencher 5.1 software (GeneCodes, Ann Arbor, MI, USA). Diagnostic genotyping was performed by fragment length analysis using a different forward primer (ATCCCTTGCAGGATGATGAC) to produce amplicons of 218 bp or 247 bp representing the wild type or variant allele, respectively. The 29 bp size difference was visualized on a Fragment Analyzer capillary gel electrophoresis instrument (Advanced Analytical Technologies, Ames, IA, USA).

#### *2.6. Availability of Data and Material*

The whole-genome data of an SNE-affected Yorkshire terrier are freely available at the European Nucleotide Archive (ENA) under sample accession number SAMEA3928145. All accession numbers of the used control genomes are available in Supplementary Table S1.

All genome positions are reported with respect to dog reference genome assembly CanFam3.1 and NCBI annotation release 105. All references to the canine *SLC19A3* gene correspond to the accessions NC\_006607.3 (NCBI accession), XM\_022409850.1 (mRNA), and XP\_022265558.1 (protein).

#### **3. Results**

#### *3.1. Homozygosity Analysis*

Based on the clinicopathological diagnosis of Leigh-like subacute necrotizing encephalopathy (SNE) in all examined Yorkshire terriers and the similarities to the recessively inherited conditions in SNE-affected Alaskan husky dogs and THMD2-affected humans, as well as the available pedigree

information of the two SNE-affected siblings, a recessive mode of inheritance was postulated. Therefore, homozygosity mapping assuming identity-by-descent (IBD) was used to determine critical genomic regions shared across four SNP array genotyped cases. This revealed five genome regions with a total of ~4.1 Mb located on five different dog chromosomes (Table 1), representing 0.17% of the canine reference sequence. Visual inspection of these regions in the WGS of the affected dog did not reveal any evidence for copy number variants or large structural rearrangements.


**Table 1.** Regions of shared homozygosity detected in four subacute necrotizing encephalomyelopathy (SNE)-affected Yorkshire terriers.

<sup>1</sup> in respect to dog reference genome assembly CanFam3.1.

#### *3.2. Identification of the Causative Variant*

Filtering the variants of a single affected Yorkshire terrier against 729 public control genomes [16], including 60 breed controls, for single-nucleotide variants (SNVs) and short indels present in the five identified IBD-regions resulted in only a single private protein-changing variant (Figure 1a). The indel affecting ~45 bp is located in exon 2 of the *thiamine transporter 2* (*LOC486151*) gene, also known as solute carrier family 19 members 3 (*SLC19A3*) gene (Figure 1b). PCR and subsequent Sanger sequencing confirmed the homozygous presence of this small structural variant in SNE-affected Yorkshire terriers and revealed the detailed features of the indel: a 35 bp insertion replacing 6 bp and thereby disturbing the correct reading frame (Figure 1c). There are 15 currently annotated transcript isoforms for the canine *SLC19A3*, which is in reverse complementary orientation with respect to the canine reference genome. While the canine SLC19A3 protein length is 495 amino acids, the human protein (NP\_001358340.1) has 496 amino acids, from which 408 (82.3%) are identical between dog and human. The Yorkshire terrier variant leads to a frameshift and a premature stop codon (c.205\_210delins35; p.Pro69Ilefs\*45) truncating ~86% of the wild type coding sequence (Figure 1c).

**Figure 1.** Subacute necrotizing encephalopathy (SNE)-associated *SLC19A3* loss-of-function variant in Yorkshire terriers. (**a**) IGV [17] screenshots of the genome region on canine chromosome 25 with the *SLC19A3*:c.205\_210delins35 variant in an affected and a control Yorkshire terrier (NC\_006607.3:40417780-40417930); The indel variant detected in the SNE-affected dog is indicated by a red arrow. (**b**) Schematic representation of the canine *SLC19A3* gene showing the location of both pathogenic variants in exon 2 (XM\_022409850.1): the herein identified indel (red arrow) and the insertion previously described in encephalopathy-affected Alaskan huskies (blue arrow) [12]. Note that the number of 5 -untranslated exons (grey) varies between transcript isoforms, whereas the five protein-coding 3 -exons (black) are more conserved; (**c**) Sanger sequencing electropherograms illustrate sequences of a homozygous SNE-affected Yorkshire terrier, a heterozygous carrier, and a homozygous wild type dog. The red arrows indicate that the 35 bases shown in red are inserted, whereas the 6 bases in blue are deleted in the mutant allele. The predicted consequence of the shift in the reading frame altering the amino acid sequence of the SLC19A3 protein and leading to a premature stop is shown above.

#### *3.3. Targeted Genotyping of the Variant*

Genotyping by fragment size analysis of the 172 available Yorkshire terriers confirmed perfect segregation of the detected *SLC19A3* variant with the observed disease phenotype. Only the ten SNE-affected dogs were homozygous for the variant allele (Table 2). Two obligate carriers and one tested normal littermate were heterozygous carriers of the variant, while 162 controls tested homozygous for the wild type allele (Table 2).

**Table 2.** Segregation of the *SLC19A3*: c.205\_210delins35 genotypes with subacute necrotizing encephalopathy in Yorkshire terriers.


<sup>1</sup> including 60 dogs with WGS data <sup>2</sup> includes 2 obligate carriers and 1 normal littermate of the affected dogs.

#### **4. Discussion**

In this study, the obtained genetic results elucidate the underlying aetiology of the previous clinical and pathological characterization of a Leigh-like subacute necrotizing encephalopathy in the affected Yorkshire terriers, which resembles the human Leigh syndrome. The *SLC19A3* variant found by a combination of SNP genotyping-based homozygosity mapping and whole-genome sequencing, confirmed by Sanger sequencing, segregated perfectly in the investigated cohort of >200 unrelated Yorkshire terriers.

Numerous homozygous as well as compound heterozygous variants have been reported before in different regions of *SLC19A3* in human patients suffering from thiamine metabolism dysfunction syndrome-2 [19]. *SLC19A3* is a member of solute carrier family 19 and encodes thiamine transporter 2. Together with thiamine transporter 1, it is necessary for transport and homeostasis of thiamine that is important in brain development. [20]. *Slc19a3*-knockout mice showed progressive wasting and lethargy leading to a premature death as well as a significant decrease in thiamin uptake, even though there were no obvious histological changes in the brain [21].

The herein-described most likely pathogenic variant (XP\_022265558.1:p.Pro69Ilefs\*45) lies within the second of 12 transmembrane domains of the SLC19A3 protein and, therefore, affects ~86% of the wild type sequence. The *SLC19A3* gene probability of loss-of-function intolerance is pLI = 0.104 [22], which indicates variants in *SLC19A3* leading to a loss of gene function are most likely recessive, where loss of a single copy is often tolerated but the loss of both copies is not. The herein-described variant leads to an insertion of a premature termination in the second out of five coding exons, suggesting that any synthesized mRNA would likely be degraded through nonsense-mediated decay, unlikely to produce a fully functional protein. Heterozygous carriers did not show a visible clinical phenotype, as they can most likely compensate due to the presence of the normal protein, albeit at a decreased amount.

In conclusion, our results provide strong evidence for a breed-specific deleterious variant in *SLC19A3* as the most likely genetic cause of monogenic autosomal recessive Leigh-like subacute necrotizing encephalopathy in Yorkshire terriers, and they enable the development of a genetic test for veterinary diagnostic and breeding decisions. Finally, this presents the second, most likely breed-specific pathogenic variant in the canine *SLC19A3* gene in SNE-affected dogs.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/10/1215/s1, Table S1: Accession numbers of 721 dog and 9 wolf whole-genome genome sequences.

**Author Contributions:** Conceptualization, K.M. and C.D.; methodology, M.D.; data curation, V.J.; formal analysis, M.D. and A.L.; investigation, M.D.; resources, K.M., D.C., M.R., K.J., S.M., T.G., S.R., A.F. and A.L.F.-P.; writing—original draft preparation, A.L., M.D., K.M. and C.D.; writing—review and editing, M.D., A.L., K.M. and C.D.; visualization, A.L. and C.D.; software, V.J.; supervision, C.D.; funding acquisition, C.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to the owners of all dogs who provided samples and shared valuable information. The Next Generation Sequencing Platform and the Interfaculty Bioinformatics Unit of the University of Bern are acknowledged for performing the WGS and providing high-performance computational infrastructure. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *NSDHL* **Frameshift Deletion in a Mixed Breed Dog with Progressive Epidermal Nevi**

**Matthias Christen 1,2, Michaela Austel 3, Frane Banovic 3, Vidhya Jagannathan 1,2 and Tosso Leeb 1,2,\***


Received: 24 September 2020; Accepted: 29 October 2020; Published: 30 October 2020

**Abstract:** Loss-of-function variants in the *NSDHL* gene have been associated with epidermal nevi in humans with congenital hemidysplasia, ichthyosiform nevi, and limb defects (CHILD) syndrome and in companion animals. The *NSDHL* gene codes for the NAD(P)-dependent steroid dehydrogenase-like protein, which is involved in cholesterol biosynthesis. In this study, a female Chihuahua cross with a clinical and histological phenotype consistent with progressive epidermal nevi is presented. All exons of the *NSDHL* candidate gene were amplified by PCR and analyzed by Sanger sequencing. A heterozygous frameshift variant, c.718\_722delGAACA, was identified in the affected dog. In lesional skin, the vast majority of *NSDHL* transcripts lacked the five deleted bases. The variant is predicted to produce a premature stop codon truncating 34% of the encoded protein, p.Glu240Profs\*17. The mutant allele was absent from 22 additionally genotyped Chihuahuas, as well as from 647 control dogs of diverse breeds and eight wolves. The available experimental data together with current knowledge about *NSDHL* variants and their functional impact in humans, dogs, and other species prompted us to classify this variant as pathogenic according to the ACMG guidelines that were previously established for human sequence variants. Therefore, we propose the c.718\_722delGAACA variant as causative variant for the observed skin lesions in this dog.

**Keywords:** *Canis lupus familiaris*; animal model; genodermatosis; dermatology; skin; CHILD syndrome; ILVEN; epidermal nevus; precision medicine

#### **1. Introduction**

Congenital hemidysplasia, ichthyosiform nevi, and limb defects (CHILD) is a well-described syndrome in humans (OMIM #308050) [1–5]. CHILD patients typically have unilateral limb defects together with the characteristic skin lesions termed as CHILD nevi [2]. The severity of the limb defects ranges from aplasia of whole limbs to dystrophy of single finger- or toenails [1,6]. CHILD syndrome is an X-linked semidominant trait. Heterozygous females show the CHILD phenotype, while embryonic lethality is observed in hemizygous males. Due to the random X-chromosome inactivation in heterozygous females, the skin lesions often follow a characteristic pattern, the so-called Blaschko's lines [7].

To date, a wide range of variants in the *NSDHL* gene have been described as suspected causative variants for the pathogenesis of CHILD syndrome in humans [3,8]. In mice, variants in this gene are responsible for the bare patches (*Bpa*) and striated (*Str*) phenotypes [9]. *Bpa* and *Str* mice do not show limb defects, but their skin phenotype closely resembles CHILD nevi [9]. In veterinary medicine, congenital epidermal nevi without limb defects and candidate causative variants in the *NSDHL* gene

have been described in two Labrador Retrievers [10], one purebred Chihuahua [11], and one cat [12] (OMIA 002117-9615, 002117-9685).

The *NSDHL* gene encodes the NAD(P)-dependent steroid dehydrogenase-like protein, which is involved in cholesterol biosynthesis and localized on ER membranes as well as on the surface of lipid droplets [8,13]. A deficiency of this enzyme leads to accumulation of metabolic products from the cholesterol biosynthetic pathway [14,15]. The combination of the toxic effect of those sterol precursors together with decreased cholesterol synthesis and secretion into the stratum corneum is ultimately responsible for the observed skin lesions in CHILD syndrome in humans [16,17].

The currently used therapy for the skin-associated changes in humans and companion animals is aimed at the incorrect cholesterol synthesis. Good results were achieved with topical application of early inhibitors of the cholesterol synthesis pathway in combination with cholesterol-ointments [11,18,19].

In a continuation of our earlier studies [10–12], the present study aimed to characterize the clinical and histopathological phenotype in a female Chihuahua cross with congenital progressive epidermal nevi and to identify the likely causative genetic variant.

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

All animal experiments were performed according to local regulations. The dog in this study was privately owned and was examined with the consent of the owner. The "Cantonal Committee for Animal Experiments" approved the collection of blood samples from control dogs (Canton of Bern; permit 71/19).

#### *2.2. Animal Selection*

A female spayed Chihuahua cross was investigated. The clinical and dermatological examinations were done by two board certified veterinary dermatologists (M.A. and F.B.). An EDTA blood sample was collected for genomic DNA isolation. Routine histopathological examination of skin biopsies was performed before and after euthanasia of the animal.

Additionally, 22 blood samples of purebred Chihuahuas, which had been donated to the Vetsuisse Biobank were used. They represented unrelated controls without skin lesions. The same 22 dogs had been used as controls in our earlier publication [11]. Furthermore, 655 publicly available canine genome sequences [20] were analyzed as control dataset (Table S1); 594 of these had been used in our earlier publication [11]. The skin phenotypes of these 655 control dogs were not consistently documented. As the investigated phenotype is rare and at the same time very striking and obvious, we assumed that these animals were all nonaffected. The available data on their phenotype are summarized in Table S1.

#### *2.3. DNA and RNA Extraction*

Genomic DNA was isolated from EDTA blood samples with the Maxwell RSC Whole Blood DNA Kit using a Maxwell RSC instrument (Promega, Dübendorf, Switzerland). Total RNA was extracted from skin biopsies using the RNeasy Mini Kit (Qiagen, Hilden, Germany). The RNA was cleared of genomic DNA contamination using the QuantiTect Reverse Transcription Kit (Qiagen). The same kit was used to synthetize cDNA, as described by the manufacturer.

#### *2.4. Gene Analysis*

The CanFam3.1 dog reference genome assembly and the NCBI annotation release 105 were used. Numbering within the canine *NSDHL* gene corresponds to the NCBI RefSeq accession numbers XM\_014111859.2 (mRNA) and XP\_013967334.1 (protein).

#### *2.5. PCR and Sanger Sequencing*

Primer pairs for the amplification of all ten exons of the *NSDHL* gene were described previously [8] and are given in Table S2. PCR products for each exon were amplified from genomic DNA using AmpliTaq Gold 360 Master mix (Thermo Fisher Scientific, Reinach, Switzerland). PCR products were visualized using a 5200 Fragment Analyzer instrument (Agilent, Basel, Switzerland), which uses capillary electrophoresis to enable accurate sizing and quantification of nucleic acids. All exons of the *NSDHL* gene were analyzed by direct Sanger sequencing of PCR amplicons. After treatment with exonuclease I and alkaline phosphatase, amplicons were sequenced on an ABI 3730 DNA Analyzer (Thermo Fisher Scientific, Reinach, Switzerland). Sanger sequences were analyzed using the Sequencher 5.1 software (Gene Codes, Ann Arbor, MI, USA). For the RT-PCR on skin cDNA, a forward primer located at the boundary of exons 6 and 7 together with a reverse primer located in exon 10 was used (Table S2). After an initial denaturation of 10 min at 95 ◦C, 35 cycles of 30 s at 95 ◦C, 30 s at 60 ◦C, and 60 s at 72 ◦C were performed, followed by a final extension step of 7 min at 72 ◦C. The RT-PCR products were analyzed on a Fragment Analyzer and sequenced as described above.

#### **3. Results**

#### *3.1. Clinical Examination, Necropsy, and Histopathology*

A female spayed Chihuahua mix was examined at 7 and 27 months of age. The dog presented with a history of band- to plaque-like cutaneous lesions affecting both sides of the body and head since adoption at 3 months of age. The lesions were sharply demarcated and characterized by alopecia, hyperpigmentation, verrucous hyperplasia with pronounced enlargement of follicular ostia, tan-colored scaling, brown to black crusting and, in some areas, peripheral erythema along with malodor (Figure 1a). Band-like lesions often followed the lines of Blaschko (Figure 1a,b). Although both sides of the body were affected, there was a striking lesion severity lateralization to the right side. A somewhat sharp linear demarcation of the lesions was noted along the ventral midline (Figure 1a,b). The paw pads of both front limbs exhibited severe hyperkeratosis, fissuring, and moderate tissue swelling.

Over the course from the first to the second examination, the severity and extent of epidermal nevi progressed to involve diffuse areas on head, both inner ear pinnae, neck, and dorsal and lateral trunk (Figure 1c,d). Pronounced lesional lateralization of the epidermal nevi on the right side of the body was still present (Figure 1c,d). Lesion-associated pruritus and right front limb lameness were initially not observed but became pronounced over time.

Multiple general physical examinations over time yielded findings within normal limits with the exception of the skin. Hematology results and biochemistry parameters showed mild anemia, leukocytosis, hypergammaglobulinemia, and hypoalbuminemia. Schirmer tear test results remained within normal limits over time. Multiple skin scrapes were negative for ectoparasites. Numerous skin cytologies at various points in time revealed mild to partially severe, chronic, secondary bacterial infections with mainly coccoid bacteria. A bacterial tissue culture grew methicillin-resistant *Staphylocoocus pseudintermedius* at one point in time.

Histopathological examination of multiple skin biopsies revealed massive parakeratotic hyperkeratosis mixed with compact lamellar orthokeratotic hyperkeratosis along with irregular, severe acanthosis with broad rete ridges extending deep into the superficial dermis (Figure 1e,f). The hyperkeratosis and acanthosis were also present in follicular infundibular epithelium and were associated with follicular distention supported by an underlying proliferative granular cell layer and stratum spinosum. There were small to large clusters of neutrophils within the stratum corneum and multilaminated, pustular crusts enclosing mixed bacteria covered portions of the epidermis. Hair follicles with sebaceous and apocrine glands were incarcerated in the dermis. Moderate numbers of melanocytes and melanin-filled macrophages were present at the dermal-epidermal junction with mild pigmentary incontinence in the superficial dermis. The dermis contained occasional periadnexal infiltrates of lymphocytes, plasma cells, and neutrophils with some dermal edema. When correlated

with the description of the linear nature of the patient's lesions, the histopathological findings were consistent with an epidermal nevus.

**Figure 1.** Clinical and histopathological phenotype of progressive epidermal nevi in a Chihuahua mix. (**a**,**b**) Clinical phenotype at 7 months of age. Pronounced lateralization of the skin lesions at the initial presentation on the right axillary area (**a**) and right side of ventral/inguinal area (**b**). (**c**,**d**) Clinical phenotype at 27 months of age. Progression of the severity and extent of epidermal nevi to involve diffuse areas in the axillae, ventral abdomen, and inguinal region (**c**) as well as head, both inner ear pinnae, neck, and dorsal and lateral trunk (**d**). (**e**) Histopathology of a skin biopsy taken at 27 months of age from the ventral neck showing severe parakeratotic hyperkeratosis and acanthosis with broad rete ridges extending deep into the superficial dermis. Follicular infundibula are similarly affected by acanthosis and hyperkeratosis and appear distended. (**f**) Details of the histopathology at higher magnification.

The patient's skin underwent histopathological evaluation of different body sites at three different points in time (at 4, 7, and 27 months of age), which all yielded similar results with variations only in regard to the severity of the observed changes and the extent and nature of secondary infections.

Treatments over time included systemic antibiotics (cephalexin, cefpodoxime, clindamycin, and amoxicillin/clavulanic acid), systemic antifungal (fluconazole), antimicrobial shampoos and sprays (chlorhexidine/miconazole), antiseborrheic shampoo and conditioner (gluconolactone-based), antibiotic ointment (mupirocin 2%), topical steroid (betamethasone dipropionate 0.05%), oral glucocorticoid (prednisone), oral vitamin A, oral isotretinoin, oral antihistamine (diphenhydramine), cryotherapy, and interleukin 31 antibody injections, which all yielded minimal to no improvement of clinical signs. Notable exceptions were the interleukin-31 antibody injections, which were accompanied by a reduction in the patient's pruritus, and oral prednisone, which correlated with significantly improved right front paw tissue swelling and lameness.

Due to the steady decline in the patient's overall quality of life, the client elected humane euthanasia at 2 years and 3 months of age. Upon necropsy, no pathologic abnormalities besides the cutaneous lesions and mild, generalized lymphadenopathy were found.

#### *3.2. Genetic Analysis*

As clinical and histopathological findings resembled previously published companion animals with congenital epidermal nevi [10–12], we hypothesized that the phenotype in the affected dog was due to a heterozygous variant in the *NSDHL* gene. Hence, *NSDHL* was investigated as the top functional candidate gene.

Sanger sequencing of all ten exons of the *NSDHL* gene of the examined dog identified a single variant in the coding sequence. The variant is a five base pair deletion in exon 9 of the *NSDHL* gene (Figure 2). The genomic variant designation is NC\_006621.3:120,752,486\_120,752,490delGAACA. It is a frameshift variant, XM\_014111859.2:c.718\_722delGAACA.

**Figure 2.** Details of the *NSDHL*:c.718\_722delGAACA variant. (**a**) Genomic organization of the canine *NSDHL* gene with 10 annotated exons. The position of the variant on the genomic DNA and the mRNA level are indicated. The position of the primers used for the RT-PCR experiment are also indicated. (**b**) RT-PCR amplification products obtained from cDNA from skin of a healthy control and lesional skin from the affected dog. The samples were analyzed on a FragmentAnalyzer capillary gel electrophoresis instrument. The predominant band in the affected dog corresponds to a correctly spliced transcript

lacking 5 nucleotides of coding sequence (r.718\_722del). The minor band that migrated slower than the main product most likely represents heteroduplex molecules consisting of one strand containing the 5-nucleotide deletion annealed to a wildtype strand. (**c**) Relative quantification of these two bands. Please note that the ratio of these end-point RT-PCR amplicons should only be seen as a semiquantitative proxy for the true ratio of transcripts. (**d**) Sanger electropherograms derived from genomic PCR products of a control dog and the affected dog. The five bases that are deleted in the affected animal are indicated with dotted lines. The arrow indicates the beginning of overlapping signals in the electropherogram of the affected dog. The intensities of the overlapping signals correspond well to the expected 50:50 allelic ratio for a heterozygous animal. (**e**) Analysis of mRNA. RT-PCR amplicons from the skin of a normal control dog and lesional skin of the affected dog were sequenced. The vast majority of *NSDHL* transcripts of the affected dog lacked 5 nucleotides corresponding to positions 718-722 in the coding sequence. Note the difference in relative signal intensity of the two alleles in the cDNA amplicon compared to the genomic sequence in the affected dog.

RT-PCR on RNA derived from a lesional skin biopsy of the affected dog yielded an amplicon of the expected size. Sanger sequencing confirmed that the transcripts were normally spliced. We observed pronounced allelic imbalance in lesional skin of the affected dog. The *NSDHL* transcripts were almost exclusively expressed from the mutant allele and lacked the five deleted nucleotides (XM\_014111859.2:r.718\_722del; Figure 2b,c,e). The frameshift deletion truncates 122 codons (34%) from the open reading frame of the wild-type *NSDHL* transcript, XP\_013967334.1:p.(Glu240Profs\*17).

We genotyped 22 unaffected Chihuahuas for the frameshift variant. None of these control dogs carried the mutant allele. The mutant allele was also absent from whole-genome sequence data of 647 control dogs of diverse breeds and eight wolves (Table S1).

#### **4. Discussion**

In this study, we identified a heterozygous *NSDHL*: c.718\_722delGAACA frameshift variant in a mixed breed dog with a severe form of progressive epidermal nevi. The *NSDHL* gene encodes the enzyme NAD(P) dependent steroid dehydrogenase-like, which is involved in cholesterol biosynthesis and has been associated with human CHILD syndrome [8], as well as canine and feline congenital epidermal nevi [10–12].

In this single case investigation, we attempted to apply the American College of Medical Genetics and Genomics (ACMG) approved guidelines for the interpretation of sequence variants in human genetics to a dog with the *NSDHL*:c.718\_722delGAACA variant [21]. Three arguments support the pathogenicity of the *NSDHL*:c.718\_722delGAACA variant.

Computational and predictive data provide a very strong argument for pathogenicity (PVS1) as the frameshift variant is predicted to result in a loss-of-function (LOF) of the *NSDHL* gene [21]. *NSDHL* LOF variants in human and animal patients represent an established mechanism of disease [8–10].

Population data indicate that the mutant allele was only present in a single affected dog, while it was absent from more than 600 control dogs. We rate this as a moderate evidence for pathogenicity (PM2) [21].

Other data finally provide additional supporting evidence for pathogenicity (PP4) [21]. The clinical and histopathological skin phenotype is highly specific and *NSDHL* is the only known candidate gene for such a phenotype. Furthermore, the distribution of the skin lesions along Blaschko's lines and the female sex of the patient strongly point to an X-chromosomal causative variant.

According to the established ACMG criteria, the combinations of one very strong, one moderate, and one supporting type of evidence is sufficient to classify the *NSDHL*: c.718\_722delGAACA variant as pathogenic.

To the best of our knowledge, the found variant is the third disease-causing *NSDHL* variant identified in dogs. Previous studies reported a 14 kb deletion spanning the last three exons of *NSDHL* gene [10] and a single missense variant p.Gly234Arg [11] in dogs with congenital epidermal nevi.

#### **5. Conclusions**

We identified a frameshift variant, *NSDHL*: c.718\_722delGAACA as the probable causative variant for progressive epidermal nevi in a single mixed breed dog. Established guidelines for the interpretation of human sequence variants justify the classification of this variant as pathogenic.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/11/1297/s1, Table S1: Accession numbers of 647 dog and eight wolf sequences; Table S2: Primer sequences for the *NSDHL* gene.

**Author Contributions:** Conceptualization, T.L.; investigation, M.C., M.A., F.B., V.J., and T.L.; data curation, V.J.; writing—original draft preparation, M.C., M.A., and T.L.; writing—review and editing, M.C., M.A., F.B., V.J., and T.L.; supervision, T.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to the dog owner who donated samples and participated in the study. We thank Nathalie Besuchet Schmutz for expert technical assistance, Jan Henkel for help with the RNA experiments, the Next Generation Sequencing Platform of the University of Bern for performing the high-throughput sequencing experiments, and the Interfaculty Bioinformatics Unit of the University of Bern for providing high-performance computing infrastructure.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **Mapping of Diabetes Susceptibility Loci in a Domestic Cat Breed with an Unusually High Incidence of Diabetes Mellitus**

**Lois Balmer 1,2,**†**, Caroline Ann O'Leary 3,**†**, Marilyn Menotti-Raymond 4, Victor David 5, Stephen O'Brien 6,7, Belinda Penglis 8, Sher Hendrickson 9, Mia Reeves-Johnson 3, Susan Gottlieb 3, Linda Fleeman 10, Dianne Vankan 3, Jacquie Rand 3,11,**‡ **and Grant Morahan 1,\*,**‡


Received: 29 September 2020; Accepted: 11 November 2020; Published: 19 November 2020

**Abstract:** Genetic variants that are associated with susceptibility to type 2 diabetes (T2D) are important for identification of individuals at risk and can provide insights into the molecular basis of disease. Analysis of T2D in domestic animals provides both the opportunity to improve veterinary management and breeding programs as well as to identify novel T2D risk genes. Australian-bred Burmese (ABB) cats have a 4-fold increased incidence of type 2 diabetes (T2D) compared to Burmese cats bred in the United States. This is likely attributable to a genetic founder effect. We investigated this by performing a genome-wide association scan on ABB cats. Four SNPs were associated with the ABB T2D phenotype with *p* values <0.005. All exons and splice junctions of candidate genes near significant single-nucleotide polymorphisms (SNPs) were sequenced, including the genes *DGKG, IFG2BP2, SLC8A1, E2F6, ETV5, TRA2B* and *LIPH*. Six candidate polymorphisms were followed up in a larger cohort of ABB cats with or without T2D and also in Burmese cats bred in America, which exhibit low T2D incidence. The original SNPs were confirmed in this cohort as associated with the T2D phenotype, although no novel coding SNPs in any of the seven candidate genes showed association with T2D. The identification of genetic markers associated with T2D susceptibility in ABB cats will enable preventative health strategies and guide breeding programs to reduce the prevalence of T2D in these cats.

**Keywords:** diabetes mellitus; Burmese cats; susceptibility; single-nucleotide polymorphism; genetic markers; LIPH

#### **1. Introduction**

Type 2 diabetes mellitus (T2D) is a polygenic disease with complex inheritance [1,2]. Studies of T2D in humans have identified over 80 genetic variants associated with disease susceptibility, but most of these contribute a small proportion of overall risk of disease, and the detailed molecular basis of how these contribute to disease is not well understood [3]. Animal models can provide additional evidence of the role of genes in T2D. For example, a polymorphism in the melanocortin 4 receptor gene (MC4R) was associated with T2D in overweight domestic cats, similar to humans [4]. Animal models can also provide a valuable opportunity to identify novel T2D genes. We demonstrated previously that Australian-bred Burmese (ABB) cats are at increased risk of developing T2D, and exhibit features typical of human T2D [5,6], including inadequate insulin secretion (dysfunctional β cells); impaired insulin action (insulin resistance) [7]; late age of onset; risk factors such as obesity and physical inactivity; islet vacuolation and amyloid deposition. Clinical features of diabetes in ABB cats are similar to atypical T2D in humans, which is most common in African-Americans, and includes presentation with very high blood glucose concentrations [8,9]. Diabetic ABB cats respond to insulin in a very similar way to human patients. Affected humans present with a history of polyuria, polydipsia and weight loss, but otherwise have a T2D phenotype and profile, and frequently have a family history of T2D. These patients respond to initial insulin treatment within days to weeks with long-lasting, insulin-free, near-normoglycaemic remission [9]. Similarly, cats with diabetes have high remission rates (67 to >80%) resulting in normoglycaemia without the need for insulin, if tight glycaemic control is instituted soon after diagnosis [10,11].

In a study of a cohort of 12,576 cats, the prevalence of diabetes was significantly higher in ABB cats than in domestic short or longhaired cats (*p* < 0.001), with the incidence rising to 10% in Burmese cats over eight years old [5]. Other Australian domestic cat breeds exhibit a diabetes incidence of ~0.25% to 1% [12]. American, European-and ABB cats have had distinct breeding histories since the 1970s [13]. Burmese cats bred in New Zealand and the UK, from where the ABB cats are believed to have originated, also demonstrate increased risk of diabetes [14,15]. However, Burmese cats in other countries such as the US do not show such increased T2D incidence [16,17]. US-bred Burmese cats did have different incidence rates for other genetic diseases [18].

Type 2 diabetes in ABB cats likely arises from a genetic founder effect. Several population bottlenecks have been experienced in cat evolution, including domestication and subsequently in the creation of breeds (mostly within the last 50–150 years) [19]. Of most relevance is the recent bottleneck associated with the establishment of the breed in Australia in 1957 from a very small number of individuals, possibly as few as five [13,20]. Founder effects can result in otherwise rare variants being established in a population and have had great utility in discovery of disease genes in domestic animal species. Thus, the genetic architecture of the ABB cat population provides an opportunity to identify novel genetic variants mediating T2D susceptibility. Identification of such genes could increase our understanding of molecular events leading to T2D and would be immediately applicable in the veterinary setting by identifying cats at risk of diabetes and improving breeding programs. Preventative management could also be implemented for at-risk individuals via dietary changes, weight loss and therapies to control glycaemia.

We hypothesize that genetic variant(s) predisposing to T2D demonstrate a higher frequency in ABB cats compared to American Burmese cats, thus providing a unique opportunity for discovery of diabetes-associated genes. We performed an exploratory case-control genome-wide association scan (GWAS), followed by validation of candidate SNPs in an expanded cohort. A similar approach identified candidate genes in other diseases [21,22].

#### **2. Materials and Methods**

This study was conducted in accordance with the Declaration of Helsinki and the protocol was approved by the Ethics Committee of The University of Queensland (AEC Approval Number: SVS/040/10/NCI/ABBOTT, date: 150910 to 150912; AEC Approval Number: SVS/200/12/NESTLE PURINA/ABBOT ANIMAL HEALTH, date: 180712 to 180715; AEC Approval Number: SVS/256/13/MEDINCELL/NCI, date: 311013 to 311016).

#### *2.1. Characterisation of Australian-Bred Burmese Cats*

Two ABB cat cohorts were studied. The first ("exploratory") cohort contained 10 cats diagnosed with T2D and 10 control non-diabetics with normal glucose tolerance. Of the diabetic cats, four were male. Samples were obtained for five of these from cat-specific practices in Brisbane, Australia, one individual from a general practice, and 4 individuals from a veterinary commercial diagnostic laboratory (*IDEXX*, Brisbane, Australia) that had received blood samples from veterinarians as part of the diagnostic work-up for cats that had clinical signs suggestive of diabetes.

Cats were classed as T2D if they were first diagnosed at 8 years of age or older based on hyperglycaemia (blood glucose concentrations over 15 mmol/L) and clinical signs of polydipsia and polyuria. At the time of sampling, 4 cats were in remission and maintained euglycaemia without insulin. Selecting cats above 8 years of age reduced the risk of including cats with type 1 diabetes, which develops at an earlier age [23] and is relatively rare [24]. Inclusion of cats with other types of diabetes causing insulin resistance, most commonly acromegaly, was unlikely as cats were treated with typical doses of insulin and were not reported to have consistent clinical and/or pathology findings of acromegaly [25]. Further, 5/7 cats had feline pancreatic lipase immunoreactivity (fPLI) below 8 μg/L [26] and two had values below 12 μg/L and therefore not consistent with a diagnosis of active pancreatitis [27]; 7/8 tested had normal total thyroxine (TT4) and one had hyperthyroidism (TT4 107), one had liver disease, 3 had mildly increased liver enzymes, and 4 had renal dysfunction.

The 10 control ABB cats (five males) were all older than 8 years and exhibited fasting and 2 h blood glucose concentrations in a glucose tolerance test within reference ranges for age-matched cats. They were clinically healthy based on history, physical examination, routine haematological, routine biochemical, feline pancreatic lipase, and TT4 testing. For the glucose tolerance test, cats were fasted overnight and 0.5 g/kg 50% glucose injection BP was injected via a cephalic catheter placed 3 h before and samples collected from the pinna or paw at 0, 2, 3 and 4 h. All cats hada2h glucose concentration <10 mmol/L (upper 95% for normal reference range for glucose concentration 2 h following 0.5 g/kg glucose) and all blood glucose concentrations had returned to <6.5 mmol/L (upper limit 95% reference range for fasting concentration) by 2 to 2.5 h. All control cats had fPLI measurements <7.6 μg/L, 8 of 8 tested had normal TT4, and all had normal haematology and biochemistry.

The second cohort consisted of 84 ABB; 37 with T2D and 47 controls. Of the T2D cases, 24 were male and 13 female. Eight were from a diabetes-specific practice in Melbourne and 4 from a cat-specific practice in Brisbane. Of these latter 12 cases, 5 were in remission and the remaining 7 were stable on insulin doses <1.5 U/kg. Evidence of renal dysfunction was present in 5/7, hyperthyroidism in 1/5, and liver disease in 1/7. Five of these 12 cats had fPLI values available, with 4 being below 8 μg/L. The remaining 25 samples were obtained Australia-wide by the veterinary diagnostic laboratory *IDEXX*. All cases were either diagnosed with blood glucose over 15 mmol/L; 8 were newly diagnosed, 2 were euthanased at diagnosis, 11 were on insulin <1.5 U/kg and 4 had unknown insulin treatment histories. Veterinarians caring for cats were contacted and the history, clinical details and therapeutic history of the cases were obtained to ensure these cats did not have other forms of diabetes, such as acromegaly. Further, 14/20 cats had fPLI <12 μg/L, 24/24 had normal total TT4. On haematology and biochemistry (and urinalysis where available), biochemical evidence of renal dysfunction was present in 6 of 26, liver disease in 7, renal and liver disease in 5, and there were no data for 1 animal.

ABB control cats in cohort two (*n* = 47) were eight years or older, with a screening blood glucose concentration (measured after entry to the clinic and any time in relation to eating) of <8.1 mmol/L measured from fluoride-oxalate samples on an automated serum chemistry analyser (*n* = 42 from *IDDEX*), or immediately after sample collection using a handheld glucose meter calibrated for cats (Abbott AlphaTRAK) (*n* = 5 from research projects) [28]. Of these, 14 were male, 33 female. Screening blood glucose concentration <8.1 mmol/L is lower than the upper limit of the screening blood glucose reference range (9.7 mmol/L) and in 45/47 cats, screening blood glucose was <6.5 mmol/L, a stringent cut-off point for aged, client-owned cats that have not been fasted and are potentially subject to stress (which may cause hyperglycaemia) [29]. Burmese cats without signs of diabetes have been shown to have fasting blood glucose concentrations on average 1.9 mmol/L higher (*p* < 0.05) than age-matched domestic cats [5]. Of these 47 cats, one had no haematology, biochemistry or urinalysis data available, evidence of renal dysfunction was present in 15, liver disease in 2, liver and renal disease in 6, and 23 had normal haematology and biochemistry (and urinalysis when available) findings.

#### *2.2. US-Bred Burmese Cats*

A third cohort comprised 84 DNA samples from American-bred Burmese cats; 37 female, 21 male and 26 of unknown sex. These cats were not tested for diabetes, but are representative of a population with a normal low risk of type 2 diabetes (~0.25% to 1%) [17].

#### *2.3. Genome-Wide Association Mapping and Analyses for Association with Diabetes: Genomic*

Genomic DNA was extracted using a QIAamp DNA Blood Mini Kit (QIAGEN, Hilde, Germany). Genotyping in the first cohort of cats was performed using the Illumina iSelect custom feline genome chip containing 58,444 SNPs. All SNPs were analysed with Illumina Genome Studio software and subjected to stringent quality control procedures. SNPs with GenTrain Scores <0.60, call rate <0.95, or failed genotype concordance in replicates were removed from analysis. SNPs were excluded if they exhibited a minor allele frequency (MAF) in cases <0.01 due to sample size, or case/control bias in missingness. In total, 38,487 SNPs remained after these cut offs. Association testing was performed in PLINK [30] using Fisher's Exact test, with *p* values corrected for multiple tests by Bonferroni correction. Correction for genetic admixture (caused by potential misrecorded matings) was assessed and corrected using Eigensoft [30]. The initial screening GWAS, with 10 cases and 10 controls, was regarded primarily as an exploratory study, with any suggestive findings to be followed in a subsequent cohort.

#### *2.4. PCR, Gel, Digest and High-Resolution Melt Analysis*

Targeted genotyping assays were developed to analyse the novel variants that we identified. Depending on the SNP, assays used were allele-specific PCR, single-strand conformation polymorphism PCR, restriction digestion or high-resolution melt analyses. Primer sequences are in Supplementary Table S1.

In general, 20 ng of genomic DNA from each sample was amplified by PCR in a 20 μL total PCR reaction mixture including 0.2 U Platinum Taq DNA polymerase (Life Technologies, Perth, Australia), 50 mM MgCl2 (Life Technologies, Victoria, Australia), 10 mM nucleoside triphosphates (dNTPs) (SIGMA-Aldrich, Castle Hill, Australia) and 50 ng primers (SIGMA-Aldrich, Castle Hill, Australia). The cycling parameters were: one cycle at 96 ◦C for 5 min, then 34 cycles of 96 ◦C for 20 s, 58–66 ◦C for 20 s, 72 ◦C for 30 s and a final extension at 72 ◦C for 3 min. The PCR products were separated on a 3% agarose gel (BIO-RAD Laboratories, Hercules, CA, USA).

Genotyping was performed where possible by differential restriction enzyme digestion. After digesting 15 μL of PCR product for 3 h with the appropriate enzyme (Table 2), and manufacturers' instructions, the products were electrophoresed on a 3% agarose gel.

If PCR products (such as SNP6, chrC2:84361252+84361412), could not be distinguished by gel separation, they were tested by High Resolution Melt analysis. Briefly, they were amplified from 20 ng of genomic DNA in a 20 μL total PCR reaction mixture containing ResoLight Dye (Roche, Indianapolis, IN, USA), 0.2 U Platinum Taq DNA polymerase (Life Technologies, Victoria, Australia), 50 mM MgCl2 (Life Technologies, Victoria, Australia), 10 mM dNTPs (SIGMA-Aldrich, Castle Hill, Australia) and

50 ng primer (SIGMA-Aldrich, Castle Hill, Australia). Following PCR, the samples were heated to 95 ◦C, rapidly cooled to 45 ◦C, reheated to 70–95 ◦C with a ramp rate of 0.2 ◦C/s, and analysis of melt curves was conducted using the Light Cycler 480 Software Release 1.5.0 SP3 (Roche, Sydney, Australia) [31].

Allele-specific PCRs were used to detect the G/A polymorphism in the *E2F6* gene. The allele-specific reaction mixture contained 20 μg of genomic DNA in a 20 μL PCR reaction mixture including 0.2 U Platinum Taq DNA polymerase (Life Technologies, Victoria, Australia), 40 mM MgCl2 (Life Technologies, Victoria, Australia), 10 mM dNTPs (SIGMA-Aldrich, Castle Hill, Australia) and 50 ng allele specific primers (SIGMA-Aldrich, Castle Hill, Australia). Primer sequences are indicated in Supplementary Table S1. Forward primers were designed with a T-A nucleotide mismatch to create primer 3 instability and increased allele specificity. The G and A 3 -terminal nucleotides of the allele-specific forward primers annealed to the G and A nucleotides of the *E2F6* gene. The reverse primer was designed to anneal to a highly conserved region of the *E2F6* gene in Australian-bred Burmese. The thermal profile consisted of one cycle at 96 ◦C for 5 min, then 34 cycles of 96 ◦C for 20 s, 58–66 ◦C for 20 s, 72 ◦C for 30 s, and a final extension at 72 ◦C for 3 min. Reactions were performed to determine AA/AG/GG status of all ABB DNA samples. The PCR products were separated on a 3% agarose gel (BIO-RAD Laboratories, Hercules, CA, USA).

#### *2.5. Sequencing of Candidate Genes*

Genes near the SNPs identified in the GWAS as loci of interest (Table 1) were inspected using the November 2017 ICGSC Felis\_catus 9.0/felCat9 and candidate genes were chosen for sequencing On chromosome A3, SLC8A1 located within 20 Kb of SNP1 and E2F6 was within 20 Kb of SNP2. On chromosome C2, SNP4 was located adjacent to 3 candidate genes; ETV5 within 20 Kb, DGKG within 200 Kb, and TRA2B within 150 Kb. On chromosome C2, SNPs 5–6 were located near IGF2BP2 within 200 Kb, and LIPH within 100 Kb.

**Table 1.** Candidate SNPs identified from the exploratory GWAS. A1 indicates the minor allele (i.e., least frequently observed based on the whole sample); F\_A and F\_U indicate the frequency of the minor allele in the affected and unaffected cats, respectively. The *p* value is shown calculated by Fisher's Exact test and Bonferroni corrected by PLINK. No odds ratio could be calculated for SNP2 because the associated allele was not found in unaffected cats. SNPs were named in the Illumina array only as custom numbers; for convenience, they are named here as SNPs 1 to 6.


Exons and the 5 and 3 untranslated regions (UTR) of all candidate genes were sequenced (primer sequences shown in Supplementary Table S2). To amplify the DGKG gene, M13 tails were added to all primers, and a touchdown PCR protocol was followed [32]. Sequencing for the DGKG gene was performed with Version 1.1 big dye terminator sequencing using M13 forward and reverse primers [33]; the remaining genes sequenced with the PCR primers. PCR products were cleaned using the PCR Clean up Kit (Qiagen), 20 ng/μL of this product was then mixed with 20 ng/μL of the forward or reverse primer (Sigma-Aldrich, Castle Hill, Australia), Big Dye Terminator (V1.3) and Buffer (×5) (Ambion). The cycling parameters were: one cycle at 96 ◦C for 3 min, then 25 cycles of 96 ◦C for 10 s, 50 ◦C for 2 s, 60 ◦C for 4 min. The sequencing reactions were cleaned and analysed on the ABI3730. Sequences were analysed using Sequence Scanner 2 software (Applied Biosystems, Victoria, Australia). Sequences were aligned using Clustal Omega 2019 (EMBL-EBI, Cambridge, UK) to identify SNPs. Variant protein

sequences were submitted to Protein Variation Effect Analyser (PROVEAN V1.0) (Rockville, MD, USA) protein prediction software [34]. No further analysis was conducted on SNPs that did not result in protein coding changes.

#### **3. Results**

#### *3.1. Search for Genetic Markers Associated with Type 2 Diabetes*

Our strategy was to conduct an exploratory GWAS with 10 diabetic and 10 control ABB cats to characterise sequence polymorphisms between diabetic and non-diabetic Burmese which might be examined in a larger population sample, and to validate any candidate polymorphisms in a larger cohort. Genotyping was performed using custom Illumina arrays [35]. As expected for a sample this small, no SNP provided a *p* value below the threshold for significance after Bonferroni correction (i.e., *<sup>p</sup>* <sup>8</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>7). However, six SNPs were identified with promising *<sup>p</sup>* values <sup>&</sup>lt;0.0005 (Table 1). These SNPs tagged regions on cat chromosomes A3, C1 and C2. The SNPs on the custom Illumina chip had custom numerical annotations. For convenience, these six candidate SNPs are designated here as SNPs 1 through 6.

#### *3.2. Validation of the Candidate Type 2 Diabetes Associated SNPs*

The six candidate SNPs were genotyped in a second cohort of ABB cats supplemented with samples from Burmese cats from the low-prevalence US population. Only one (SNP6) was not significant in the validation cohort. SNPs 1–4 demonstrated significant association after Bonferroni correction while a further SNP (SNP5) had a *p* value of 0.01 (Table 2).

**Table 2.** Testing candidate SNPs in a validation cohort. DNA samples from a further 37 diabetic ABB cats, 47 normoglycaemic Australian-bred Burmese cats, and 84 Burmese cats from the low-prevalence US population were genotyped at the SNPs implicated by the GWAS. Results were analysed using Plink. A1 indicates the minor allele (i.e., least frequently observed based on the whole sample); F\_A and F\_U indicate the frequency of the minor allele in the affected and unaffected cats, respectively. The *p* value is shown calculated using Fisher's Exact test and Bonferroni correction. No odds ratio could be calculated for SNP2 because the associated allele was not found in the non-diabetic cats.


Haplotype analyses demonstrated that SNPs on chromosome C2 tagged intervals that were strongly associated with type 2 diabetes (Table 3). Even the SNP that was not significant as a singleton (SNP6), demonstrated a significant association when combined with SNP5, tagging a chromosome interval that showed association. As these SNPs span a genomic interval of well over 500 kb, these results suggest that they may tag a founder haplotype that increases the risk of type 2 diabetes, or that the individual variants near these SNPs interact to increase susceptibility.

**Table 3. Haplotype analyses**. Data from the validation cohort were analysed using the haplotype association test implemented in PLINK. Results for the most associated haplotype are shown.


#### *3.3. Selection of Candidate Type 2 Diabetes Genes*

Genes near the associated SNPs were inspected using the November 2017 ICGSC Felis\_catus 9.0/felCat9. The two associated SNPs on chromosome A3, SNP1 and SNP2, were located within introns of *SLC8A1* and *E2F6*, respectively. There are no other known genes within 50 Kb 5 or 3 of these SNPs. Neither of these genes has previously been implicated in type 2 diabetes susceptibility. The most highly associated SNP, SNP3, resided in an area with no known coding genes within at least 300 Kb.

SNP4 is located in an intron of *ETV5*. The human orthologue of *ETV5* was reported to be associated with type 2 diabetes and obesity in several human populations [36,37]. SNPs 5 and 6 are within the same intron of the *LIPH* gene, which encodes a membrane-bound triglyceride lipase. The genomic region on chromosome C2 located between SNPs 4 and 6, contains only four other known genes based on sequence homology or conserved synteny with other species. These are: *TRA2B*, *CHCHD4*, *IGF2BP2* and *SENP2*. Of these, *TRA2B* (also known as SFRS10) has been linked in humans and mice to obesity and so to related metabolic syndromes including type 2 diabetes [38]. Studies have indicated the involvement of *IGF2BP2* in type 2 diabetes in humans [39].

#### *3.4. Search for Polymorphisms within Candidate Genes*

To identify genetic variants in the genes identified above, all coding regions of these genes were sequenced, including the 3 and 5 UTRs of each. Information about the genomic regions sequenced, any discovered SNPs, the nucleotide position of each SNP, any predicted amino acid change and the effect of the mutation as estimated by the PROVEAN algorithm is shown in Supplementary Table S3. Association results of selected SNPs with the diabetes trait are shown in Supplementary Table S4. Eleven SNPs were discovered in *SLC8A1*, of which two were found to contain a potentially deleterious amino acid substitution. Twelve SNPs were identified in *E2F6*, with one predicted to be deleterious. However, none of these SNPs were significantly associated with the disease phenotype, suggesting that they were commonly occurring in the Burmese cat breed but not the cat(s) used to establish the genome sequence of this species. No likely causal mutations, including non-synonymous substitutions, polymorphisms in splice junctions, introduction of premature stop codons or insertions/deletions resulting in frame-shifts were found in either *IGF2BP2* or *DGKG*. Fourteen SNPs were identified in *ETV5* with two predicted to be deleterious, but these also were not significantly associated with the disease phenotype. Two SNPs and one deletion were identified in the *TRAF2B* gene; all were found to reside within introns.

Finally, eight novel SNPs were identified in *LIPH*, none of which were in exons. One of these SNPs was significantly (*p* = 1.6 <sup>×</sup> 10−5) associated with diabetes, conferring a relative risk of 2.6. No feline transcripts of this gene have been sequenced, but by alignment of the cat genomic and human cDNA sequences, the significant SNP is located in the presumed 3'UTR region. This is reminiscent of the 3'UTR SNP in the *IL12B* gene which affects gene expression and is associated with susceptibility to type 1 diabetes and other diseases in humans [40,41]. This novel SNP was confirmed and was significantly associated with type 2 diabetes in the total cohort of ABB cats.

#### **4. Discussion**

Our exploratory GWAS found an association of six SNPs with type 2 diabetes. Seven candidate genes were investigated for the presence of variants that could contribute to type 2 diabetes. Of these candidates examined, only *LIPH* exhibited a novel disease-associated variant. Two genes near SNP 4, *ETV5* and *DGKD*, did not have disease-associated sequence variants. The region on Chr C2 near SNPs 5 and 6, contained four known genes (*TRA2B*, *CHCHD4*, *IGF2BP2* and *SENP2*) but none of these genes contained disease-associated protein-coding changes, despite *TRA2B* and *IGF2BP2* having association with type 2 diabetes in other species [39]. Thus, these gene variants compared to the reference cat genome sequence are found in the Burmese cat lineage but are not associated with T2D susceptibility.

SNP3, the most highly associated SNP, is in a genomic region of low gene density which has no known coding genes within at least 300 Kb. Many SNPs associated with disease risk occur in such "gene deserts". They may reside in enhancer elements [42] and regulate expression of distant genes, as shown for type 1 diabetes risk genes [43].

This study demonstrates the potential of ABB cats to identify novel genes involved in type 2 diabetes susceptibility. The ABB cat population originated from a small number of founders, with subsequent inbreeding leading to extensive linkage disequilibrium (LD) [44]. In US-bred Burmese cats, a more outbred population, average haplotype blocks have been reported to be 0.5 Mb [44], along with an average inbreeding coefficient of 0.22 and observed heterozygosity of 0.4 [45]. Extended LD observed in many domesticated animal breeds can be created by population bottlenecks, small effective population sizes and use of popular sires, and has proved highly beneficial in genetic mapping [45]. Thus, in this study, the ABB cat, with its genetic architecture and high prevalence of type 2 diabetes, allowed identification of genetic elements of large effect size using smaller sample sizes than required for outbred populations.

Up to 7% of Australian cats seen by veterinarians are Burmese [12], and up to 10% of ABB cats over 8 years old develop diabetes [5], so ABB cats comprise a significant proportion of the population of domesticated cats requiring veterinary care. Diabetes has major quality of life implications for cats and also major cost implications for owners [46]. Our results have clinical relevance for Burmese cats at increased risk of diabetes in Australia, and in other countries like the UK [14] and New Zealand [15], and possibly for other domestic cat populations.

Our results could reduce the burden of care on ABB cat owners by identifying cats carrying risk alleles. Identifying at-risk cats allows early intervention to identify pre-diabetic cats and prevent development of overt diabetes. Early identification of pre-diabetic cats is important, once cats are pre-diabetes (fasting blood glucose >7.5 mmol/L and 3 h glucose >14 mmol/L), they have an 88% probability of being diabetic within 9 months [47]. Potential interventions include weight control, a low carbohydrate diet to reduce postprandial increases in glycaemia, and insulin-sensitising drugs. Currently such drugs are not commonly used because of the difficulty in identifying at-risk pre-diabetic cats. Further, general dietary and body weight recommendations for cats are not always adopted by owners but may have increased uptake if owners knew their cat was genetically at risk. Identification of pre-diabetic and subclinical diabetic cats could result in these cats being managed without the need for insulin. Further, improved monitoring would expedite early diagnosis of diabetes and quick implementation of tight glycaemic control. This has benefits for diabetic cats and their owners because remission rates in excess of 80% can be achieved compared to 30–40% remission rates if tight glycaemic control is delayed [48].

This study, by identifying risk loci in the ABB cat population, particularly allows cat breeders to make informed choices to avoid type 2 diabetes risk in their breeding pedigrees. Cats are bred between 1 and 6 years, well before clinical signs of diabetes are evident, so identification of carriers of risk loci would help reduce type 2 diabetes prevalence. Further investigations may also show whether these loci may also contribute to diabetes in other cat breeds.

Type 2 diabetes in ABB cats is clearly not monogenic; while we have mapped a few loci, their susceptibility alleles are not present in all diabetic cats genotyped, and doubtless others with smaller effects are also involved. Despite this, our identification of risk alleles is a step toward a greater understanding of the metabolic basis of diabetes in cats. Identification of risk alleles in cats may also have relevance for the human disease [49], and this cat population could provide an animal model for studying the role of these genes in the disease process.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/11/1369/s1. Supplementary Table S1: Nucleotide sequence of primers used for genotyping candidate SNPs, Supplementary Table S2: Primers used to amplify and sequence fragments of candidate genes, Supplementary Table S3: Search for polymorphisms within candidate genes, and Supplementary Table S4: Testing Novel SNPs for Asociation with Type 2 Diabetes.

**Author Contributions:** L.B. and G.M. sequenced the genes, conducted the analysis and GWAS analysis and helped write the manuscript. L.B., G.M., C.A.O. and J.R. wrote the manuscript and devised the project, M.M.-R., V.D., S.O., B.P., S.H., M.R.-J., S.G., L.F. and D.V. helped with the design of the project, collection of data, DNA extraction and results coordination. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported [in part] by the Intramural Research Program of the NIH, Frederick National Laboratory, Center for Cancer Research, National Cancer Institute. L.B. and G.M. were supported by Diabetes Research Western Australia (DRWA) and NH&MRC program grant.

**Acknowledgments:** The authors thank C Bloom and T Evans, both of whom provided a sample of feline blood.

**Conflicts of Interest:** C.A. O'Leary, none; L. Balmer—none; M. Menotti-Raymond—none; V.A. David—none, S.J. O'Brien—none; B. Penglis—none; S. Hendrickson—none; M. Reeves-Johnson—none; S. Gottlieb—none; L. Fleeman—none; D. Vankan—none; G. Morahan—none; J. Rand—none.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **The Genetic Basis of Obesity and Related Metabolic Diseases in Humans and Companion Animals**

#### **Natalie Wallis and Eleanor Ra**ff**an \***

Anatomy Building, Department of Physiology, Development and Neuroscience, University of Cambridge, Downing Street, Cambridge CB2 3DY, UK; njw64@cam.ac.uk

**\*** Correspondence: er311@cam.ac.uk

Received: 12 October 2020; Accepted: 16 November 2020; Published: 20 November 2020

**Abstract:** Obesity is one of the most prevalent health conditions in humans and companion animals globally. It is associated with premature mortality, metabolic dysfunction, and multiple health conditions across species. Obesity is, therefore, of importance in the fields of medicine and veterinary medicine. The regulation of adiposity is a homeostatic process vulnerable to disruption by a multitude of genetic and environmental factors. It is well established that the heritability of obesity is high in humans and laboratory animals, with ample evidence that the same is true in companion animals. In this review, we provide an overview of how genes link to obesity in humans, drawing on a wealth of information from laboratory animal models, and summarise the mechanisms by which obesity causes related disease. Throughout, we focus on how large-scale human studies and niche investigations of rare mutations in severely affected patients have improved our understanding of obesity biology and can inform our ability to interpret results of animal studies. For dogs, cats, and horses, we compare the similarities in obesity pathophysiology to humans and review the genetic studies that have been previously reported in those species. Finally, we discuss how veterinary genetics may learn from humans about studying precise, nuanced phenotypes and implementing large-scale studies, but also how veterinary studies may be able to look past clinical findings to mechanistic ones and demonstrate translational benefits to human research.

**Keywords:** obesity; genetics; companion animals; metabolic disease; comparative genomics; dogs; cats; horses

#### **1. Introduction**

Obesity presents a major health problem in humans and companion animals alike. An estimated 39% of people were overweight or obese in 2016, a value nearly triple that recorded in 1975 and equating to over 2 billion adults [1]. Mirroring that trend are increases in pet obesity, with as many as 63% of cats [2] and 59% of dogs [3] reported to be overweight or obese. Obesity was declared an epidemic by the World Health Organisation (WHO) in 1997 [4] and similarly identified as a major threat to pet health by BSAVA and WSAVA [5].

Precise definitions of obesity are debated, but it is generally accepted that obesity is the unhealthy accumulation of body fat. What is not controversial is that obesity is a consequence of energy intake chronically exceeding energy output. Consequently, obesity has commonly been considered a ramification of poor self-control in people or of inept management by animal owners. However, considerable evidence now shows obesity is better regarded as a disease of disordered energy homeostasis in which a multitude of genetic and environmental factors can contribute to increasing body fat.

This review first examines the pathophysiology of obesity and the role of genetics in the disease, focussing on the wealth of evidence from human and rodent studies. We then review the genetics of obesity and related metabolic disease in companion animals. Finally, we consider opportunities

for future research in companion animals that may improve understanding of both animal and human obesity.

#### **2. Factors Contributing to Obesity**

Recent, rapid increases in the prevalence of obesity have been caused by changes in activity and diet in the human population, and the same is likely true in animal populations [6,7]. However, it is clear that, although much of the human population is exposed to an "obesogenic environment" (with ready access to high calorie food and increasingly sedentary lifestyles), only a subset become overweight or obese. It is now well established that multiple factors including socioeconomic status, education level, and genetics are associated with whether a person is likely to develop obesity [8]. In companion animals, similarly diverse risk factors have been identified, with biological factors such as age, sex, and breed recognised alongside owner management factors in dog, cat, and horse obesity [9–12].

Acknowledgement of the multiple obesity risk factors is important, because it informs efforts to reduce obesity. Stigma against overweight people and parents of overweight children is well recognised [13,14], and the same is true for owners of overweight companion animals [15,16]. Such stigmas arise from the widespread view that weight gain is due to lax efforts to regulate food intake and exercise, either by a person, parent, or animal owner. It is important to acknowledge the risk factors beyond an individual's control (such as genetics) in order to improve the effectiveness of obesity prevention and treatment programmes.

#### *Obesity Susceptibility Is Highly Heritable*

Humans within the same environment, be that energy surplus or energy scarcity, display a highly heritable variance in body condition [17]. A wealth of data from twin and adoption studies, bolstered by later estimates of chip heritability in the era of high-density genotyping, supports that human obesity, indicated by body mass index (BMI), is a highly heritable trait. Heritability estimates range from 71–81% [17–19]. However, despite intensive efforts, the genes and mutations responsible for the majority of this heritability remains to be elucidated.

#### **3. Studies of Monogenic Obesity Have Been Highly Informative**

Obesity is usually a complex trait, with many genomic loci contributing incrementally to modulate an individual's susceptibility. However, monogenic forms of obesity also exist with patients usually coming to attention, because they develop severe obesity early in life. Interrogating the genetics of these rare patients, in combination with research in rodent and cellular models, has been hugely informative in elucidating the molecular basis of the regulation of energy homeostasis and body weight. Early studies of patients with monogenic forms of obesity focussed on candidate genes chosen based on information from rodent models of obesity [6].

#### *3.1. The Discovery of Leptin*

In 1950, a rodent model demonstrating severe obesity was identified, and the gene responsible named the obesity (*ob*) gene [20]. A similar obesity phenotype due to a different gene was discovered in diabetic mice, named the diabetes (*db*) gene [21,22]. Subsequent parabiosis experiments in which mice of contrasting genotype were surgically joined to share a circulation led to the conclusion that *ob*/*ob* mice lacked a circulating factor that controls eating behaviour, whereas *db*/*db* mice possessed such a factor but were not able to respond to it [23,24]. In 1994, the circulating factor was pinpointed and its function delineated; a hormone called leptin, which is secreted from fat cells [25]. The *ob* and *db* genes were subsequently identified as the genes for leptin (*Lep*) and its receptor (*Lepr*), respectively.

Shortly afterwards, a frameshift mutation in the human leptin gene (*LEP*) [26] was identified in children with severe, early-onset obesity. The mutation caused congenital leptin deficiency that was successfully treated with recombinant leptin therapy (See Figure 1). Since this, many more severely obese individuals have been identified with mutations in *LEP* [27–29] and provided help with recombinant leptin. Such studies clearly justify the importance of genetic research and its translational significance for prevention and treatment of disease.

Patients with mutation in the human leptin receptor gene (*LEPR*) [30–32] who display extreme obesity have also been identified. However, recombinant leptin treatment in these patients is entirely ineffective, since they suffer from leptin resistance as opposed to leptin deficiency [33–35].

**Figure 1.** Effect of recombinant leptin treatment on child with congenital leptin deficiency. Photographs of a 3-year-old child before leptin treatment weighing 42 kg (**left**) and the same child weighing 32 kg (**right**) after 4 years of treatment with recombinant leptin therapy. Figure from Farooqi et al. [36].

#### *3.2. The Leptin–Melanocortin Pathway*

We now understand that the primary effector mechanism for leptin's action is in the central nervous system (CNS), where it activates the hypothalamic leptin–melanocortin pathway [37], a neuroendocrine signalling mechanism, which transmits a signal about the status of the body's energy reserves to the brain and translates it into effector signals to promote optimal energy balance. In summary (See Figure 2), the hormone leptin is produced and released from adipose tissue. In the insulin-dependent fed state, insulin stimulates leptin release, and in greater amounts when energy reserves (in the form of fat stored in adipocytes) are larger [38,39]. In the brain, leptin acts on receptors in the arcuate nucleus (ARC) of the hypothalamus to activate proopiomelanocortin (POMC) neurons to produce the pre-pro-protein proopiomelanocortin (POMC). POMC undergoes proteolytic cleavage to produce a number of neuroactive peptides, the most important of which in regulating energy homeostasis are α- and β- melanocyte stimulating hormone (α-MSH, β-MSH). MSH peptides act on melanocortin receptor 4 (MC4R) expressed on second-order neurons of the paraventricular nucleus (PVN) of the hypothalamus, resulting in reduction in food intake and modified energy metabolism by interaction with multiple other pathways [17,40].

**Figure 2.** The leptin–melanocortin pathway. Simplified schematic of the leptin–melanocortin signalling pathway, which has a critical role in energy homeostasis by acting as a nexus through which information about energy status in the periphery can be relayed to the central nervous system (CNS) and integrated to control food intake and energy expenditure. ARC, arcuate nucleus of the hypothalamus; POMC, proopiomelanocortin; LepR, leptin receptor; PVN, paraventricular nucleus of the hypothalamus; MSH, melanocyte-stimulating hormone (α-MSH, β-MSH, γ-MSH); MC4R, melanocortin-4 receptor; SIM1, single-minded 1.

#### *3.3. Disruption of Leptin–Melanocortin Signalling Leads to Obesity*

Mutations that disrupt the leptin–melanocortin pathway are associated with severe obesity in both rodents and humans [41,42]. Dominant mutations in the proopiomelanocortin gene (*POMC*) are reported to cause severe, early-onset obesity in affected patients, often with other neuroendocrine features consistent with the diverse physiological roles of the multiple POMC-derived peptides [43–46]. Similarly, obesity has been reported due to mutations affecting prohormone convertase 1 gene (*PCSK1*), one of the enzymes responsible for proteolytic cleavage of POMC to its neuroactive derivatives, either singly or as part the more complex genetic condition Prader–Willi Syndrome [47].

Variants residing in the gene *SIM1* have also been associated with severe obesity and Prader–Willi-like syndromes [48,49], an effect attributed to this transcription factor's integral role in development of the PVN, a hypothalamic nucleus, which is most notable as the major site of MC4R expression [48].

More common are mutations in the *MC4R* gene, which have been shown to cause both dominant and recessive forms of monogenic obesity [50–53] and are responsible for up to 6% of severe, early-onset obesity cases [46,54–56]. More recently, *MC4R* variants with less severe effects on receptor function have been shown to be major modifiers of obesity risk in the wider human population [57–59]. Notably, those data show that the genetic background against which the *MC4R* mutation occurs has a large influence on the penetrance of the obesity phenotype.

#### *3.4. Other Causes of Monogenic Obesity*

Variants affecting other biological pathways have also been identified as monogenic causes of obesity, with the majority related to CNS regulation of energy metabolism. Semaphorin 3 gene (*SEMA3*) variants are rare causes of severe early-onset obesity and affect energy balance through their role in melanocortin neuron development [60]. Brain-derived neurotrophic factor (BDNF) acts on its receptor

tropomyosin receptor kinase B (TRKB), and there is increasing evidence that this signalling plays a significant role in sustaining equilibrium of energy balance in the brain [61]. Mutations affecting the concordant protein-coding genes (*BDNF*/*TRKB*) have been reported as causes of monogenic human obesity.

It is noteworthy that all the aforementioned variants primarily modify eating behaviour. The exception to this came with the discovery that rare variants in the Kinase suppressor of Ras 2 gene (*KSR2*) appear to cause obesity by affecting both energy intake and energy expenditure (although primarily via affecting the central control of energy metabolism) [62]. Reminiscent of this, a *CREBRF* variant that is common in Samoans but very rare in other populations, appears to be a major modifier of obesity risk by altering energy use in the body [63].

#### **4. Common Human Obesity**

Although cases of monogenic obesity exist, they account for only 5–7% of severe obesity cases [17, 54,64] and much less of the obesity that develops later in life. Common obesity is a complex trait, caused by the additive effect of hundreds, possibly even thousands, of common genetic variants [65]. Genome-wide association studies (GWAS) have provided valuable insight into which loci, genes, and variants are responsible. Large, consortium-based studies involving hundreds of thousands of human subjects have been performed on quantitative traits such as BMI (the best available indicator of body fat percentage for large scale studies) and obesity-related traits such as waist-to-hip ratio (WHR, an indicator of where an individual's body fat is stored) and measures of insulin resistance (and related metabolites), leading to identification of hundreds of quantitative trait loci [66] for obesity and related traits across the genome [66–71].

Despite those successes, just 3–5% of obesity heritability is explained by existing GWAS data [72]. The "missing heritability" of obesity [72,73] is hypothesised to be due to large effect rare variants yet to be identified: many loci of small/moderate effect too common to find with GWAS; non-additive genetic effects; and copy number variants (CNVs); among others [73–76]. Of those, the large-effect rare variant is particularly relevant to those studying veterinary species—selective breeding may lead to enrichment of variants of large effect within a breed, which are otherwise rare across the species overall [77]. Notably, an expanding body of evidence implicates the microbiome in obesity pathophysiology, and obesity-associated gut microbiota populations have been shown to be heritable [78]. Thus, an individual's microbiome may also account for a fraction of the heritable component.

Nonetheless, identified obesity-associated loci can be used to generate estimates of individuals' risk of developing obesity, known as polygenic risk scores (PRS), whereby a weighted effect score is generated as the product of allele count and effect score for risk variants. PRS have been proposed as clinical tools although their application is currently limited [79,80], in part because there is evidence that their validity may be limited in ethnic groups other than that in which they were originally derived. Nonetheless, PRS have been shown to be effective measures with which to stratify genetic obesity risk across a population, as in Figure 3. They have also improved our understanding of how background polygenic risk can alter trait penetrance in the presence of mutations with moderate–large effect on the trait, such as in *MC4R* [57], mentioned above.

#### *From GWAS to Function*

Notably, a large proportion of genes implicated by human obesity GWAS are preferentially expressed in the CNS [67,81], providing further evidence for the critical role of the brain in controlling energy homeostasis. Many genes implicated in monogenic obesity have also been flagged in these large-scale obesity GWAS, including *LEPR*, *MC4R*, and *POMC*, although they are usually not the most prominent findings [67,82].

Moving from GWAS-associated locus to identify causative variants and their mechanisms of action is difficult. Some of the most common obesity-associated loci identified lack clear mechanisms by which variants exert their effects. A notorious example of this is the repeated and robust association of

BMI with genetic variants lying within an intron of FTO α-ketoglutarate-dependent dioxygenase gene (*FTO*) also known as fat-mass- and obesity-associated gene (in mice). Those findings led to extensive efforts to characterise the role of this gene and provide evidence that *FTO* plays a part in regulating food intake. However, it later transpired that the focus on the closest gene was inappropriate and that obesity association for this locus is, at least in large part, due to altered regulation of a neighbouring gene, iroquois homeobox 3 (*IRX3*), which has an impact on peripheral adipocyte metabolism [83]. Although those findings have been challenged [84,85], it remains an example of how complex it has been to move from BMI locus to novel obesity biology.

**Figure 3.** Histogram of cumulative effect of BMI risk alleles in a large human GWAS. Mean BMI for each bin is shown by the black dots (with standard deviation), corresponding to the right-hand y axis and is compared to simple PRS, based on unweighted risk allele counts. These data demonstrate how the cumulative effect of multiple polygenic risk alleles is strongly associated with BMI across the human population. Figure from Locke et al. [67].

#### **5. Genetic Insight into Obesity Comorbidities and Metabolic Syndrome**

Obesity exerts a huge toll on human health due to its attendant comorbidities including type 2 diabetes mellitus (T2D), cardiovascular disease (CVD), stroke, and non-alcoholic fatty liver disease (NAFLD), among others [86–88]. Development of these comorbidities is largely mediated by a constellation of clinical features commonly referred to as "metabolic syndrome" that include raised blood pressure, dyslipidaemia, increased triglycerides, and insulin resistance [89,90].

However, it is recognised that the severity of obesity does not necessarily correlate with the onset or severity of metabolic syndrome, and that some people remain metabolically healthy despite severe obesity [91–93]. Those differences occur not just between individuals but also between ethnic groups, leading to different recommended cut-offs for what is regarded as a "healthy" BMI in different populations. For example, a WHO expert consultation considered alternative BMI "cut-offs" in Asian populations given their high risk of T2D and CVD at BMIs considerably lower than previously recommended [94].

How obesity leads to its comorbidities has been intensively studied in humans and model organisms, informed by genetics. For instance, the relationship between obesity and T2D is reinforced by the presence of vertical pleiotropy with shared loci implicated in GWAS for both T2D and BMI. The causal relationship between obesity and its complications has been confirmed by Mendelian randomisation studies [95]. Similar suggestions and confirmation of cause and effect come from GWAS and Mendelian randomisation studies of BMI/insulin resistance and NAFLD [96].

More nuanced mechanistic insight has come from the study of rare patients with severe, early-onset, monogenic forms of insulin resistance. Studies of rare patients with mutations affecting the insulin receptor (INSR) or which interrupt signalling downstream of it have provided insight into the pathogenesis of insulin resistance and consequent dyslipidaemia or hepatic lipid accumulation in

common obesity [93]. For instance, common variants in the gene encoding insulin receptor substrate 1 (*IRS1*), an important component of the intracellular signalling pathway, are associated with T2D, most likely by making affected individuals more susceptible to developing obesity-associated insulin resistance [97].

Similar value came from the study of patients with lipodystrophy, a rare condition in which there is complete or partial absence of adipose tissue. Recognition that the absence of fat was associated with early and severe insulin resistance and features of metabolic syndrome was central to establishing that the ability to store fat is critical to maintenance of normal metabolic function [98]. Those conclusions were later supported by analysis of genetic determinants of insulin resistance in the wider population, which provided evidence that variation in susceptibility to obesity co-morbidities are in part due to variation in individuals' ability to develop and maintain healthy fat storage depots [70].

#### *Overview—Molecular Mechanisms Underlying Obesity Co-Morbidities*

Together, human genetic studies and those in cell and animal models mean we now understand much more about how obesity causes disease. Much of those data support the theory of "limited adipose expandability". In brief summary, it appears that the body can accumulate fat in a healthy manner only until it reaches the limit of adipose tissue's ability to expand. Subsequently, adipose tissue dysfunction leads to inflammation and increased lipid flux, which in turn leads to ectopic lipid accumulation in non-adipose organs such as the liver and muscles. Combined with local and systemic inflammation, this leads to widespread insulin resistance, which in turn, results in increased insulin demand as the body attempts to maintain glucose homeostasis. Initially, insulin-producing β cells in the pancreas meet this demand but ongoing insulin resistance and ectopic lipid deposition in the pancreas itself can ultimately lead to β cell failure, hyperglycaemia, and T2D. Simultaneously, the altered lipid flux leads to dyslipidaemia and multiple mechanisms converge to cause hypertension [99,100].

Whether an individual develops obesity-associated co-morbidities may therefore be influenced by genes at many levels: a genetic predisposition to weight gain increases the chance of obesity developing; the ability to develop and maintain healthy fat stores appears to differ between individuals and across ethnic groups; if variants that perturb insulin signalling are present, they may promote the development of insulin resistance; and genes may also influence the ability of β cells to respond with adequate insulin where resistance develops. Consequently, in common disease, it is the net effect of many variants affecting such systems, combined with environmental influences, which determines whether obesity-related comorbidities develop.

#### **6. Applying Current Knowledge to Study Companion Animal Disease**

It is well established that obesity and its sequelae are heritable traits in humans and laboratory animals. The majority of obesity genes affect central regulation of food intake but notable examples whereby genes affect energy expenditure do exist. Genetic variance in susceptibility to obesity co-morbidities is well established and can be due to effects on a variety of body systems and processes.

Importantly, this wealth of information is available to inform investigations into companion animal disease. However, there is much yet to learn, and it is equally true that study of spontaneously occurring, companion animal models of obesity is a potentially valuable way to learn new biology relevant to all species. To date, studies of the molecular basis of obesity and associated disease remain limited in companion animals, with the predominance of small-scale and candidate genes studies somewhat limiting their value.

Available evidence suggests much commonality in the pathophysiology of obesity between humans, dogs, cats, and horses, but there are important differences in the incidence and severity of obesity-related co-morbidities between companion animal species, which likely reflect fundamental differences in physiology. Below, we consider in sequence the studies performed in dogs, cats, and horses to investigate the genetic basis of weight gain and, where relevant, obesity-related disease.

#### **7. Canine Obesity Genetics**

An estimated 34–59% of pet dogs in developed countries are overweight or obese [3,9,101–104]. Multiple comorbidities are associated with canine obesity, most notably: orthopaedic disease, exacerbation of breathing difficulties, and urinary incontinence, along with significantly decreased life expectancy [105–108]. Canine obesity is associated with increased blood pressure, dyslipidaemia, and insulin resistance, and this has been characterised as "obesity-related metabolic dysfunction" [109,110]. There is also some epidemiological evidence that obesity is a risk factor for developing canine, but this is not a recognised problem in clinical practice diabetes [111]. Other human obesity-related complications such as CVD and NAFLD do not have equivalent recognised presentations in canine practice. As with humans, studies suggest that weight loss leads to normalisation of many of the associated metabolic perturbations in dogs [109].

Evidence for the role of genetics in governing canine obesity susceptibility comes from the fact that obesity prevalence differs dependent on breed, with some breeds predisposed (e.g., Labrador retrievers, pugs, golden retrievers) and others resistant (e.g., greyhounds, whippets) [10,112,113]. Concordantly, differences in eating behaviour and food preference between breeds have been shown [114,115].

Although some of the variation in obesity susceptibility between breeds could be down to differences in "fashion" for each breed, the repeated finding of breed as a risk factor for obesity means genetic determinants are more likely the cause. Such variation in trait susceptibility between breeds is common in dogs due to the species' unusual population architecture in which there is high diversity as a whole combined with high homogeneity within breeds, the product of population bottlenecks at breed formation and subsequent intensive selection for breed-specific traits [116–118]. Heritability of canine obesity is yet to be determined, but some have speculated that it will be reminiscent of high estimates in humans and other animal species [119–122].

The genetic basis of body mass has been investigated using GWAS in dogs, but those results cannot be regarded as indicative of obesity because of the wide variability in body morphology breeds. Instead, those study results likely better reflect the combination of height, length, musculature, and limb: trunk ratio [123].

To date, no GWAS of canine fat mass or related traits have been reported. In part, this may be related to the difficulty of accurately phenotyping obesity in dogs; the aforementioned wide variety of skeletal size and shape in dogs precludes using simple and commonly available measures such as weight to estimate fat mass. Quantification of adiposity can be done using sophisticated imaging or biochemical techniques, electrical impedance or ultrasound measurement of fat deposits, but those cannot practically be applied in large cohorts [124]. The clinically widely used and well-validated alternative is to assign dogs a "body condition score" (BCS). This should be performed by a veterinary professional who uses a variety of haptic and visual clues to assign dogs to one of a series of ordinal scores. The most widely used 9-point scale [125] has been thoroughly validated [124,126] such that each one point increase in BCS above "normal" equates to approximately 10% increase in fat mass. A 5-point scale is also available and has been reported in some of the studies below.

#### *7.1. Genes Investigated in Canine Obesity*

Several candidate gene studies have been performed in domestic dogs and closely related, farmed *Canidae* species. The merit of looking at genetic factors in other *Canidae* species is arguable, but they are reported here, since previous studies have found them relevant. The best characterised of the later mentioned studies are the ones that have considered genes implicated in the leptin–melanocortin signalling pathway. All gene variants discussed can be found in Supplementary Table S1.

#### 7.1.1. POMC

As discussed above, POMC has an integral role in leptin–melanocortin signalling, a neuroendocrine pathway highly conserved across species. Raffan et al. [7] identified various canine *POMC* variants including a 14 bp deletion at position 17:19431807-19431821, present in Labrador retrievers and flat-coated retrievers (FCRs) but absent from a wide variety of other breeds tested. The deletion was shown to be a major modifier of weight, adiposity (measured as BCS on a 9-point scale), and appetite (measured using a validated, owner-reported measure of eating behaviour [114]) in ≥210 Labrador retrievers, with findings for weight and appetite replicated in ≥196 FCRs. Appropriate confounders were accounted for by adjustment for age, sex, neuter status, and colour, thus association results can be confirmed with high confidence. Notably, the mutation was more common in Labrador retrievers working as assistance dogs compared to the pet population. For each additional allele carried by either breed, there was an approximately 0.5 BCS/2 kg increase in adiposity/body mass, and a similar incremental increase in food motivation.

The authors showed the *POMC* deletion results in altered production of β-MSH and β-endorphin, two neuroactive peptides derived from *POMC*, and confirmed that canine MC4R receptors have similar affinity for and response to α-MSH and β-MSH as the equivalent human pairings. Whilst the reduced action of β-MSH on MC4R in the leptin–melanocortin signalling pathway is the most likely mechanism resulting in this phenotype, β-endorphin may also be implicated.

Subsequent studies have confirmed the presence of this mutation [120,127,128]. In the British Labrador retriever population, Davison et al. [127] found there was no association for the deletion with occurrence of diabetes mellitus (DM) in the breed (a finding that is not surprising given that obesity-associated insulin resistance is not clinically recognised as predisposing dogs to diabetes, despite there being some epidemiological evidence that it may play a role population-wide).

The data had relevance to those studying human obesity biology, because canine *POMC* is highly homologous to human *POMC*. The same cannot be said for rodents that lack a proteolytic cleavage site in the pro-protein and so do not produce β-MSH. Since humans with mutations causing β-MSH absence are very rare [129,130], the canine model in this case provided important insight into the role of this particular peptide in energy homeostasis, as well as establishing a clear role for genetics in governing obesity susceptibility in dogs.

#### 7.1.2. MC4R

Canine *MC4R* gene polymorphisms have been identified in several studies, with some attempting to analyse their effects on obesity-related phenotypes. *MC4R* is a suitable candidate gene due to MC4R's central position in the leptin–melanocortin pathway and given *MC4R* variants are among the most common cause of monogenic human obesity. Humans with *MC4R* mutations are also modestly taller than unaffected individuals of the same age and sex, or similarly obese unrelated subjects without *MC4R* mutations [53,131]. The majority of the identified variants lack functional data for characterisation of their consequences, but Yan et al. [132] reported that the missense variant c.637G>T (p.Val213Phe, rs852614811) had no significant effect on cAMP-mediated signalling downstream of the receptor.

Skorczyk et al. [133] conducted a study on canid *MC4R*, mapping and characterising it for the first time. The study cohort incorporated 31 dogs of 19 breeds and 35 farmed red foxes (*Vulpes vulpes*). They identified three polymorphic sites in the canine *MC4R* gene and four in the red fox. The three canine *MC4R* polymorphisms identified were a non-synonymous variant (c.637G>T, p.Val213Phe, rs852614811); a synonymous coding variant (c.777T>C, p.Ala259Ala, rs851987283); and a 3 UTR variant (c.\*33C>G, rs851539399). No association testing was conducted.

Van den Berg et al. [134] studied *MC4R* in 32 golden retrievers. They identified the same polymorphic sites as Skorczyk et al. [133] (c.637G>T, p.Val213Phe, rs852614811; c.777T>C, p.Ala259Ala, rs851987283; c.\*33C>G, rs851539399), plus an additional synonymous coding variant (c.868C>T, p.Leu290Leu, rs851062983). After linkage disequilibrium (LD) filtering, they performed an association study for the three remaining polymorphisms (c.637G>T, c.777T>C, c.868C>T) in a larger cohort of golden retrievers (*n* = 187). Several morphological measures were tested: weight, length, height, and body index score (BIS), but not BCS. For the association, they appropriately corrected for sex, age, and

polygenic effects. However, they do not mention any account for neuter status (a well-recognised risk factor), and there was no measure of relatedness in the sample. No statistically significant association was found for any of the phenotypes. Whilst this may be genuine, the study was underpowered to find a small effect size, and these variants may warrant further investigation.

Zeng et al. [135] conducted an *MC4R* candidate gene study in beagles from a research colony. By sequencing two beagle dogs, they identified the previously reported synonymous coding mutation (c.777T>C, p.Ala259Ala, rs851987283) and a novel missense substitution (c.302C>A, p.Thr101Asn). (Note that the way this mutation is referred to in the original paper and a subsequent review article is somewhat confused. Zeng et al. [135] mis-name the mutations as C895T (for c.777T>C) and A420C (for c.302C>A). The missense mutation (c.302C>A) results in an amino acid substitution from threonine to asparagine (p.Thr101Asn = T101N), which can be clarified by figures within the paper. However, they misname this as N101T. Mankowska et al. [136], in a subsequent review, note the mistake but further misname it c.301A>C (p.Asn101Thr). The authors' (unpublished) capillary sequencing data means we can confirm the correct annotation as c.302C>A, p.Asn101Thr. Variant nomenclature corrections can be found in Supplementary Table S1). In 120 beagles, c.302C>A was significantly associated with body weight (*p* < 0.05). In contrast, c.777T>C was only significantly associated with body weight in the heterozygous (but not homozygous) state and in bitches only. Whilst it is plausible at least one of the single nucleotide polymorphisms (SNP) are associated with adiposity, there are important limitations to the study. Specifically, weight does not equate to adiposity but could indicate changes in height, length, or muscle mass. There was no adjustment for confounding factors such as age, sex, or neuter status, and no assessment of whether the variants were in LD and therefore whether these associations constitute independent signals.

Raffan et al. [7] sequenced the *MC4R* region in 33 Labrador retrievers of which 15 were obese and 18 were lean. Two novel *MC4R* variants were identified (c.989G>T, p.Ser330Ile; c.\*227C>T) plus three previously reported variants (c.637G>T, p.Val213Phe, rs852614811; c.\*33C>G, rs851539399; c.777T>C, p.Ala259Ala, rs851987283). None were distributed differently between the small lean and obese groups, and the variants were not pursued further. This group also sequenced *AGRP* in this cohort, another valid candidate gene which codes for agouti-related protein—a neuropeptide known to modulate food intake in the ARC [137]. No *AGRP* variants were identified.

Mankowska et al. [136] investigated *MC4R* variants in 270 dogs of four breeds in which they identified six known polymorphic sites (c.637G>T, p.Val213Phe, rs852614811; c.777T>C, p.Ala259Ala, rs851987283; c.868C>T, p.Leu290Leu, rs851062983; c.\*33C>G, rs851539399; c.\*227C>T; and c.-435T>C, rs852471376). Of the 270 dogs, they had full phenotypic data for 164. They concluded that none of the identified variants displayed differential association with BCS (5-point scale) or body weight. The study population was dominated by Labrador retrievers (*n* = 187) and no potential confounding factors were accounted for.

#### 7.1.3. FTO

The top association signal on most human GWAS for BMI, *FTO*, is of debatable merit as a candidate gene for canine obesity, due to the controversy as to whether *FTO* itself or neighbouring genes are the true effector pathway causing the association signal (see above). Grzes et al. [138] investigated *FTO* SNPs in four *Canidae* species including the dog. In 39 dogs of 14 breeds, sequencing identified six polymorphic sites in *FTO*: one missense variant (c.23C>T, p.Thr8Met, rs852870212 described as "23 C/T, Thr1Met"), two intronic variants, and three 3 flanking variants. Presence of the missense mutation differed by dog breed, but the small number of dogs genotyped meant association tests were (appropriately) not performed. However, of seven *FTO* polymorphisms in the red fox, one (a 5 flanking region variant) was tested for association with body weight and pelt weight (which the authors report is affected by subcutaneous fat mass) in a larger cohort of 390 red foxes. No significant association was identified.

Grzemski et al. [139] used targeted next generation sequencing (NGS) in 32 Labrador retrievers in a region incorporating *FTO* and *IRX3*. Several polymorphisms were identified and tested for association with BCS (5-point scale) in a larger Labrador retriever cohort (*n* = 165), suitably adjusted for sex, age, and multiple testing. No association was found. The authors reported 56% nucleotide identity between human and canine *FTO* and 72% for *IRX3*. Additionally, they compared methylation status of CpG islands between lean and obese dogs in a smaller cohort (*n* = 28)—no differences were found.

#### 7.1.4. MC3R

Skorczyk et al. [140] conducted a candidate gene study on the gene encoding melanocortin receptor 3 (*MC3R)*. The *MC3R* gene is a reasonable candidate, although it is notable that rodent and human phenotypes associated with *MC3R* variants are more subtle than for the related *MC4R,* with *MC3R* most likely affecting the maintenance of body mass within the homeostatically controlled upper and lower limits [141–144].

The study cohort used by Skorczyk et al. [140] consists of four canid species: 31 dogs of 19 breeds, 35 red foxes, 30 arctic foxes (*Vulpes lagopus*/*Aloplex lagopus*), and 30 Chinese raccoon dogs (*Nyctereutes procyonoides procyonoides*). Multiple polymorphisms were found in *MC3R* in three of the species; no variants were identified in *MC3R* of the arctic fox. Two polymorphic sites were identified in the dog: short tandem repeat (STR) in 5' flanking region (c.-90delT, rs853092001) and a synonymous substitution (c.142C>T, p.Leu48Leu, rs8916554). Association analysis for *MC3R* variants was then performed in cohort of 376 male red foxes. Two SNPs were in LD and both were associated with a small, statistically significant increase in body mass. There was limited information provided on potential environmental confounders and population stratification, with no apparent adjustment. It therefore cannot be ruled out that any associations observed may be as a result of these factors.

#### 7.1.5. INSIG2

Insulin-induced gene 2 protein (INSIG2) has a role in lipid metabolism and has been linked in human GWAS studies to various measures of circulating blood lipids. Its coding gene (*INSIG2*) has been implicated in human GWAS for BMI [145], but that association has not been reliably replicable [146], so its validity as a candidate obesity gene is somewhat limited. Grzes et al. [138] analysed polymorphic sites in *INSIG2* in four *Canidae* species including the dog. In 32 dogs of 14 breeds, seven polymorphic sites were identified: two 5 UTR variants referred to as 5 flanking region variants, (c.-90G>A (referred to as–91 G/A), rs852813691; c.-1C>T, rs852335828), one missense variant (c.40C>A, p.Arg14Ser, rs850773724), and four intronic variants. No association test was conducted in dogs, but one intronic SNP identified in the red fox was significantly associated with pelt weight in 390 red foxes. Skin weight may be influenced by subcutaneous fat mass, but this is not a valid measure of adiposity and such findings in the red fox may not be translatable to dogs.

#### 7.1.6. GPR120/FFAR4

G-protein coupled receptor 120 (GPR120), also known as free fatty acid receptor 4 (FFA4/FFAR4), functions as a receptor for unsaturated long-chain free fatty acids. The gene encoding this protein is known as *GPR120*/*FFA4*/*FFAR4* [147,148], and coding mutations in the gene have been associated with human obesity [149]. The following variants are described based on alignment to the canine *FFAR4* gene. Miyabe et al. [150] found nine *FFAR4* polymorphisms in a cohort of 141 dogs of 21 breeds: five synonymous substitutions (c.252C>G, p.Ala84Ala, rs852631320; c.282C>G, p.Pro94Pro (referred to as p.Asp94Asp), rs851850900; c.702A>G, p.Thr234Thr; c.726G>A, p.Thr242Thr; c.984T>C, p.Asn328Asn, rs852472019) and four non-synonymous substitutions (c.287T>G, p.Leu96Arg; c.307G>A, p.Ala103Thr; c.446G>C, p.Gly149Ala; c.595A>C, p.Thr199Pro, rs853030954 (referred to as c.595C>A, p.Thr199Pro)). SNPs were tested for association with BCS, the frequency of c.595A>C (referred to as c.595C>A) was significantly higher in dogs with a higher BCS. However, several confounding factors are unaccounted for and the study failed to correct for population stratification. Therefore, such identified associations

may have been due to confounders or, given the unequally represented multi breed cohort, due to unaccounted-for population stratification.

#### 7.1.7. PPARs

Peroxisome proliferator-activated receptors (PPARs) are a group of nuclear receptor proteins that function as transcription factors and have roles in the regulation of cellular differentiation, development, and metabolism [151,152]. Nishii et al. [153] sequenced the genes encoding PPARβ and PPARγ (*PPARB* and *PPARG*) in two dogs and observed tissue-specific expression of various PPARs. They also investigated the presence of polymorphic sites in *PPARG* and identified a single polymorphism. No association test with phenotype was conducted. Since a relatively small canine cohort were genotyped, it is possible that other *PPARB* polymorphisms were missed.

#### 7.1.8. Adipokines

Adipose tissue communicates with the rest of the body in part by release of a range of molecular signals known as adipokines. Adipokines include leptin, pro-inflammatory molecules such as TNFα and IL6, and the much-debated resistin, the latter three of which were investigated in dogs by Mankowska et al. [154]. The choice of these genes is a little surprising in that, although there is evidence that each has a role in the development of insulin resistance secondary to obesity, there is little evidence that variants promoting inflammation are causal in obesity [155] and considerable evidence to the contrary (https://www.ebi.ac.uk/gwas/) [156]. Nevertheless, in 77 dogs of 17 breeds, Mankowska et al. identified multiple variants, including 13 in *TNF*, four in *IL6*, and eight in *RETN*. Three of these variants were missense substitutions: one in *TNF* (c.548A>T, p.Glu183Val) and two in *RETN* (c.19C>T, p.Leu7Phe, rs852470997; c.236C>G, p.Ser79Cys, rs851766760—referred to in the paper as a synonymous coding variant). The five most common variants (*TNF*: c.-40A>C, rs22216187; c.233+14G>A; *IL6*: c.309+215T>C; *RETN*: c.194-69T>A, rs853182485; c.75G>A, p.Glu25Glu, rs852185407) were genotyped in 260 dogs and tested for association with BCS, using breed-specific sub-groups to do the association analysis, including an "others" category for the poorly represented breeds. No association was found for the *IL6* and *RETN* variants with BCS in any breed group. The two *TNF* SNPs (c.-40A>C and c.233+14G>A) were significantly associated with body condition in Labrador retrievers but not in any other breed group. Whilst those associations may be meaningful, the failure to include recognised risk factors such as age, sex, and neuter status or to detect or correct for population stratification in the sample could have affected these results.

#### **8. Feline Obesity and Associated Disease**

Estimates suggest 12–63% of pet cats are overweight or obese [2,157–162]. Feline obesity is associated with multiple comorbidities, most notably DM and hepatic lipidosis [12,158,163–167]. This arguably makes human and feline obesity comorbidities closer compared to dogs [168,169]. However, human comorbidities such as hypertension and atherosclerosis are not commonly observed in obese cats [107,170,171].

Obesity's relationship to DM in cats is well characterised with evidence suggesting a similar pathophysiology to human T2D [107,171–173] in which obesity leads to insulin resistance and ultimately β-cell dysfunction and diabetes [166,173–177]. Cats are therefore considered a suitable animal model for study of obesity associated T2D, with translational significance to humans. In contrast, whilst obesity is a well-recognised risk factor for feline hepatic lipidosis, the pathophysiological link between the two is less well characterised and subject to some debate [176,178,179].

#### *8.1. Evidence for the Role of Genetics in Feline Obesity and Related Disease*

Multiple studies suggest breed as a risk factor of feline obesity, with certain breeds displaying predisposition to becoming overweight/obese [12,180–182]. Although breed-specific obesity risk differs by study, domestic shorthaired cats (DSH) are consistently found to be at high risk, whilst longhaired breeds are generally at lower risk. Persian cats also have low obesity risk in most studies [12,181,182], but one study found them to be at high risk [180]. Together, these data suggest obesity may be at least in part a heritable trait in cats.

Notably, pet cats are made up of a majority of mixed breed cats (commonly referred to as DSH and domestic long-haired (DLH) breeds but are most commonly outbred and relatively genetically diverse) and a minority of pedigree cats, in which genetic architecture is reminiscent of dogs, with low within-breed diversity and evidence of population bottlenecks and genetic selection [183,184].

Until recently, genetic studies in cats have been stymied by the absence of a well-annotated, complete feline genome and lack of a commercial feline SNP genotyping array. Consequently, there are fewer genetic studies reported in cats overall, and none of the common forms of feline obesity. However, one group has reported candidate gene studies for feline DM [185], a related trait and one that is likely to be subject to similar genetic influences (including overlap with obesity predisposition) as described above for humans.

#### *8.2. Familial Obesity in a Feline Colony*

In a population of well-characterised, related research cats, a familial form of obesity has been reported [186]. Some cats displayed a clear predisposition to obesity, and segregation analysis suggested a single major gene was likely responsible. The report is more akin to human monogenic obesity or mutations of large effect against a variably obesogenic polygenic background.

A subsequent study found that the cats predisposed to being obese had a lower energy requirement and higher food intake than cats that did not tend to gain weight [187]. However, the energy requirement measurements were (for the obesity-prone group) performed not long after a period of restricted feeding and weight loss, interventions well recognised to cause reduced energy expenditure irrespective of baseline status [188]. In a subsequent generation of the same cohort, food intake and energy expenditure were investigated [189], and obesity-prone cats had higher food intake early in life but not lower energy expenditure.

In the same feline cohort, Keller et al. [190] investigated metabolic responses to different diets in cats predisposed to obesity vs. lean cats. No difference in metabolic response was found between the two groups. Additionally, in conference proceedings [191] the same group report attempts to map obesity genes were made and identified plausible candidate genes. Further comment is not possible given the scant information reported.

#### *8.3. Genetics of Diabetes Mellitus in Pet Cats*

Forcada et al. [185] investigated whether *MC4R* polymorphisms were associated with diabetes in DSH and Burmese cats. For each of a non-diabetic control group and a diabetic case group, there were 60 lean and 60 overweight cats, making 240 in total. The authors report that one *MC4R* polymorphism (c.92 C>T, p.Leu31Pro, rs783632116) was significantly more common in obese diabetic cats than obese non-diabetic cats, a finding not replicated in lean cats. Heterozygous and homozygous carriers in the non-diabetic subgroup were merged (assuming a dominant mode of inheritance). The authors do not report comparison of allele frequencies between lean and overweight cats (irrespective of DM status).

The authors speculate that this variant may act independently of an effect on body weight, but these reviewers suggest that to be a bold statement given the limitations of the data presented and the weight of evidence from other species about the role of MC4R in controlling food intake and obesity predisposition. For example, although *MC4R* is significantly associated with T2D in human GWAS, that association disappears when the analysis model corrects for BMI [192,193]. Consequently, it seems equally or more plausible that the association reported may in fact be a result of vertical pleiotropy rather than a direct effect.

The same group in conference abstracts report the results of a GWAS for feline DM in a cohort of 581 DSH cats [194,195], later adding Burmese cats [196]. However, the data are not reported in sufficient detail to reiterate here.

#### **9. Obesity and Related Metabolic Disease in Horses**

As in other companion animal species, equine obesity is common with an estimated 20–70% of horses overweight/obese [11,197–200]. This is a significant clinical problem, because obesity is a risk factor for the development of laminitis, a common, crippling disorder of the equine hoof. That association is thought to be mediated predominantly via a collection of risk factors known as equine metabolic syndrome (EMS). EMS was first described in 2002 [201], and the pathophysiological links between obesity and laminitis have been extensively studied since, although it is acknowledged that the current working understanding requires refinement [202].

EMS is defined by the presence of insulin dysregulation, characterised by clinical features including hyperinsulinaemia (either at baseline or in response to glucose challenge), hyperglycaemia, and/or evidence of peripheral insulin resistance [202]. Insulin resistance is commonly present in EMS, although there has been some debate as to whether that is always true and if alternative routes by which insulin dysregulation may develop exist [203]. Notably, although obesity and EMS are common, it is rare for horses and ponies to become diabetic [202].

Not all overweight equines develop EMS, nor does EMS always cause laminitis. Similarly, not all horses that have clinical features of EMS are overweight. Those paradoxes exist between individuals and across breeds, with some breeds apparently particularly prone to developing laminitis despite only moderate weight gain [198]. This is reminiscent of the situation described above in humans, whereby there is variability between individuals and ethnic groups concerning whether, and at what point, obesity-associated complications occur. It is clear, therefore, that a better understanding of equine obesity and its related conditions is required.

#### *9.1. Genetics Influence Equine Obesity, EMS, and Laminitis*

Evidence that genetics influence the development of equine obesity come from recognition that breed and "type" are clear risk factors with ponies at highest risk, followed by cob type breeds [11,204–208]. In Pura Raza Español horses, the heritability of BCS was found to be 14–24% [209]. Importantly, a within-breed study can only estimate the variance due to genetic variation within that (homogeneous) breed, so this heritability figure should not be generalised to the species as a whole. To date, there are no reports of GWAS or candidate gene studies in equine obesity.

The genetics of EMS have been more intensively interrogated. Breed is a well-recognised risk factor and breeds largely overlap with breeds at high obesity risk, unsurprising given the association between the conditions [202]. A similar variability between breeds has been demonstrated for insulin sensitivity and related biochemical parameters [210–213]. A recent study used genome-wide SNP data to estimate heritability of several traits known to be perturbed by EMS (glucose, insulin, measures of insulin sensitivity, and dyslipidaemia) and found they were moderately to highly heritable [214]. Again, such within-breed comparisons are valuable but are likely to underestimate the true heritability in the equine population.

Genetic studies of laminitis are worthy of report, given the common co-segregation of obesity, EMS, and laminitis. In an early study of crossbreed ponies, the authors concluded laminitis was a dominant trait in the pedigree with variable penetrance due to sex, age, and epigenetic-related variables [215]. Subsequent GWAS are mentioned below.

#### *9.2. GWAS for EMS and Related Traits*

Lewis et al. [216] performed a GWAS for laminitis and related traits in a population of 64 Arabian horses. A locus on chromosome 14 was associated with laminitis and insulin concentration. In a second cohort of 50 horses of the same breed, the identical phenotypes were not available, but the same region was associated with BCS and an alternative measure of insulin resistance. The closest gene was the poorly characterised *FAM174A*, which the authors sequenced. They found two closely linked variants, present in multiple breeds. The authors suggested an 11-guanine polymorphism near *FAM174A*

might have potential as a predictive test for horses at risk of obesity/EMS/laminitis. Subsequently, an Australian group found no association with metabolic traits for that marker in 20 (non-Arabian) ponies [217]. Similarly, a larger study in multiple breeds, including Arabians, failed to replicate the BCS, laminitis, or insulin resistance associations, although assuming a dominant model did identify a significant association with the adipokine adiponectin [218].

Selective breeding can lead to enrichment of alleles within a breed, which are scarce in the wider population, meaning it is plausible that a genetic variant may exist in Arabians at this locus, genetically linked to the 11-guanine allele but which has yet to be identified. Thus, a real finding in Arabians might not be replicated in other breeds. However, replication in a larger, independent cohort of Arabians would be advisable. There is no evidence that this allele would be a suitable genetic test for obesity, EMS, or laminitis predisposition in other breeds.

Norton et al. [219] recently performed a GWAS focussed on EMS in a larger equine cohort (*n* = 550) representing two high-risk breeds (Welsh ponies and Morgan horses). By collecting rich metabolic data, they were able to test for association with multiple "endophenotypes", more precise biochemical markers of insulin-related traits. Using endophenotypes may be more informative than performing a GWAS for a "convergent" phenotype such as laminitis, which can represent the clinical endpoint of multiple pathophysiological processes. The group appropriately adjusted for relevant confounding factors and for population stratification although not for BCS, meaning the results may identify loci associated with obesity rather than EMS per se. Hundreds of loci were identified across the multiple genome scans and the authors prioritised those which appeared to affect more than one breed.

The authors expressed surprise that there was not more overlap between breeds. That may be because major genetic determinants may be private to individual breeds, but it is perhaps more likely that the study was underpowered to find variants of small effect, or which were rare or invariate in one population. Protein coding genes in prioritised regions were enriched for involvement in pathways of inflammation, glucose metabolism, and lipid metabolism, all plausible as contributing to EMS pathology.

A final equine study considered the gene *HMGA2*, which had previously been identified as associated with short stature in pony breeds. In humans, short stature is a risk factor for insulin resistance. Concordantly, a study found pony breeds (which are smaller than horses) are more likely to get EMS [220]. The authors hypothesised a causal link between *HMGA2* and increased EMS susceptibility in ponies. They found a strong association with height in 264 Welsh ponies and a lesser association with several metabolic traits, which they reported as a pleiotropic effect.

#### **10. From Humans to Animals and Back Again**

#### *10.1. Lessons for Animal Genetics*

This review has summarised the wealth of research into the genetics of human obesity and its related metabolic perturbations, and the relative paucity of efforts to date in dogs, cats, and horses. Consequently, those studying veterinary species have much to learn from human geneticists and those studying laboratory animal models. By familiarising ourselves with those fields, there is much scope to fast-track animal studies to provide maximum insight into animal disease.

In particular, we note how careful delineation of clinical phenotypes has led to both genetic diagnoses and mechanistic insight in human patients, such as with familial partial lipodystrophy. Historically, such lean patients presenting with metabolic syndrome or T2D were considered inexplicable outliers. By recognising groups with shared clinical features, lipodystrophy has been recognised, patient care has improved, and we have a matured understanding of how healthy fat depots are essential to metabolic health [98]. Might similar clinical groupings exist, unrecognised, as underlying atypical presentations of EMS or determining why only some overweight cats become diabetic?

Human genetic epidemiological studies have, for some time, attempted to address the issue of cause and consequence between obesity and related pathologies. Performing GWAS not only for endpoints such as T2D or NAFLD, but also for those phenotypes corrected for BMI was a start in teasing out causal relationships between them, a process that has been improved with the advent of Mendelian randomisation [95]. Such approaches operate best at scale, which may limit their application in veterinary studies. Even so, an acute awareness that multiple routes can converge on a single phenotype (e.g., obesity, adipose dysfunction, and insulin signalling impairment converge to produce signs of metabolic syndrome) should inform our design and interpretation of veterinary GWAS.

Finally, human genetics also provides a lead in how to maximise the benefit of genetic studies. One notable contrast is that results of human genetic studies are rarely left "hanging"—mapped loci are further investigated, mouse models made to test the function of unknown genes, patients with specific genetic diagnoses are further studied, and precision treatments developed. This reflects a more established, larger, and better-funded research environment but provides a model for veterinary researchers to emulate.

Genetic findings have informed best practice to flag patients with uncommon presentations of common disease who may benefit from precision treatments [98,221]. In common, polygenic human disease, PRS are being used to counsel human patients about disease risk [80] and to stratify patients in research studies to understand better variable penetrance of variants of large effect [57]. In companion animals, there is a genetic test available for the retriever *POMC* mutation that can warn owners they have a dog at high risk of obesity and prompt them to effectively institute appropriate preventative measures [7], but other genetic findings in the field of obesity have yet to reach the veterinary clinic. Might we in future years detect cats with insulin signalling defects as candidates for insulin sensitising drugs in the pre-diabetic stage? Or might genetic profiling identify ponies at highest risk of laminitis? If so, genetic testing might prove a valuable clinical tool although one to be used only after robust validation of their utility in populations (e.g., breeds) different to that in which they were initially proposed.

#### *10.2. Lessons from Animal Genetics*

Fortunately, the flow of information between species need not be along a one-way street. Genetic studies of animal disease have already been informative to human research [7,222,223], and there is much potential to discover more. As genomic tools and better-annotated genomes become available, veterinary studies will be better able to provide insight into not just disease-associated loci but particular genes, mutations, and mechanisms too. That will clearly be a benefit to the species studied but also has the potential to benefit human health.

The current quality of genome builds and genomic tools available for companion animal studies vary by species. In dogs, the third-generation genome build (CanFam 3.1) has been available since 2014 [224], based on data from a Boxer dog, Tasha. Sequencing efforts integrating data from long read technologies produced genome sequences from a Great Dane (UMICH\_Zoey\_3.1, GCA\_005444595.1, PacBio RSII) [225] and a basenji, (Basenji\_breed-1.1, GCA\_004886185.2, Sequel), meaning there are now three potential reference sequences available in the species. Canine geneticists were also the first to drive production of commercial companion animal genotyping arrays, with the earliest containing 49k SNP markers [226], later increasing to 173k, 220k (CanineHD BeadChip, Illumina, San Diego, CA, USA), and 710k densities (CanineHD Array, Thermofisher, Waltham, MA, USA). Today, there are available arrays that cover >1 million genetic markers (Canine Genotyping Array, Thermofisher, Waltham, MA, USA), a subset of which were selected for use on the 460k/670k arrays (Canine Genotyping Array A/B, Thermofisher, Waltham, MA, USA). Additionally, high quality imputation from a lower to higher density SNP array level and to genome sequencing density has been successful in dogs [227,228].

In cats, the newest genome build (Felis\_catus\_9.0) was made available in 2017 [184], and one relatively low-density genotyping array of 63k SNP markers exists [229] but is not commercially available. No use of imputation software in cats has been reported. In horses, the most recent genome build EquCab3.0 has been available since 2018 [230], and the first equine array was reported in 2014 at a density of 54k (Equine BeadChip, Illumina, San Diego, CA, USA) [231]. Nowadays, two more SNP

arrays exist, both of which are commercially available, with marker densities of 65k (Equine BeadChip, Illumina, San Diego, CA, USA) and 670k (MNEc670K, Affymetrix, Santa Clara, CA, USA) [232]. Notably, a test equine array containing 2 million markers (MNEc2M) was created to inform the creation of MNEc670K and used successfully for imputation [232,233]. Imputation has also been performed in an equine population, from 54/65k density up to 670k [234].

Companion animal models of disease have the potential to "add value" to understanding broader biology for many reasons. Our animal companions share our homes and environments, spontaneously develop diseases similar to human conditions, for which they are often diagnosed and treated similarly but over a shorter time course. In some cases, disease processes that occur commonly in companion animals may be hard to study in humans. For instance, the retriever *POMC* mutation occurs in approximately 25% of Labrador retrievers, but the equivalent human mutations are very rare. The same molecule, β-MSH, is absent in rodents. Hence, in this situation, dogs provide an excellent animal model to shed light on a previously hard-to-illuminate area of human biology [7].

Similarly, cats are arguably more suitable than traditional animal models for the study of obesity associated T2D. Diabetes in cats bares closer resemblance to human T2D than that of the rodent model; it occurs naturally in cats, whereas T2D research in rodents is most often an induced disease state [235]. This means using the feline model of diabetes offers benefits for improved understanding of human T2D development, particularly for polygenic models [236,237].

Notably, there are advantages to studying such inbred populations, particularly dogs. Recent population bottlenecks at breed formation mean dogs have a very different genetic architecture to humans that makes complex trait mapping uniquely tractable in the species [238]. This means complex traits can be mapped in a breed with fewer individuals and fewer markers than in human populations. Dogs are therefore a compelling model for studying human metabolic disease [77,239–241]. Although, trait mapping in such inbred species means mapping to much larger loci than in humans [242,243], making causative variant identification within an associated locus more difficult.

Although cats do not display the same high level of LD as dogs, they have also been selectively bred and are proposed as potentially valuable models of human hereditary disease [183,244,245]. Horse genomes display lower levels of LD than both dogs and cats [183], but breed structure means there is potential for disease alleles being enriched and relatively easily studied within a discrete, definable population, meaning they too have potential to be valuable models of human disease [246].

#### **11. Conclusions**

The obesity epidemic is a major health concern in both human and companion animals, and there is a lot more to be discovered regarding the molecular basis of obesity and associated metabolic conditions. Despite clear evidence that obesity and related traits are highly heritable in companion animals, there are only limited studies to date investigating which genes are responsible and how they exert their effect. As this field matures, it promises tangible benefits for animal populations and, where considered as non-traditional animal models of obesity, has the potential to offer translational benefits too.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/11/1378/s1, Table S1: Candidate genes and gene variants for obesity and related metabolic disease in companion animals.

**Author Contributions:** N.W.: writing—original draft preparation, writing—review and editing; E.R.: writing—review and editing, funding acquisition, supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** E.R.'s post is funded by the Wellcome Trust Clinical Research Career Development Fellowship award 205187/Z/16/Z and N.W.'s post by the Biotechnology and Biological Sciences Research Council (BBSRC) DTP programme, BBSRC BB/M011194/1.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A** *CNTNAP1* **Missense Variant Is Associated with Canine Laryngeal Paralysis and Polyneuropathy**

**Anna Letko 1,\*,**†**, Katie M. Minor 2,**†**, Steven G. Friedenberg 3, G. Diane Shelton 4, Jill Pesayco Salvador 4, Paul J. J. Mandigers 5, Peter A. J. Leegwater 5, Paige A. Winkler 6, Simon M. Petersen-Jones 6, Bryden J. Stanley 6, Kari J. Ekenstedt 7, Gary S. Johnson 8, Liz Hansen 8, Vidhya Jagannathan 1, James R. Mickelson 2,**‡ **and Cord Drögemüller 1,**‡


Received: 28 October 2020; Accepted: 26 November 2020; Published: 27 November 2020

**Abstract:** Laryngeal paralysis associated with a generalized polyneuropathy (LPPN) most commonly exists in geriatric dogs from a variety of large and giant breeds. The purpose of this study was to discover the underlying genetic and molecular mechanisms in a younger-onset form of this neurodegenerative disease seen in two closely related giant dog breeds, the Leonberger and Saint Bernard. Neuropathology of an affected dog from each breed showed variable nerve fiber loss and scattered inappropriately thin myelinated fibers. Using across-breed genome-wide association, haplotype analysis, and whole-genome sequencing, we identified a missense variant in the *CNTNAP1* gene (c.2810G>A; p.Gly937Glu) in which homozygotes in both studied breeds are affected. *CNTNAP1* encodes a contactin-associated protein important for organization of myelinated axons. The herein described likely pathogenic *CNTNAP1* variant occurs in unrelated breeds at variable frequencies. Individual homozygous mutant LPPN-affected Labrador retrievers that were on average four years younger than dogs affected by geriatric onset laryngeal paralysis polyneuropathy could be explained by this variant. Pathologic changes in a Labrador retriever nerve biopsy from a homozygous mutant dog were similar to those of the Leonberger and Saint Bernard. The impact of this variant on health in English bulldogs and Irish terriers, two breeds with higher *CNTNAP1* variant allele frequencies, remains unclear. Pathogenic variants in *CNTNAP1* have previously been reported in human patients with lethal congenital contracture syndrome and hypomyelinating neuropathy, including vocal cord palsy and severe respiratory distress. This is the first report of contactin-associated LPPN in dogs characterized by a deleterious variant that most likely predates modern breed establishment.

**Keywords:** *Canis familiaris*; whole-genome sequencing; rare disease; contactin; neurological disorder; Leonberger; Saint Bernard; Labrador retriever

#### **1. Introduction**

Laryngeal paralysis (LP) can result from trauma or neoplasia involving the recurrent laryngeal nerves, peripheral nerve disease, or a primary or secondary disease affecting the muscle or neuromuscular junction. Loss of normal function of the larynx leads to breathing difficulties, reduced exercise and heat tolerance, as well as an increased risk of aspiration pneumonia [1]. Laryngeal nerve disease results in degeneration and atrophy of intrinsic laryngeal muscles followed by decreased or absent movement of the attendant laryngeal cartilages. During breathing, these cartilages control airflow into and out of the trachea. Affected dogs have stridor, may have a change in vocalization, and difficulty breathing due to the flaccid laryngeal vocal folds and corniculate processes of the arytenoid obstructing the lumen of the airway [1]. Normal laryngeal function protects the airway by closing off the lumen to prevent aspiration of food or water. In LP-affected dogs, the vocal folds remain in a paramedian position, causing airway resistance and turbulence, instead of abducting, as they normally would, to open the airway during inspiration. Frequently, affected dogs suffering from LP are treated by crico- or thyro-arytenoid laryngoplasty surgery, to improve breathing and, therefore, quality of life [2]. As the recurrent laryngeal nerve axons are some of the longest in the body [3], LP is often reported as part of a more generalized length-dependent polyneuropathy (PN) complex, which manifests with additional signs including proprioceptive and motor abnormalities, slowly progressing pelvic limb weakness, and loss of limb muscle mass [4].

Various mostly breed-specific canine inherited neuropathies form a heterogeneous group of degenerative diseases affecting motor and/or sensory and autonomic peripheral nerves. This group includes mixed forms of LP and PN [5], i.e., the laryngeal paralysis and polyneuropathy complex (LPPN), which has variable ages of onset among and across several dog breeds (OMIA 001206-9615, OMIA 001292-9615). Late-onset forms, e.g., geriatric onset laryngeal paralysis polyneuropathy (GOLPP), are also observed in various breeds including Labrador retrievers [6]. Leonberger dogs are known to be susceptible to LPPN; recently, a short list of potentially pathogenic variants for neurological disorders in this breed derived from whole-genome sequencing has been presented [7]. To date, variants in *ARHGEF10* [8] and *GJA9* [9] have already been associated with certain forms of the disorder and designated with breed-specific names Leonberger polyneuropathy type 1 (LPN1; OMIA 001917-9615) and Leonberger polyneuropathy type 2 (LPN2; OMIA 002119-9615), respectively. These two variants, however, do not explain all the phenotypically described cases in Leonbergers [7]. The *ARHGEF10* variant has also been reported in the related Saint Bernard breed, but again it did not explain all LPPN cases [8]. Alaskan huskies, black Russian terriers, and Rottweilers with PN including LP and respiratory distress are known to have deleterious variants in the *RAB3GAP1* gene, a member of the RAB3 protein family implicated in regulated exocytosis of neurotransmitters and hormones (OMIA 001970-9615) [10–12]. Another major risk factor for canine LP recently described in miniature bull terriers and bull terriers is a variant in the *RAPGEF6* gene encoding a widely expressed nucleotide exchange factor whose function is not well understood (OMIA 002222-9615) [13].

In general, there are limits to precisely diagnosing neurological diseases in dogs in the clinic. For example, in a previous study [14], we noticed Leonbergers that were initially clinically diagnosed as polyneuropathy-affected, although, in fact, they were suffering from leukoencephalomyelopathy, a juvenile-onset neurodegenerative disorder of the CNS white matter with distinctive pathological features, caused by a recessive variant in the *NAPEPLD* gene (OMIA 001788-9615).

Our aim in this study was to identify additional causative genetic variants associated with younger-onset laryngeal paralysis and polyneuropathy (LPPN), by focusing on two closely related giant dog breeds [15], namely the Leonberger and Saint Bernard.

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

All animal experiments were performed according to local regulations, and all animals in this study were examined with the consent of their owners. The study was approved under IACUC protocol 1903-36865A at the University of Minnesota, the Michigan State University Institutional Animal Care and Use Committee (AUF number 01/11-009-00), and by the Cantonal Committee for Animal Experiments (Canton of Bern; permit 71/19) at the University of Bern.

#### *2.2. Animal Selection*

Data on 15,378 dogs from 243 breeds, 321 dogs of mixed or unknown heritage, and 62 wild canids were collected in three different sets for this study (Table 1 and Supplementary Table S1). The discovery cohort included 426 Leonbergers either showing signs of LPPN with an age of onset ≤5 years or healthy control dogs at ≥8 years of age, and 91 Saint Bernards either showing signs of LPPN with an age of onset ≤5 years or population control dogs with genome-wide association study (GWAS) data available from unrelated studies [15,16]. All 517 Leonbergers and Saint Bernards were genotyped for the *ARHGEF10* variant [8], and the Leonbergers were also genotyped for the *GJA9* variant [9], in order to include only dogs homozygous wild type for these loci. In addition, the Leonbergers were genotyped for a previously reported leukoencephalopathy-associated *NAPEPLD* variant [14], in order to rule out another known underlying neurological disease with similar clinical phenotype.


**Table 1.** Number of canids in each studied cohort.

<sup>1</sup> Additional details including the source of all data are available in Supplementary Materials Table S1. <sup>2</sup> Includes affected dogs with previously described laryngeal paralysis and polyneuropathy (LPPN)-associated variants in *ARHGEF10* [8] and *GJA9* [9]. <sup>3</sup> Includes dogs homozygous for the previously described leukoencephalomyelopathyassociated variant in *NAPEPLD* [14]. <sup>4</sup> Includes seven mixed-breed dogs enrolled in the Michigan State University's GOLPP study [6].

A validation cohort used for targeted genotyping of the newly discovered variant consisted of 1070 dogs with known LPPN phenotypes (Table 1 and Supplementary Table S1). There was no age of onset restriction for the cases in the validation cohort. Included in the validation cohort were 193 Labrador retrievers and seven mixed-breed dogs from an ongoing geriatric onset laryngeal paralysis polyneuropathy (GOLPP) study; these were used as an independent validation group (Table 1).

Finally, a population cohort consisting of 14,112 dogs, 58 wolves, two golden jackals, one Andean fox, and one dhole (Table 1 and Supplementary Table S1), with no available information about their health status, was used to determine the absence/presence and frequency of the described variant-associated haplotype across canids.

The information about age of onset of clinical signs in the LPPN-affected dogs was available for a subset of 770 dogs from the discovery and validation cohorts with detailed health information from three breeds (596 Leonbergers, 28 Saint Bernards, and 146 Labrador retrievers) The statistical significance of the differences between groups was evaluated with Student's *t*-test and *p* < 0.05 was considered as significant.

#### *2.3. Sample Preparation*

Genomic DNA was isolated from EDTA blood samples, buccal swabs, or archived muscle biopsies by using either the Gentra PureGene kit (Qiagen, Hilden, Germany) or the Maxwell RSC Whole Blood DNA kit (Promega, Dübendorf, Switzerland).

Clinical cases of polyneuropathy and laryngeal paralysis were evaluated in three dogs homozygous for the studied *CNTNAP1* variant with available nerve biopsies; these included a 3-year-old Saint Bernard, a 3-year-old Leonberger, and a 9-year-old Labrador retriever (Supplementary Table S1). Three normal adult dog samples from Labrador retrievers were used as controls. The ages for the control dogs were 8–10 years of age. All the archived nerve specimens were obtained years prior to the identification of the *CNTNAP1* variant. Peroneal nerve specimens, pinned on cork discs to maintain length and orientation, were immersion-fixed in 2.5% glutaraldehyde in 0.1 M phosphate buffer before shipment. Upon receipt, the nerves were postfixed in 1% aqueous osmium tetroxide for 3 to 4 h before dehydration in a graded alcohol series and propylene oxide. After infiltration with a 1:1 mixture of propylene oxide and araldite resin for 4 h, nerves were placed in 100% araldite resin overnight before embedding in fresh araldite resin. Thick sections (1 μm) were cut with glass knives and either stained with toluidine blue prior to light microscopic examination, or stained with paraphenylenediamine prior to morphometry.

#### *2.4. Axonal Size Frequency Distributions and G-Ratios*

Axonal size-frequency distributions of myelinated fibers were performed on transverse sections of selected peroneal nerve biopsies determined to be adequately fixed, free from artifact, and with an intact perineurium. Images were obtained from a single section of each nerve biopsy using the Photoshop image analysis system. Profiles containing paranodal regions or Schmidt–Lanterman clefts were not included. Using a ×60 objective, the final magnification of the digitized image was equivalent to 1 pixel = 0.091 μm. Myelinated fibers were individually identified and selected prior to being sorted with an automated process into bins based on axonal area. G-ratios were calculated as the ratio between the diameter of the axon itself and the outer diameter of the myelinated fiber.

#### *2.5. Single Nucleotide Polymorphism Array Genotyping and Imputation*

The discovery cohort (426 Leonbergers and 91 Saint Bernards) was genotyped by using either the Axiom Canine Set A or HD arrays (Thermo Fisher Scientific, Waltham, MA, USA) or the Illumina CanineHD BeadChip array (Illumina, San Diego, CA, USA). Samples genotyped on lower single nucleotide polymorphism (SNP) density arrays were imputed with Beagle 4.1 [17,18], using a diverse reference dataset containing 526,045 variants in 49 wolves and 2871 dogs (including 65 Leonbergers and 23 Saint Bernards; Supplementary Table S1). The data were filtered to include only biallelic SNPs with a minor allele frequency ≥0.02, a per-SNP genotyping rate ≥95%, and a per-individual genotyping rate ≥95%. The reference dataset was phased on a per-chromosome basis, using Beagle 4.1 with default parameters of 10 iterations and an effective population size of 200. Next, a target dataset containing approximately 174,000 variants in 402 Leonbergers and 66 Saint Bernards (Supplementary Table S1) was filtered to include only biallelic SNPs with a minor allele frequency ≥0.02, and then checked for concordance with the filtered and phased reference dataset, using the Beagle 4 utility conform-gt [19]. Conforming sites of the target dataset were imputed to the reference dataset on a per-chromosome basis, using Beagle 4.1 with the following settings: window size 50 kb, overlap 3 kb, effective population size 200, and 10 iterations. The per-chromosome imputed data were concatenated

and sorted, using VCFtools 0.1.13 [20]; variants with a Beagle 4.1 dosage R-squared (DR2) ≥0.7 were retained for downstream analysis. In total, we used imputed SNP data for 468 Leonberger and Saint Bernard dogs in this study.

To evaluate haplotypes across breeds, approximately 126,000 SNPs common across genotyping platforms were extracted from non-imputed SNP genotype data for 12,931 canids, which were either generated during this study or publicly available (Supplementary Table S1) and phased with Beagle 4.1 as described above.

#### *2.6. Genome-Wide Association Study and Fine-Mapping*

The discovery cohort from the two breeds combined contained 517 dogs (144 LPPN-cases and 373 controls). Quality control filtering steps of the imputed SNP array genotyping data were carried out by using PLINK v1.9 [21]. The dataset was pruned for low minor allele frequency (0.05) and failure to meet Hardy–Weinberg equilibrium (0.0001) and consisted of 289,553 markers. An across-breed genome-wide association study (GWAS) was performed with GEMMA v0.98 [22], using a linear mixed model including an estimated kinship matrix from centered genotypes to correct for the genomic inflation. The significance threshold was estimated by Bonferroni correction. Manhattan and Q–Q plots of the corrected *p*-values were generated in R environment v3.6.0 [23], using the qqman package [24]. Haplotypes around the significantly associated locus obtained from GWAS were constructed by using Beagle 4.1 for all canids with available SNP genotype data (*n* = 13,399) (Supplementary Table S1). All genome positions refer to the CanFam3.1 reference assembly.

#### *2.7. Whole-Genome Sequencing*

Whole-genome sequence (WGS) data of 716 publicly available dogs of 131 different breeds, and nine wolves [25] (Supplementary Table S2) were studied in order to identify the causative variant in the disease-associated region obtained by the GWAS. This set included 34 Leonbergers and two Saint Bernards diagnosed with a form of LPPN unexplained by the previously known variants in *ARHGEF10* [8], *GJA9* [9], *RAB3GAP1* [10–12] and *RAPGEF6* [13], as well as seven Leonbergers and one Saint Bernard used as controls. The sequence data analysis and calling of single nucleotide variants and small indels (SNVs), including the prediction of functional effects, were described previously [25]. The Integrative genomics viewer (IGV) software 2.8.2 [26] was used for visual inspection and screening for structural variants in the region of interest in the affected dogs' WGS.

#### *2.8. Targeted Genotyping*

Polymerase chain reaction (PCR) and Sanger sequencing were used to validate and genotype the variant identified from WGS. PCR products from genomic DNA were amplified by using AmpliTaqGold360 MasterMix (Thermo Fisher Scientific), and the purified PCR amplicons were directly sequenced on an ABI3730 capillary sequencer (Thermo Fisher Scientific). The *CNTNAP1* missense variant (XM\_548083.6:c.2810G>A) was genotyped, using the following primers: TCCCTTGCCCTCCCTATATC (forward) and AGTCCTAATGCCCTCTGCTG (reverse). The sequence data were analyzed by using Sequencher 5.1 software (GeneCodes, Ann Arbor, MI, USA).

#### *2.9. Protein Predictions*

The MutPred2 [27], PROVEAN [28] and PON-P2 [29] in silico prediction tools were used to predict biological consequences of the discovered variant on the encoded protein. All references to the canine *CNTNAP1* gene correspond to the accessions NC\_006591.3 (NCBI accession), XM\_548083.6 (mRNA), and XP\_548083.3 (protein). The Genome Aggregation Database (gnomAD) [30] was searched for the corresponding variant in the human *CNTNAP1* gene (NP\_003623.1).

#### *2.10. Availability of Data and Material*

The WGS are freely available at the European Nucleotide Archive (ENA). All accession numbers of the used genomes are available in the Supplementary Table S2. The sources of SNP array genotyping data published before are detailed in the Supplementary Table S1, and the dataset generated for this study is available from the corresponding author on reasonable request. All genome positions are reported with respect to the dog reference genome assembly CanFam3.1 and NCBI annotation release 105.

#### **3. Results**

#### *3.1. Phenotype*

The herein studied affected dogs showed generic signs of LPPN (Supplementary Table S1) with the key feature across breeds being breathing difficulty, often described as noisy or raspy breathing (Supplementary Video S1). Due to LP, 247 dogs (121 Leonbergers, 10 Saint Bernards, 114 Labrador retrievers, and 2 mixed-breed dogs) underwent an arytenoid lateralization surgery, 25 of which (ten Leonbergers, seven Saint Bernards, six Labrador retrievers, and two mixed-breed dogs) tested homozygous for the studied *CNTNAP1* variant (Supplementary Table S1). Additional clinical signs, which were noted variably among the dogs, included difficulty swallowing, changes in barking frequency and quality, high-stepping and uncoordinated gait, stumbling and tripping, exercise intolerance, and limb muscle atrophy.

#### *3.2. Neuropathological and Morphometric Findings*

Peroneal nerve biopsies were evaluated from three archived normal adult Labrador retriever (8–10 years) samples (representative image in Figure 1a), and three LPPN-affected dog samples: a nine-year-old male Labrador retriever (Figure 1b), a three-year-old male Leonberger (Figure 1c), and a three-year-old male Saint Bernard (Figure 1d), all of which tested homozygous for the studied *CNTNAP1* variant. Compared to control nerve, pathological changes were similar among affected dogs of all three breeds and included a subjective decrease in the number of myelinated nerve fibers compared to control nerve (Figure 1b–d) with scattered inappropriately thin myelin sheaths for the axon diameter (best demonstrated in Figure 1b,c). The inappropriately thin myelinated fibers were not found in the nerves of control dogs. Myelin splitting and ballooning, onion-bulb formations, and axonal degeneration were not observed in any of the biopsies.

A histogram of axonal size-frequency distribution of the relative percentage of small (<5 μm) and large (>5 μm) myelinated nerve fibers is shown for the three control Labrador retrievers, the LPPN-affected Leonberger, and the LPPN-affected Labrador retriever (Supplementary Figure S1) described above. As only a partial nerve fascicle was available for the LPPN-affected Saint Bernard, those data were not included. The large and small nerve fibers are determined by the axon diameters and this does not refer to the thickness of the myelin sheath. Compared to the average values for the control Labrador retrievers, the affected Labrador retriever showed an increased population of small caliber nerve fibers and a decreased population of large caliber nerve fibers. In contrast, the affected Leonberger showed a decreased population of small fibers and an increased population of larger fibers. Calculated G-ratio, a quantitative measure of myelin thickness, was 0.586 ± 0.031 (range 0.552–0.609) for the control Labrador retrievers, 0.543 for the affected Leonberger, and 0.575 for the affected Labrador retriever.

#### *3.3. Genome-Wide Association Study and Fine-Mapping*

The across-breed GWAS using the discovery cohort of *ARHGEF10-* and *GJA9-* negative Leonbergers and Saint Bernards (144 cases vs. 373 controls) revealed a single genome-wide significantly associated region for LPPN (Figure 2a). The 15 best-associated markers were used to define a 4.6 Mb region of interest between 19.1 and 23.7 Mb on chromosome 9 (Supplementary Table S3). Fine-mapping of this region, using the available haplotypes from non-imputed SNP array genotyping data, included 87 markers centered on the best associated SNP (chr9:20,271,681). This revealed one homozygous haplotype present most frequently in LPPN-affected dogs (*n* = 21) of both breeds and not present in homozygosity in any of the controls (Supplementary Table S4). Therefore, the disease-associated region was narrowed to ~0.98 Mb (bp position 19,393,936 to 20,371,611), by a combination of sharing in the 21 homozygous cases from the discovery cohort (14 Leonbergers and 7 Saint Bernards), coupled with recombination events. Based on this analysis, we hypothesized that the causative variant explaining the GWAS hit was localized on this specific haplotype occurring in both breeds from the discovery cohort. Subsequent haplotype analysis of all 13,399 canids with available SNP genotype data provided evidence for the presence of this haplotype in total of 25 dog breeds (Supplementary Table S4).

#### *3.4. Identification of the Candidate Causative Variant*

In total, 38 protein coding genes and five lncRNAs were annotated within the disease-associated ~0.98 Mb critical interval on chromosome 9 (Figure 2b). Visual inspection of this region in the WGS of five LPPN-affected dogs homozygous for the associated haplotype (four Leonbergers and one Saint Bernard) revealed no evidence for the presence of structural variants. Filtering variants for homozygous alternative genotypes shared in these five dogs within the critical interval yielded 872 intronic or intergenic, 12 synonymous, and 18 protein-changing variants (Supplementary Table S5). In addition, WGS data were available for three Leonbergers out of the 55 dogs from the discovery cohort (Supplementary Table S1) carrying a single copy of the identified disease-associated haplotype (Supplementary Table S4). Filtering for heterozygous variants in complete linkage disequilibrium with this haplotype in three dogs reduced the number of putative variants to 93 intronic or intergenic, one synonymous, and one protein-changing variant (Supplementary Table S5).

**Figure 1.** Paraphenylenediamine-stained resin sections from the peroneal nerves of four dogs. (**a**) An adult normal control Labrador retriever, (**b**) a nine-year-old LPPN-affected Labrador retriever, (**c**) a three-year-old LPPN-affected Leonberger, and (**d**) a three-year-old LPPN-affected Saint Bernard. All three LPPN-affected dogs were homozygous for the *CNTNAP1* variant. Arrows in (**b**) and (**c**) point to nerve fibers that are inappropriately thin for the axon diameters. Bar in lower right image indicates 25 μm and is valid for all images.

**Figure 2.** Identification of a new LPPN-associated locus and variant in Leonbergers and Saint Bernards. (**a**) Manhattan plot for the two-breed genome-wide association study (GWAS) using 144 LPPN-affected dogs and 373 normal control dogs indicates a signal with multiple associated single nucleotide polymorphisms (SNPs) on chromosome 9. The -log *p*-values for each SNP are plotted on the *y*-axis versus each canine chromosome on the *x*-axis. The red line represents the Bonferroni corrected significance threshold (−log(*p*-value) = 6.76). Inset: Corrected Q–Q plot confirms that the observed *p*-values of the best-associated markers have stronger association with the trait than expected by chance (null hypothesis, red line). (**b**) Gene content in the ~0.98 Mb region of interest. Blue arrows show the best-associated markers from GWAS. Green bars represent the different genes and violet bars represent lncRNAs. (**c**) Schematic representation of the *CNTNAP1* gene showing the variant (XM\_548083.6:c.2810G>A) location in exon 18. (**d**) Schematic representation of the CNTNAP1 protein with its domains: coagulation factor 5/8 C-terminal domain (FA58C), laminin G domains (LamG), calcium-binding EGF-like domain (EGF), fibrinogen-related domain (FReD), and putative band 4.1 homologues' binding motif (4.1m). The CNTNAP1 amino acid substitution (XP\_548083.3:p.Gly937Glu) position is shown together with the amino acid multiple alignment indicating the residue is highly conserved across species.

*CNTNAP1* represents a functional candidate gene due to its involvement in human congenital hypomyelination, where vocal cord palsy is a common clinical finding [31], so we pursued the lone remaining missense variant in the *CNTNAP1* gene (CFA9:g.20298261C>T; c.2810G>A; p.Gly937Glu) further. This *CNTNAP1* variant was predicted to be deleterious by several prediction tools (MutPred2

score: 0.884, PROVEAN score: −7.667, PON-P2 probability for pathogenicity: 0.848). It is located in exon 18 of the *CNTNAP1* gene (Figure 2c) and affects a highly conserved amino acid residue at the end of the third laminin G domain of the CNTNAP1 protein (Figure 2d). Two missense variants (rs905697967:p.Gly938Arg and rs763033339:p.Gly938Glu) in the human *CNTNAP1* coding region at the corresponding position were found in the gnomAD [30]. Both variants were reported with allele frequency 7.95 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and no homozygous individuals were detected [30]. The canine missense variant in *CNTNAP1* was present in 5 of the 688 control canid WGS with a frequency of 0.004 (Supplementary Table S5), including one homozygous (English bulldog) and four heterozygous dogs (golden retriever, Labrador retriever, English bulldog, and Kerry blue terrier).

#### *3.5. The CNTNAP1 Variant Occurs in Several Breeds*

Available SNP array genotype data of 13,337 dogs and 62 wild canids (Supplementary Table S1) were inspected for the identified *CNTNAP1*-associated haplotype (Supplementary Table S4). Out of this group, targeted genotyping by PCR and Sanger sequencing was performed in 2469 canids and demonstrated perfect concordance between the *CNTNAP1*-associated haplotype and the *CNTNAP1*:c.2810G>A genotype (Supplementary Table S4). In addition, 2362 dogs without SNP array data were directly genotyped for the variant (Supplementary Table S1); this included 557 dogs from the validation cohort and 1805 from the population cohort.

In total, the variant was found in 25 different breeds and no wild canids. Homozygotes for the missense allele in breeds other than the Leonberger, Saint Bernard, and Labrador retriever were identified in 46 English bulldogs, six Irish terriers, two boxers, one bullmastiff, one Peruvian hairless dog, one Yorkshire terrier, and one golden retriever (Supplementary Table S1), all with unknown precise health history. Additionally, two LPPN-affected mixed-breed dogs enrolled in the Michigan State University's GOLPP study [6] were also homozygous for the missense allele (Supplementary Table S1).

Analysis of the validation cohort that included the three breeds with available health information (Leonberger, Saint Bernard, and Labrador retriever) demonstrated that the *CNTNAP1* variant is not present in a homozygous state in any dog apparently non-affected with LPPN (Table 2). The 18 homozygous LPPN-affected Leonbergers represent 4.1% of all as yet unexplained cases with any age of onset that were not carrying the previously identified disease-causing variants in *ARHGEF10* and *GJA9*. For the Saint Bernard, 10 out of 24 (41.6%) diagnosed dogs were homozygous mutant, whereas only 4.7% of GOLPP-affected Labrador retrievers carried two copies of the *CNTNAP1* variant. The homozygous A/A *CNTNAP1* genotype also occurred rarely in single dogs of the population controls from each of the three studied breeds (Table 2). Altogether, the mutant allele frequency was estimated as 6.6% in the studied Leonbergers (*n* = 2738 dogs), 13.9% in Saint Bernards (*n* = 305 dogs), and 5.2% in Labrador retrievers (*n* = 1524 dogs). Among the 22 other breeds segregating for this variant, the allele frequency was highest in the English bulldogs (*n* = 193 dogs) and Irish terriers (*n* = 184 dogs), estimated at 46.6% and 17.1%, respectively (Supplementary Table S1).

Mean age of onset in the LPPN-affected dogs was investigated for a subset of 770 dogs with detailed health information from three breeds (596 Leonbergers, 28 Saint Bernards, and 146 Labrador retrievers) and showed a marked difference between the cases depending on the different underlying genetic causes (Figure 3). The average age of onset of clinical signs in the limited number of dogs homozygous for the herein described *CNTNAP1* variant was 3.4, 2.1, and 7.5 years in Leonbergers, Saint Bernards, and Labrador retrievers, respectively. In comparison, the age of onset of clinical signs in the previously characterized *ARHGEF10*-associated polyneuropathy [8] was seen in Leonbergers and Saint Bernards with average ages of 2.2 and 1.6 years, respectively. Additionally, affected Leonbergers with the *GJA9* frameshift variant [9] had average age of onset of their clinical signs of 6.2 years. Interestingly, the LPPN-affected Labrador retrievers that do not carry the herein identified *CNTNAP1* variant and come from the GOLPP study [6] showed a higher average age of onset of 11.5 years (Figure 3). The difference in age of disease onset between the dogs with the identified *CNTNAP1* variant and the cases without known disease-causing mutation was statistically significant in all

three breeds (Leonberger *p*-value = 0.000001002, Saint Bernard *p*-value = 0.01681, Labrador retriever *p*-value = 0.002662). The difference between dogs with the *CNTNAP1* variant and the *ARHGEF10* variant was significant in Leonbergers (*p*-value = 0.002538) but not significant in Saint Bernards (*p*-value = 0.3095). The difference between Leonbergers with the *CNTNAP1* variant and the *GJA9* variant was statistically significant (*p*-value = 0.00000001797).

**Table 2.** Segregation of the *CNTNAP1*:c.2810G>A genotypes with laryngeal paralysis and polyneuropathy (LPPN) in three breeds with available health information.


<sup>1</sup> Includes the three herein described histopathologically confirmed cases. <sup>2</sup> LPPN-affected dogs that tested negative for *ARHGEF10* [8] and *GJA9* [9] mutations. <sup>3</sup> LPPN-affected dogs homozygous for the previously described polyneuropathy-associated variant in *ARHGEF10* (LPN1) [8], and homozygous or heterozygous for the variant in *GJA9* (LPN2) [9]. <sup>4</sup> Includes cases homozygous for the previously described leukoencephalomyelopathy-associated variant in *NAPEPLD* [14].

**Figure 3.** Age of onset of clinical signs forlaryngeal paralysis and polyneuropathy (LPPN) differs depending on underlying genetic variants and across the three breeds. Comparison of the age of onset of clinical signs in the LPPN-affected dogs (*n*=770), which were genotyped homozygous for the polyneuropathy-associated variants in *CNTNAP1*:c.2810G>A, or *ARHGEF10*:c.1955\_1958+6delCACGGTGAGC [8], and homozygous or heterozygous for the variant in *GJA9*:c.1107\_1108delAG [9], as well as the age of onset in the yet unexplained cases (other), is shown for Leonbergers (gold bars; *n* = 596), Saint Bernards (pink bars; *n* = 28), and Labrador retrievers (yellow bars; *n* = 146). Note that the *ARHGEF10* variant is only present in Leonbergers and Saint Bernards, and the *GJA9* variant only in Leonbergers.

#### **4. Discussion**

This study has revealed strong evidence for a new potentially pathogenic variant associated with laryngeal paralysis and polyneuropathy (LPPN) initially observed in two closely related giant dog breeds, the Leonberger and Saint Bernard. Interestingly, the variant also explains some cases of GOLPP in Labrador retrievers and segregates at different frequencies in 22 other unrelated dog breeds, including English bulldogs and Irish terriers, suggesting that the derived allele predates modern breed formation. Apparently, the variant was not purged by either selection or drift.

The affected *contactin-associated protein 1* (*CNTNAP1*) gene has been previously implicated in human autosomal recessive neurological diseases with a broad spectrum of clinical phenotypes and neonatal and childhood onsets: congenital hypomyelinating neuropathy type 3 (OMIM 618186), lethal congenital contracture syndrome 7 (OMIM 616286), and childhood-onset Charcot–Marie–Tooth disease [32]. *CNTNAP1* is essential in the formation of paranodal axoglial junctions in myelinated axons and is also involved in regulating neural progenitor cells and the development of the cerebral cortex [33]. Pathological variants in the *CNTNAP1* gene may lead to defective or absent proteins critical to development of central or peripheral nervous systems. Even though, based on the current gnomAD [30] database, the corresponding glycine to glutamic acid exchange occurs very rarely in humans, so far, there is no evidence reported for disease association. Human patients with other *CNTNAP1* homozygous frameshift or nonsense variants show a more severe disorder with early-onset neurological disease, including severe respiratory compromise and early lethality, while those carrying missense variants can survive beyond infancy [34]. This suggests that the missense alleles affecting the myelination and development of paranodal junctions may be hypomorphic and have some residual function. Although the precise role of the protein domains in ligand binding is not fully understood, several missense variants were predicted to impact the domain structure and protein folding [34].

Pathological changes in semi-thin transverse resin sections of the peroneal nerves of three affected dogs had similarities to those in published human cases, including reduced myelinated nerves and inappropriately thin myelin sheaths for the axon diameters [35]; however, in the affected dogs, thinly myelinated fibers were scattered and fewer in number. This may reflect the severity of disease in neonatal onset human cases and milder disease with an adult onset in dogs. The increased population of small fibers in the affected Labrador retriever, as compared to the affected Leonberger, may reflect attempts at regeneration in the Labrador retriever or genetic differences in other modifying genes. There are several limitations to the pathological studies in these dogs: The number of affected dogs with available peripheral nerve biopsies was small; detailed study of the paranodal areas of the peroneal nerves was limited by the retrospective nature of the study, the use of archived nerve specimens obtained many years prior to the identification of the variant, and the necessity of preparation of the nerves at the time of original processing for teased fibers and for longitudinal evaluation; and the standard processing for diagnostic specimens in the laboratory of one of the authors (GDS) is in transverse section. Future in-depth prospective studies of the peripheral nerves, including laryngeal nerves in more cases of confirmed LPPN with the *CNTNAP1* gene variant in each breed, are necessary to fully evaluate the observed pathological changes.

The neurological diseases identified in humans associated with variants in *CNTNAP1* support our recent speculation, based on the enrichment of this allele in Leonbergers [7], that the herein-described missense variant predicted in silico to be deleterious represents a promising candidate causative mutation for inherited neurological disorders in dogs. The striking genetic association data implicate that this mutation affects the function of the encoded protein, although we have not studied this further. Homozygosity for the missense variant in the *CNTNAP1* gene is significantly associated with the development of LPPN in large and giant-sized dogs, indicating recessive inheritance in all three studied breeds (Leonberger, Saint Bernard, and Labrador retriever). As yet we do not have convincing evidence for causality in smaller dog breeds segregating for this variant, such as the English bulldog or Irish terrier, and further study with reliably phenotyped populations is needed. However, we hypothesize that the apparently higher allele frequency in English bulldogs may be a result of

underdiagnosed LP due to breathing difficulties related to brachycephalic airway syndrome, including laryngeal collapse [36] obscuring a neurodegenerative LP. The observed higher allele frequency in Irish terriers, although without available health information, suggests that the association between the variant and LPPN phenotype is breed-specific and may not be pathogenic in some breeds although this needs to be evaluated. We also hypothesize that the observed later age of onset in the Labrador retriever group, compared to the Leonbergers and Saint Bernards, might be either due to the different genetic breed background and/or their smaller stature and correspondingly shorter laryngeal nerve length; the latter was previously suggested by correlation between growth (specifically height) and laryngeal neuropathy in horses [37].

#### **5. Conclusions**

In conclusion, we identified a potentially causative genetic variant in *CNTNAP1* associated with autosomal recessive younger-onset LPPN in large and giant dogs, specifically Leonbergers, Saint Bernards, and Labrador retrievers. Our results represent the first large animal model for a *CNTNAP1*-related neurodegenerative disease. The developed genetic test enables veterinary diagnostics and selective breeding against this deleterious variant across breeds to reduce the occurrence of LPPN. Therefore, selecting based on this additional disease-associated variant, which we have designated LPPN3, will enable dog breeders to make even greater strides in controlling the propagation of this devastating disorder and maintaining the health of Leonberger, Saint Bernard, and Labrador retriever populations. However, the fact that not all LPPN cases from the three intensively studied breeds carried the described variant, together with the broad range in age of onset of the clinical signs for the as yet unexplained cases, indicates that still unknown genetic heterogeneity of different forms of canine LPPN need to be studied in future.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/12/1426/s1. Figure S1: Histogram showing the axonal size-frequency distribution of myelinated nerve fibers in peroneal nerve sections from control Labrador retrievers (average of three archived adult control dog samples), an LPPN-affected Labrador retriever, and an LPPN-affected Leonberger. Both affected dogs were homozygous for the *CNTNAP1* variant. Axons are distributed into bins determined by axonal diameter based on perimeter. Table S1: Sample designations and detailed information of all dogs used in the study. Table S2: Sample designations and affection status of all whole-genome sequences used for filtering variants. Table S3: The 10,000 most significant markers sorted by *p*-value obtained from the across-breed genome-wide association study. Table S4: Haplotype diversity for the LPPN-associated genome region on chromosome 9. Table S5: Variants in the region of interest obtained from GWAS on chromosome 9. Video S1: Video illustrating the clinical phenotype of two LPPN-affected dogs: a Leonberger and a Saint Bernard.

**Author Contributions:** Conceptualization, J.R.M. and C.D.; data curation, B.J.S. and V.J.; formal analysis, A.L. and K.M.M.; funding acquisition, J.R.M. and C.D.; investigation, A.L., K.M.M., G.D.S., and J.P.S.; methodology, A.L. and K.M.M.; resources, P.J.J.M., P.A.J.L., P.A.W., S.M.P.-J., K.J.E., G.S.J., and L.H.; software, S.G.F. and V.J.; supervision, J.R.M. and C.D.; visualization, A.L., K.M.M., and G.D.S.; writing—original draft, A.L., K.M.M., G.D.S., J.R.M., and C.D.; writing—review and editing, A.L., K.M.M., S.G.F., G.D.S., K.J.E., J.R.M., and C.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The University of Minnesota received a gift from the Leonberger Health Foundation to support research in this breed. Kari J. Ekenstedt (K.J.E.) was supported by the Office of the Director, National Institutes of Health (NIH), under award number K01-OD027051. Other authors have no external funding to report.

**Acknowledgments:** The authors are grateful to the breeders and owners of all dogs who provided samples and shared valuable information about their dogs. We thank Nathalie Besuchet Schmutz for expert technical assistance. The Next Generation Sequencing Platform and the Interfaculty Bioinformatics Unit of the University of Bern are acknowledged for performing the WGS and providing high-performance computational infrastructure.

**Conflicts of Interest:** Both the University of Minnesota and the University of Bern are offering genotyping tests for polyneuropathy- and leukoencephalomyelopathy-associated variants in their respective laboratories, and proceeds from these tests go toward ongoing canine genetic research. S.G.F is guest editor of this special issue of *Genes*, but has not in any way been involved in or interacted with the journal's review process or editorial decision-making. The authors declare that they have no other competing interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Genes* Editorial Office E-mail: genes@mdpi.com www.mdpi.com/journal/genes

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18