Phenotype

On gross examination, the phenotypic appearance of the malformation of the three affected calves (bulldog calves 1–3) was highly similar. The calves were born at term and fully developed, but were dead at birth or died shortly after birth. Bulldog calf 1 was born after 251 days and weighed 16.3 kg. Bulldog calf 2 weighed 27.5 kg after 278 days of gestation and bulldog calf 3 weighed 18 kg after 273 days of gestation. The hair coats of the affected calves were fully developed. Wrinkles of the skin were not obvious. The short and compressed bodies with bulky abdomen, severe bilateral symmetrical shortening of the limbs and severe facial dysmorphisms resembled the “bulldog calf” phenotype with an extremely disproportionate chondrodysplasia (Fig. 1). The faces were characterized by frontal bossing and marked shortening of the muzzles resulting in a protrusion of the tongues (Fig. 2). Bilateral palatoschisis besides the midline with an opening of the hard and soft palate was evident in all three cases. The limbs were severely reduced in length but had normally sized digits.

Computed tomography (CT) was done for all three bulldog calves. The morphology of the skeletons of the three cases was very similar to each other. A severe shortening of the nasal and the incisive bones was evident. The calves had brachygnathia superior with a prolonged lower jaw (prognathism). The splanchnocrania were dysplastic with caudodorsal displacement of the orbitae in opposite to the prominent neurocrania (Fig. 3a, b). These altered ratios provided a spherical appearance of the skull, which is particularly evident in comparison with a non-affected animal (Fig. 3c, d). The chondrodysplastic animals also had a bilateral cleft palate extending through the hard and soft palate. In the comparative examination of the body skeleton of bulldog calf 2 with a phenotypically normal calf, particularly the ribs, the lumbar spine, the pelvis and the skeleton of the limbs had striking malformations (Fig. 4). The ribs were short and distances among them were irregular but the numbers were normal. The opening angle was larger than in the normal calf. The ribs extended more horizontally instead of vertically. All vertebrae were developed and their number did not deviate from normal calves. Especially the vertebrae of the lumbar spine were severely malformed. They were flat (platyspondylia) and the processus transversi of the lumbar spine were completely missing. The vertebrae of the cervical spine were flatter than in the control animal. The pelvis was severely retarded in its development. Its shape was small and stocky. The long bones of the limbs above the phalanges were short and irregularly shaped. The bones of the limbs were poorly mineralized and consisted mainly of cartilage.

Gross and microscopic findings

Post mortem examination demonstrated very similar findings in all three affected calves as exemplified for bulldog calf 2. A diffuse subcutaneous edema (anasarca) was seen in the three present cases. Skeletal abnormalities were characterized by markedly shortened and partially rotated limbs (Fig. 5) and reduced length of the vertebral column. The thorax was dorsoventrally flattened. The head was dome shaped with a protruded mandible and cleft palate. Histology of malformed limbs showed irregular proliferations of densely packed and hypertrophic chondrocytes without distinct growth plates (Fig. 6), indicative of disturbed enchondral ossification. Islands of retained cartilage with thin osteoid layers were found in the epiphysis, metaphysis and diaphysis (Fig. 7) of malformed bones. In addition, articular cartilage was thickened and showed irregular proliferations and disorganization of chondrocytes.

Pedigree analysis

Screening of the veterinary records of the herd of the LuFG Ruthe gave no indication for the occurrence of a similar malformation in the past 40 years. All newborn calves and all aborted calves are under veterinary control and any congenital abnormalities are recorded in a computerized herd management program. The dams, their progeny and any ancestors which were still at LuFG Ruthe were examined but no indications for skeletal abnormalities were found.

The pedigree revealed that all bulldog calves were sired by one polled Holstein bull used for AI. All females inseminated with semen of Energy P got pregnant and no abortion was observed in these animals. In total, 8 calves were sired by this bull and 3/8 calves were bulldog calves (Fig. 8). Inbreeding in the affected calves on ancestors of Energy P could not be seen, whereas two dams of unaffected calves were related with Energy P via male ancestors (sires A and B in the pedigree). Energy P was inbred on a male ancestor (sire B in the pedigree).

DNA analysis

