Phenotypic manifestation of the defect

Seven calves between the age of seven and 17 weeks with severe skin lesions and poor general health status were admitted to the Clinic of Ruminants. Clinical findings such as scaling and crusting of skin and secondary adhesion of hair were most evident around the muzzle, the eyes, above the sternum and on the extremities (Figure 1A,B). Skin of the inguinal region was seborrhoeic and all palpable lymph nodes were enlarged. The calves suffered from erosions in the interdigital spaces and erosive or ulcerative lesions of the oral mucosa (Figure 1C). All animals were underdeveloped in height and weight and had a history of recurring diarrhoea and pneumonia. Two of the calves were supplemented with dietary zinc (750 mg/day), but did not respond to the treatment. Due to the advanced state of the disease and with no prospect for improvement, all animals were euthanized and subjected to necropsy.In hematoxylin and eosin (H&E) stained sections, the skin of the affected animals displayed a mild to severe chronic dermatitis with serocellular crusting and partial superficial bacterial colonization. Multiple intracorneal and intraepidermal accumulations of serum as well as multiple ulcerations and epidermal necroses were present. The epidermis was severely oedematous in multiple sections (Figure 2). Dermis and epidermis were diffusely infiltrated by a moderate to high number of neutrophils. Besides moderate oedema, the dermis exhibited infiltration with lymphocytes and plasma cells to a limited extent. In the cases where the thymus was examined, decreased thymic cellularity and poorly demarcated cortex-medulla-border were observed. Further, the clinical diagnoses of enteritis and pneumonia were confirmed during the pathological examination.

In addition to the seven cases described above, another Fleckvieh calf (eleven weeks old) with severe alterations of the skin and diarrhoea was reported by a veterinarian. The described lesions were very similar to the phenotype of the previously reported seven calves. The calf received dietary zinc (750 mg/day) for a period of two weeks. However, the overall condition of the animal did not improve and the calf died at the age of 13 weeks.

Inspection of the pedigrees of the affected calves revealed a common ancestor suggesting a genetic background. On the basis of the patients’ history, the clinical and pathological findings, the calves were tentatively diagnosed to suffer from BHZD.

Analysis of SLC39A4 –the gene causing BHZD in Holstein-Friesian

Mutations in SLC39A4 are known to cause defects resembling the phenotypic appearance of the eight affected Fleckvieh calves in various species including cattle [8]. SLC39A4 is located at the proximal region of bovine chromosome 14 (BTA 14: 1,719,732 bp – 1,724,221 bp). The gene was re-sequenced in a case–control panel consisting of all affected animals, all available dams and sires and randomly selected, unaffected control animals. Totally ~7 kb of genomic sequence was screened, resulting in the detection of ten SNPs (Additional file 1). The mutation causing BHZD in Holstein-Friesian was not present in the diseased animals and none of the detected polymorphisms was associated with the disease phenotype, nor was any of the polymorphic sites compatible with the supposed pattern of recessive inheritance.

Identification of the disease-associated region

Since the analysis of SLC39A4 did not reveal a potentially causal mutation, we applied an array-based approach to identify the underlying genomic region. The eight affected calves together with 1,339 unaffected Fleckvieh bulls were genotyped with the Illumina BovineHD BeadChip. A genome-wide association study using genotypes of 644,450 SNPs revealed a strong association signal on BTA 21. Eighty-two SNPs located within an 18.19 Mb interval from 53,140,245 bp to 71,333,740 bp were significantly associated (P < 7.88 × 10^-8) (Figure 3A). The most prominent association signal (P = 5.87 × 10^-89) was observed for BovineHD4100015383, located at 69,873,257 bp. As expected, the proximal region of BTA 14 harbouring SLC39A4 did not show association at all.

Autozygosity mapping within the distal region of BTA 21 revealed a common 1,023 kb segment of extended homozygosity in the eight affected animals (70,550,045 bp – 71,573,501 bp) (Figure 3B). However, one out of 1,339 animals of the control group was homozygous for the disease-associated region as well. Signal intensities obtained from high-density genotyping revealed no indication for the presence of large structural variants (deletions, copy number variations) within the associated region (Additional file 2).

The frequency of the associated haplotype was estimated in a sample of 10,355 unaffected Fleckvieh animals. Among them, 380 were heterozygous carriers and three were homozygous, yielding a frequency for the associated haplotype of 1.86%. The haplotype distribution shows no deviation from the Hardy-Weinberg equilibrium (P = 0.75).

Analysis of CRIP1 and CRIP2 – genes with similar function as SLC39A4

The segment of extended homozygosity contains 17 transcripts and protein encoding genes, respectively (UMD3.1 annotation [10], Figure 3C). Two of these genes encode cysteine-rich proteins (CRIP1 and CRIP2). CRIPs are highly expressed in the intestine [11] and are essential for zinc absorption as well as for immune system related cytokine regulation [12, 13]. Thus, they represent excellent candidate genes for the observed disease phenotype.

