Animals
Genomic DNA was isolated from EDTA-blood, hair root or tissue samples. We sampled 42 standard Munchkin and 15 non-standard Munchkin cats from several subtypes (Additional file 12). Standard Munchkin cats included the subtypes Genetta (n = 10), Napoleon (n = 2), Dwelf, (n = 1), Bambino (n = 10), Lambkin (n = 1), longhaired (n = 11) and shorthaired (n = 7). Non-standard Munchkin cat samples were available for the subtypes Genetta (n = 2), Bambino (n = 6), longhaired (n = 1) and shorthaired (n = 6). Cats of the subtype Genetta belonged to a family segregating into long-legged and short-legged types (Additional file 3). The sire and three of his kittens (two standard and one non-standard Munchkin cats) underwent a detailed clinical examination. The parents and two full sibling kittens from this segregating litter were chosen for WGS. For controls, we sampled 203 non-Munchkin cats representing British shorthair, Bengal, Ragdoll, Exotic Shorthair, Don Sphinx, Birman, Chartreux, Maine Coon, Norwegian Forest, Persian, Russian Blue, Scottish Fold, Siamese, Siberian, Turkish Angora, Oriental Shorthair, and domestic shorthair.
Clinical examination
For the clinical investigation, members of the segregating Munchkin cat family of the subtype Genetta were available. A cat breeder provided a 4-year old standard Munchkin tomcat and three of his progeny including two male standard and a female non-standard Munchkin kitten for a detailed examination. The 4-year old standard Munchkin tomcat underwent a CT scan in sternal recumbency using a multislice helical CT scanner (Brilliance 64-CT, Philips Medical Systems, Best, The Netherlands). A slice thickness for the head of 0.67 mm (0.026 in) and settings of 120 kV/ 250 mAs and for the body a slice thickness of 0.67 mm (0.026 in) and settings of 140 kV/ 250 mAs were set for analysis. A 1024 × 1024 matrix was used.
In addition, absolute total length, diaphyseal length and the diaphyseal diameter of humerus, radius, ulna, femur, tibia, metacarpalia and metatarsalia bones were measured according to Jones and Jolly [31]. The relative length of individual bones was determined through the ratio of length of the respective individual bone to the length measured from the proximal end of the humerus to the distal end of metarcarpus III.
Haplotype analysis
To refine the critical genome-wide associated region on FCA B1, Kompetitive Allele Specific PCR (KASP) assay was used (LGC Genomics, Middlesex, UK). We selected 14 SNPs (PRJEB30080) of which four SNPs were located distally to the critical region and ten SNPs within the genome-wide associated region at a distance among each other of about 430 to 7600 kb on FCA B1. These 14 SNPs were genotyped in 9 standard Munchkin cats, one non-standard Munchkin cat and 79 controls. These animals were unrelated. KASP genotyping was performed according to standard protocols (LGC) using an ABI7300 real-time system for 96 well plates (Additional file 13). A case-control analysis was carried out using SAS/Genetics, version 9.4 (Statistical Analysis System, Cary, NC, USA). Haplotype analysis was done with the procedure HAPLOTYPE in SAS.
Whole genome sequencing
Animals for WGS were from the Munchkin cat family of the subtype Genetta including the standard Munchkin parents (sire and dam) and two members of their litter with a standard Munchkin cat male and a non-standard Munchkin cat female (Additional file 3). Libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, USA) and run on an Illumina NextSeq500 in a 2 × 150 bp paired-end mode. Quality control was performed using fastqc 0.11.5 [32] and reads were trimmed using PRINSEQ (V 0.20.4) [33].
WGS data was mapped to the cat reference genome Felis catus 8.0 (Ensembl) using BWA 0.7.13 [34]. Sorting, duplicate marking and indexing was performed using Picard tools (http://broadinstitute.github.io/picard/, version 2.9.4) and SAMtools 1.3.1 [35]. For final variant calling GATK version 3.7 [36] with the Base Quality Score Recalibrator (BQSR), Haplotype Caller [37] and Variant Recalibrator were used. Variant calling data were compared to WGS data from 16 different individuals including Bengal (SAMN05980341, SAMN05980358), Siamese (SAMN05980342), Peterbald (SAMN05980359, SAMN14103114), Oriental shorthair (SAMN14103123, SAMN14103124, SAMN14103127), and domestic shorthair (SAMN14103116, SAMN14103119, SAMN14103120, SAMN14103121, SAMN05980374, SAMN05980340, SAMN05980323, SAMN05980325).
A read depth of 8–999 and quality score values > 20 were applied for variant detection. Those variants, which were heterozygous in the three standard Munchkin cats, homozygous wild type in the non-standard Munchkin cat and homozygous wild type in all controls and concordant to an autosomal monogenic dominant inheritance, were filtered using SAS (version 9.4, Statistical Analysis System, Cary, NC, USA) for further investigation. Only variants with high or moderate effects according to SNPEff version 4.3 t (2017-11-24, SNPEff database Felis catus8.0) [38] were considered for further analysis. A candidate gene list from previous studies about chondrodysplasia was used to find candidate genes in the filtered variant list [24].
Structural variant detection
For structural variant detection, LUMPY software, version 0.2.13 [39], integrating multiple structural variation signals jointly across multiple samples, was run for Bam files of all four Munchkin cats and 16 controls. The generated VCF file was filtered for structural variants on FCA B1, heterozygous in the standard Munchkin cats and homozygous wild type in the non-standard Munchkin cat and the 16 controls.
Mutation analysis
To test whether the 108 bp insert is transcribed or not transcribed into mRNA, we isolated RNA from hair roots stabilized in RNAlater reagent (Qiagen, Hilden, Germany) of two standard Munchkin cats of the subtype Genetta and two domestic shorthair cats. All samples were transcribed into cDNA according to standard protocols (Maxima First Strand cDNA Synthesis Kit, Thermo Fisher). Primer pairs were designed using Primer3 tool (version 0.4.0, http://bioinfo.ut.ee/primer3-0.4.0/). One primer pair was designed to produce a PCR-product in mutant and wild type allele (PCR-type 1), one primer pair to produce a PCR-product only if the mutant allele is present and the insertion is transcribed (PCR-type 2), and one control primer pair to produce a PCR-product in the presence of wild type allele (PCR-type 3, Additional file 14). The first primer pair (PCR-type 1), located across the deletion from exon 10 (ENSFCAT00000009602.6) or exon 9 (ENSFCAT00000055794.2) to the 3′-UTR within the non-deleted region, yields an expected amplicon size of 510 bp for the mutant allele if the insertion is transcribed and 402 bp for the mutant allele if the insertion is not transcribed, as well as 2514 bp for the wild type allele (Additional file 7). The second primer pair (PCR-type 2) from the proximal region of exon 10 (ENSFCAT00000009602.6) or exon 9 (ENSFCAT00000055794.2), to a region within the 108 bp insert produces an expected amplicon of 110 bp only in the mutant allele when the insertion is transcribed. This PCR-product is not present in the wild type or mutant allele when the insertion is not transcribed. The third primer pair (PCR-type 3) was designed to bridge exon 10 and 11 (ENSFCAT00000009602.6) or exon 9 and 10 (ENSFCAT00000055794.2) with an expected amplicon size of 185 bp for the wild type allele in standard Munchkin cats and control cats with wild type alleles on both chromosomes. Furthermore, we designed five primer pairs (PCR-type 4–7) for cDNA sequencing the region spanning exon 10 (ENSFCAT00000009602.6) or exon 9 (ENSFCAT00000055794.2) and those parts of the 3’UTR present in the mutant and wild type allele (Additional file 15). PCR included 0.25 μM of the respective primers, 20 ng cDNA, 12.5 μl UCP HiFidelity MasterMix (Qiagen, Hilden, Germany) and was prepared as recommended by Qiagen, Hilden, Germany. Sanger sequencing was done for the standard Munchkin dam from the Munchkin cat family of the subtype Genetta (Additional file 3) for PCR-type 1 and in one domestic shorthair cat for PCR-type 4–7.
In addition, open reading frames were predicted using ORF Finder (http://www.ncbi.nlm.nih.gov/projects/gorf/). Protein domains predicted by InterProScan (http://www.ebi.ac.uk/interpro/search/sequence-search) for PANTHER [40] and PIRSF [41] were obtained by Ensembl database (protein summary) [42]. For comparative species alignments Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) was used [43]. Primer blast was done using NCBI nucleotide blast (https://blast.ncbi.nlm.nih.gov/Blast.cgi? PROGRAM = blastn&PAGE_TYPE = BlastSearch&LINK_LOC = blasthome). The sequence of the 108 bp insertion was compared to the reference genomic sequences Felis catus 8.0 and Felis catus 9.0. The sequence similarity, highest alignment score and the expected value were calculated.
Validation
A structural variant composed of a 3303 bp deletion (ss5015497294) and a 108 bp insertion (ss5015497295) within the critical region on FCA B1 was Sanger sequenced and visualized using a duplex PCR in individuals including standard Munchkin, non-standard Munchkin and control cats. Duplex PCR was performed in 260 individuals including 15 non-standard Munchkin, 42 standard Munchkin and 203 controls from breeds other than Munchkin cats. Sanger sequencing was done in the standard Munchkin sire, dam and the standard Munchkin kitten. A forward primer proximally of the detected structural variant and two different reverse primers located within the deleted region and distal of this deletion were designed using Primer3 tool (version 0.4.0, http://bioinfo.ut.ee/primer3-0.4.0/) (Additional file 16 and Additional file 17). Duplex PCRs were performed in 22-μl reaction volumes containing 2 μl DNA, 1.5 mM deoxyribonucleoside triphosphates, 5 pmol of primers MK_del_R and MK_wt_R, 10 pmol of primer MK_wt_F, 1.5 U of Taq polymerase in the reaction buffer supplied by the manufacturer (MP Biomedicals, Eschwege, Germany) and 4.2 μl enhancer reagent (MP Biomedicals). After a 5 min initial denaturation at 95 °C, 40 cycles of 30 s at 94 °C, 30 s at 59 °C, and 45 s at 72 °C were run on a Thermocycler TProfessional 96 (Biometra, Göttingen, Germany). All samples were visually evaluated on a 1% agarose gel using Gel iX20 Imager (Intas Science Imaging Instruments, Göttingen, Germany).
Break point validation was performed in the three standard Munchkin cats of the subtype Genetta with WGS data using the primers MK_wt_F and MK_del_R accoding to the duplex PCR protocol. Evaluation of Sanger sequences was performed using Sequencher 4.8 software (GeneCodes, Ann Arbor, MI, USA).