We performed a runs of homozygosity (ROHs) in windows of 15 SNPs for the three bulldog calves using 761,461 SNPs from the Illumina bovine high density beadchip (Illumina, San Diego, CA) using PLINK, version 1.07 (http://pngu.mgh.harvard.edu/purcell/plink/). We were not able to identify consensus homozygosity regions for all three affected calves. The three bulldog calves had 356, 370 and 389 ROHs with an average size of 350, 458 and 344 kb, respectively. We rerun this analysis after having removed all SNPs not being in Hardy-Weinberg-equilibrium (P < 0.00000001). However, we were not able to find consensus ROHs for the three bulldog calves. Subsequently, we employed Illumina whole genome sequencing for one bulldog calf, sperm of its sire Energy P and a normal progeny of Energy P to search for variants shared by a male bulldog calf and the sperm sample of its sire but not shared by a normal male paternal half sibling of the bulldog calf. In addition, we used 19 bovine whole-genome sequences from four different breeds as controls. Reads were mapped to the UMD3.1 bovine reference genome resulting in a mean coverage of 12.4X for the bulldog calf, 13.2X for sperm of Energy P, 8.6X for the normal progeny of Energy P and 10-15X for the 19 controls. Raw data including the Energy P family and the 19 controls contained 29,848,582 variants whereof 2,039,166 were unknown in dbSNP (https://www.ncbi.nlm.nih.gov/snp/).

Filtering for variants with high, moderate or low effects according to SNPEff predictions on the protein structure and variants associated with chondrodysplasia-related genes left only one variant within COL2A1 on BTA5 (UMD3.1). We were not able to identify mutations predicted to have a moderate or high-impact effect on protein structure in any other chondrodysplasia-related gene (Additional file 1). The variant within the bovine gene COL2A1 was a heterozygous missense mutation (g.32476082G > A, c.2986G > A, ss2019324576) within exon 42 based on the transcript variant COL2A1–201 (Ensembl). This mutation causes an exchange of glycine to serine (p.996G > S). Using SIFT [12] and PolyPhen-2 [13] this missense mutation was predicted to be deleterious with a score of 0 (SIFT) or damaging with a score of 1.0 (PolyPhen-2). A search of the NCBI and Ensembl databases for SNPs could not retrieve this missense mutation.

Sanger sequence analysis of the COL2A1 gene for two bulldog calves and an unrelated Holstein Friesian control animal from another herd as shown in Fig. 9 confirmed the ss2019324576 variant in bulldog calves and showed the absence in the control. The semen sample of Energy P revealed the presence of the mutant allele A with a much smaller size peak than the samples from the bulldog calves. The estimated ratio among the area under allele A in comparison to the area of both alleles was approximately 5%. We developed a PCR (mismatch)-restriction fragment length polymorphism (PCR-RFLP, Additional file 2) and validated this PCR-RFLP in the three bulldog calves and the semen sample of Energy P (Additional file 3). The dams of the three bulldog calves, an unrelated control animal (Additional file 3) and all five normal progeny of Energy P were homozygous wildtype (Additional file 4).

Subsequently, we genotyped 105 polled Holstein bulls and 105 normal Holstein calves whereof 48 were descendants from cows of the dairy herd of the LuFG Ruthe and in addition, 66 animals from other cattle breeds using the PCR-RFLP for the ss2019324576 variant (Additional file 2). We were not able to identify this mutated allele in any other individuals (Additional file 4). Therefore, we performed a genome-wide association study (GWAS) for the three affected calves, 96 polled Holstein bulls and 96 normal Holstein calves using the SNPs of the Illumina bovine high density beadchip and the ss2019324576 variant. After quality control and exclusion of beadchip SNPs with a minor allele frequency (MAF) <0.05, 586,979 SNPs were left for the GWAS. The GWAS gave a highly significant peak for the ss2019324576 variant (P-value = 6.436e-30, after Bonferroni correction for multiple testing), whereas all other SNPs were not significant (P-value > 0.05, after Bonferroni correction for multiple testing). These GWAS results support that this COL2A1-associated mutation is very likely a spontaneous germline mutation and chondrodysplasia causing mutations in other genomic and non-coding regions are unlikely.

We calculated COL2A1 containing haplotypes for all 195 animals using beadchip SNPs and the ss2019324576 variant (SAS/Genetics, Statistical Analysis System, Cary, NC). These haplotypes were restricted to 15 SNPs spanning 46,558 bp from 32,441,205 to 32,487,763 bp on BTA5. We searched for the allelic pattern of the haplotype surrounding the ss2019324576 variant and associated with the mutant allele A in affected calves. After having identified this haplotype only shared by the three bulldog calves (A-A-C-G-A-A-A-G-A-G-[A]-G-G-G-G) we searched all normal animals for the proximally and distally of the ss2019324576 variant identical haplotype (A-A-C-G-A-A-A-G-A-G-[G]-G-G-G-G). In 14.5% of the normal animals, this latter haplotype was present. This might indicate that the haplotype with the wildtype allele G of the ss2019324576 variant may be common among Holsteins, particularly in polled Holstein sires and furthermore, the three bulldog calves share the ss2019324576 A-variant surrounding haplotype with some of the normal calves.

Animals

Genomic DNA was obtained from blood and tissue samples of five normal and three affected progeny of the Holstein bull Energy P, their Holstein dams and sperm from Energy P. In addition, 48 EDTA blood samples from Holstein calves born on the LuFG Ruthe, further 57 randomly collected Holstein calves and 105 polled Holstein bulls were available. Furthermore, 66 blood samples from five different cattle breeds and a crossbred were employed for the present study.

Computed tomography

The CT image acquisition was performed in sternal recumbency of the affected calves. CT scans were acquired with a multislice helical CT scanner (Brilliance-TM-CT-16 BigBore, Philips Healthcare, The Netherlands) using conventional settings (90 kV, 100 mAs) and 2 mm slice thickness.

Histopathological examination

Bone and articular tissues were prepared for the affected calves with a diamond saw and formalin fixed for 48 h. Prior to paraffin embedding samples were decalcified in HNO₃ solution (5%) for 24 h. Subsequently, 3 μm-thick sections were cut on a microtome, mounted on glass slides and stained with hematoxylin and eosin.

Genotyping

Genomic DNA was extracted from the EDTA-blood samples using routine procedures. Concentration of extracted DNA was determined using the Nanodrop ND-1000 (Peqlab Biotechnology, Erlangen, Germany). DNA concentration was adjusted to 50–80 ng/μl. We performed a high-density genome-wide association analysis employing the Illumina bovine high density beadchip comprising 777,962 SNPs with a mean spacing of 3.4 kb across the entire bovine genome. Raw data were analysed using the genotyping module version 1.9.4 of the GenomeStudio V2011.1 software (Illumina). All SNPs underwent quality control. SNPs with a call rate < 0.95 were excluded from further analyses. Due to these specifications, 761,461 SNPs were left for runs of homozygosity (ROHs). Calculation of ROHs was done using PLINK, version 1.07 (http://pngu.mgh.harvard.edu/purcell/plink/). After exclusion of SNPs with a minor allele frequency (MAF) <0.05, 586,979 SNPs were left for the GWAS. In addition to this set of beadchip SNPs, the COL2A1:g.32476082G > A SNV was added for the analysis. The positions of all SNPs and the SNV were according to the UMD 3.1 bovine genome assembly.

Whole-genome sequencing

DNA samples from one male bulldog calf, sperm of the bull Energy P and a normal male offspring of Energy P were prepared for sequencing using the NEBNext Ultra II DNA Library Prep Kit for Illumina for high quality libraries (NEB, New England BioLabs, Ipswich, MA) according to standard NEB-protocols. Size selection and indexing was quality controlled on an Agilent 2100 Bioanalyzer system with a High Sensitivity DNA kit (Agilent Technologies, Santa Clara, CA). Finally, indexed libraries were sequenced on the Illumina NextSeq 500 for 300 cycles in 2 × 150 paired-end modes. After quality control of fastq files using fastqc 0.11.3 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) reads were mapped to the UMD3.1 bovine reference genome using BWA 0.7.12 [18], converted into binary format [19] and underwent variant calling using GATK [20]. Further variant effect prediction and classification of effects was done using SNPEff version 4.1 b [21, 22]. In addition, variant calling and bioinformatic analysis included 19 whole-genome sequences of different cattle breeds (Holstein, Galloway, Limousin, Vorderwald cattle) as controls. All fastq data can be obtained through the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra), BioProject PRJNA350593 (BioProject submission: SUB2036163, SRA experiment IDs: SRX2268677, SRX2268678, SRX2268680).