CRIP1 and CRIP2 (71,390,596 bp – 71,392,110 bp and 71,365,858 bp – 71,382,666 bp, respectively) were re-sequenced in the case–control panel. We screened ~6 kb of the genomic sequence of CRIP1, resulting in the detection of 18 SNPs and two insertion/deletion (InDel) polymorphisms (Additional file 1). Approximately 8 kb of the genomic sequence of CRIP2 were screened, resulting in the identification of 30 SNPs (Additional file 1). None of the variants in CRIP1 and CRIP2 was compatible with the disease phenotype. Thus, variation in the two positional and functional candidate genes is likely not causal for the observed disease.

Identification of the underlying mutation by exploiting whole-genome sequencing data

The compatible InDel-polymorphism (rs381259516) is also segregating among 191 non-Fleckvieh animals (128 Holstein, 15 Jersey, 48 Angus), which have been sequenced in the context of the 1000 bull genomes project [15]. No regulatory or functional consequence was predicted for the compatible SNP (rs385301007) in intron 30 of LOC100299595. Manual re-annotation of INF2 revealed that the presumed missense mutation (rs384306864) is not located within the coding region of the gene (Additional file 3). Therefore, these three variants were excluded as being causal for the described phenotype. Only a point mutation in exon 6 of PLD4, resulting in a premature stop codon (c.G645A, p.W215X, BTA 21:71,001,232 bp, rs378824791) (Figure 4, Additional file 3), was retained as candidate causal mutation. The resulting protein is shortened by 273 amino acids and lacks essential domains for enzymatic activity [16, 17]. The mutation was confirmed by Sanger sequencing (Additional file 4).

Validation of the p.W215X-mutation in PLD4

Survival rate of progeny descending from bulls carrying the p.W215X-mutation

In order to study the effect of the deleterious allele on the population level, we estimated the survival rate of the descendants from different mating types with regard to the p.W215X-genotype in sires and maternal grandsires. If both, sire and maternal grandsire are heterozygous for the p.W215X-mutation, the probability of the resulting calf being homozygous is 12.5%. The survival rate of calves descending from such matings (N = 1,213) is significantly lower (P = 2.971 × 10^-8) compared to the survival rate of calves descending from non-risk matings across all age classes, as expected in the case of recessive inheritance of the deleterious allele (Figure 5). At day 300, 16.63% of the calves resulting from risk matings have perished, while the mortality of calves descending from non-risk matings is 9.66% only (Table 3).

Animal ethics statement

Semen samples were collected by approved commercial artificial insemination stations as part of their regular breeding and reproduction measures in cattle industry. The collection of blood samples was carried out by trained veterinarians during treatment of affected animals following standard veterinary protocols in Germany. No ethical approval was required for this study.

Animals and DNA extraction

In order to confirm a genetic predisposition to the observed phenotype, a case–control panel consisting of the eight affected calves, their dams and sires and twelve unaffected control animals was set up. Seven affected animals were identified and examined by veterinarians of the Clinic for Ruminants. One affected calf was reported by a herd veterinarian. Blood samples of the affected calves and their dams were collected by trained veterinarians following standard procedures and DNA was extracted using proteinase K digestion and salt-out extraction. For the control group, genomic DNA was prepared from semen straws following standard protocols using proteinase K digestion and phenol-chloroform extraction.

Annotation and polymorphism screening of SLC39A4, CRIP1 and CRIP2

The GENOMETHREADER software tool [33] was used to predict the genomic structure and localization of SLC39A4, CRIP1 and CRIP2 based on the University of Maryland UMD3.1 assembly of the bovine genome sequence [10] and the Dana–Farber Cancer Institute bovine gene index release 12.0 [34] together with the annotated RNA sequences of the UMD3.1 assembly [10]. The GENOMETHREADER output was viewed and edited using the Apollo sequence annotation editor [35]. The genes were PCR amplified (the primers are listed in Additional file 5), including exons, introns and the flanking 3'- and 5'-regions. Sequencing reactions were done using BigDye^® Terminator v1.1 Cycle Sequencing Kit (Applied Biosystems, Foster, CA., USA; Life Technologies Corporation, USA). Electrophoresis of purified sequencing reactions was performed on the ABI 3130x1 Genetic Analyzer (Applied Biosystems, Foster, CA., USA; Life Technologies Corporation, USA). The Phred/Phrap/Polyphred software suite [36–38] was used for base calling, sequence alignment and polymorphism detection. Sequences were viewed with consed[39].

High- density genotyping and quality control