Variant filtering

Variants from whole-genome sequencing analysis were filtered for high, moderate or low effects according to SNPEff predictions using SAS, version 9.4 (Statistical Analysis System, Cary, NC). Private variants exclusively found in the bulldog calf and Energy P, but not in the normal male half-sib of the bulldog calf and further reference cattle data were filtered out and further controlled for dbSNP entries. Predicted effects of these private variants were investigated using the Variant Effect Predictor [22] for SIFT [12] predictions and PolyPhen-2 [13]. All variants heterozygous for the mutant allele in the bulldog calf and Energy P, but homozygous wildtype in the normal male half-sib as well the control animals were evaluated for their potential impact on chondrodysplasia using a candidate gene list [6] and GO terms. All candidate genes with potential protein structure influencing mutations were visually searched for mutations and deletions >10–15 bp using the Integrative Genomics Viewer (IGV, 2.3.x, http://software.broadinstitute.org/software/igv/) [23, 24].

Sanger sequencing of COL2A1

Primer pairs were designed to amplify exon 42 of COL2A1 (Additional file 5). These primer pairs had product sizes of 560 bp, 785 bp and 793 bp and included the respective exon with their adjacent introns. PCR was performed in 22-μl reaction volumes containing 2 μl genomic DNA, 70 μM deoxyribonucleoside triphosphates, 5 pmol of each primer, and 3.5 U of Taq polymerase in the reaction buffer supplied by the manufacturer (Qiagen, Hilden, Germany). After a 5-min initial denaturation at 94 °C, 36 cycles of 30 s at 94 °C, 1 min at 60 °C, and 45 s at 72 °C were performed in a T Professional Thermocycler (Whatman Biometra, Goettingen, Germany). The amplicons of bulldog calves 1–3, a semen sample of Energy P, three normal progeny of the sire Energy P and a further unrelated control animal were sequenced with forward and reverse sequencing primers on a Genetic Analyzer 3500 (Applied Biosystems by Life Technologies, Darmstadt, Germany). Sequences were analysed with Sequencher 4.8 (Genes Codes, Ann Arbor, MI, USA). Sequences were screened for mutations using the bovine reference sequence from NCBI (Bos_taurus_UMD_3.1.1, GCF_000003055.6) and Ensembl (UMD3.1, GCA_000003055.3).

Validation

A PCR (mismatch)-restriction fragment length polymorphism (PCR-RFLP) was developed for the SNP g.32476082G > R. Primer pairs, amplicon size in base pairs (bp), annealing temperature, restriction enzyme and incubation temperature are given in Additional file 2.