The eight affected calves and 1,339 unaffected Fleckvieh bulls were genotyped with the Illumina BovineHD BeadChip comprising 777,962 SNPs. Genotype calling was performed using default parameters of Illumina's BeadStudio. The chromosomal position of the SNPs was determined based on the UMD3.1 assembly of the bovine genome [40]. 343 mitochondrial, 1,224 Y-chromosomal and 1,735 SNPs with unknown chromosome position were not considered for subsequent analyses. Quality control was carried out with PLINK v1.07 [41]. 19,734 SNPs for which genotyping failed in more than 5% of the individuals and 110,746 SNPs with minor allele frequency < 0.5% were omitted. Sporadically missing genotypes were imputed using Beagle genetic analysis software [42]. The final dataset comprised 1,347 animals (eight affected calves and 1,339 controls) and 644,450 SNPs.

Genome-wide association study

To account for population stratification and the resulting inflation of false positive associations, a mixed model based association analysis was performed. We used GEMMA[43] to fit the linear mixed model y = μ + xb + Zu + e, where y denotes the affection status (coded as 1 and 2 for affected and unaffected animals, respectively), μ is the intercept, x is a vector of marker genotypes, b is the SNP effect, Z is an incidence matrix, u is a vector of random polygenic effects ~ N(0, G${\sigma }_a^2$), where ${\sigma }_a^2$ is the additive genetic variance and G is the genomic relationship matrix (GRM) among the 1,347 animals built based upon 627,627 autosomal SNPs following VanRaden's approach [44] and e is a vector of error terms.

Haplotype analysis

Haplotypes of 1,347 animals of the initial genome-wide scan were inferred using default parameters of Beagle genetic analysis software (see above). To assess the frequency of the associated haplotype in a larger sample of the Fleckvieh population, array-derived genotypes of another 9,016 unaffected animals were analysed. Of the 9,016 animals, 7,000 were genotyped with the BovineSNP50 BeadChip (50 K) and 2,016 were genotyped with the Illumina BovineHD BeadChip (777 K). High-density genotypes and haplotypes of the 50 K data set were inferred based on haplotypes of the 777 K data set using a combination of Beagle (see above) and Minimac[45], which yields high imputation accuracy in cattle [46].

Exploiting whole-genome re-sequencing data for mutation screening

The genomes of 43 animals of the Fleckvieh population were sequenced to an average coverage of 7.46x, among them an affected calf and a healthy animal being homozygous for the disease-associated haplotype with coverages of 7.8x and 6.2x, respectively. Sequencing on Illumina GA IIX and Hiseq 2000 instruments was performed as detailed by Jansen et al.[14]. Paired-end reads were obtained and mapped to the bovine reference sequence (see above) using the Burrows-Wheeler Aligner (BWA) [47]. PICARD (http://picard.sourceforge.net) was used to mark PCR-duplicates. Subsequent multi-sample variant calling with mpileup[48] yielded genotypes at 17.17 million sites. Beagle phasing and imputation was applied to improve the primary genotype calls. A detailed overview of the entire variant calling pipeline and all obtained variants is presented in Jansen et al.[14]. Of 17.17 million sites, 7,086 SNPs and 574 InDels were located within the 1,023 kb segment (70,550,045 bp – 71,573,501 bp) of extended homozygosity on BTA 21. To account for inaccurately genotyped variants due to the low-coverage sequence data (e.g., mis-calling of heterozygous genotypes for rare variants [49, 50]), we filtered for variants that were segregating (heterozygous or homozygous for the non-reference allele) in the affected animal, heterozygous or homozygous for the reference allele in the healthy animal being homozygous for the disease-associated haplotype and homozygous for the reference allele in 41 healthy animals.

5'-exonuclease diagnostic assay of the c.G645A-PLD4mutation

A 5'-exonuclease assay was developed to obtain genotypes for the nonsense mutation (c.G645A, p.W215X, Chr21:71,001,232 bp, rs378824791), using 5'-GGG CGG CAC ATC TAC GT-3' and 5'-CCA GGG CGG ACG AAC TC-3' as PCR primers, and 5'-CAT GGA CTG GCG GTC C-3' (wild type G allele) and 5'-CAT GGA CTG ACG GTC C-3' (mutant A allele) as probes (TaqMan^®, Life Technologies Corporation, USA). Reactions were carried out on an ABI7500 Real-Time PCR system (Life Technologies Corporation, USA) using standard procedures and analysed using the allelic discrimination endpoint analysis mode of the 7500 software package v2.0.5. Genotypes for the polymorphism were obtained in 3,650 animals representing three different breeds (Fleckvieh, Holstein-Friesian, Braunvieh).

Survival analysis

The survival rate of calves from two different mating types was estimated using a Kaplan-Meier estimator as implemented in the R package ('survival') [51]. For the control group the survival rate of 2,552 calves descending from matings where the sire does not carry the p.W215X-mutation whereas the maternal grandsire is a carrier S5(non-risk matings) was estimated. The case group comprised 1,213 calves descending from matings where both sire and maternal grandsire carry the p.W215X-mutation (risk-matings).