A 0.5-Mbp deletion on bovine chromosome 23 is a strong candidate for stillbirth in Nordic Red cattle

Abstract

Background

A whole-genome association study of 4631 progeny-tested Nordic Red dairy cattle bulls using imputed next-generation sequencing data revealed a major quantitative trait locus (QTL) that affects birth index (BI) on Bos taurus autosome (BTA) 23. We analyzed this QTL to identify which of the component traits of BI are affected and understand its molecular basis.

Results

A genome-wide scan of BI in Nordic Red dairy cattle detected major QTL on BTA6, 14 and 23. The strongest associated single nucleotide polymorphism (SNP) on BTA23 was located at 13,313,896 bp with $- \log_{10} ({\text{p}}) = 50.63$. Analyses of component traits showed that the QTL had a large effect on stillbirth. Based on the 10 most strongly associated SNPs with stillbirth, we constructed a haplotype. Among this haplotype’s alleles, HAP_QTL had a large negative effect on stillbirth. No animals were found to be homozygous for HAP_QTL. Analysis of stillbirth records that were categorized by carrier status for HAP_QTL of the sire and maternal grandsire suggested that this haplotype had a recessive mode of inheritance. Illumina BovineHD BeadChip genotypes and genotype intensity data indicated a chromosomal deletion between 12.28 and 12.81 Mbp on BTA23. An independent set of Illumina Bovine50k BeadChip genotypes identified a recessive lethal haplotype that spanned the deleted region.

Conclusions

A deleted region of approximately 500 kb that spans three genes on BTA23 was identified and is a strong candidate QTL with a large effect on BI by increasing stillbirth.

Background

Previously, the molecular causes of genetic defects in cattle were mainly detected by examining dead or malformed calves [1–9]. The recent availability of high-throughput genotyping platforms and next-generation sequencing technologies has substantially accelerated the discovery of lethal genetic factors [10–14]. However, polymorphisms that are not associated with characteristic phenotypic manifestations, such as physically malformed calves, remain difficult to identify. A complementary approach is quantitative trait locus (QTL) mapping. For example, a QTL with a large effect on fertility in Nordic Red dairy cattle was first reported on Bos taurus autosome (BTA) 12 [15]. The molecular basis for this QTL was subsequently identified as a 660-kb deletion that causes embryonic death [16]. A whole-genome association scan in Nordic Red dairy cattle detected a major QTL for birth index (BI) on BTA23. BI is a composite index that describes the contribution of additive genetic effects to a calf’s genetic potential to be born.

In this study, we analyzed two component traits that contribute to BI (i.e. stillbirth and calving ease) for the QTL on BTA23. Our objectives were: (1) to identify which component traits are affected by the BI QTL on BTA23, and (2) to understand its molecular basis.

Methods

Animal Care and Use Committee approval was not required for this study because the phenotype and genotype data used were routinely collected as part of a breeding program, and no live animal experiments were performed.

Animals and traits

The study was carried out on three populations of Nordic Red dairy cattle (RDC) from Denmark (RDCDNK), Finland (RDCFIN) and Sweden (RDCSWE). A total of 4631 progeny-tested bulls with breeding values for BI and its component traits and genotype information were used for QTL mapping.

As dependent traits, we analyzed the de-regressed breeding values (DRP) for calving-related traits. We scanned the breeding values for BI to identify QTL. BI is a compound index that predicts a sire’s total direct additive genetic effect on calving by combining direct (calf) effects of stillbirth (SB) in the first (SBF) and later (SBL) calvings and calving ease (CE) in the first (CEF) and later (CEL) calvings. Calf size (CS) is not included in the prediction of BI by the Nordic Cattle Genetic Evaluation (NAV; http://www.nordicebv.info). However, birth weight is positively correlated with calving assistance/difficulty [17]. Therefore, we analyzed CS in the first (CSF) and later (CSL) calvings. Calves that were alive or dead within 24 h after birth were recorded as 1 or 0 (for SB), respectively. Farmers subjectively assessed and classified CE into two categories in Sweden and four categories in Denmark and Finland) and CS into four categories but this was recorded only in Denmark. We analyzed de-regressed breeding values [18] that were produced as part of the routine breeding value evaluation by the NAV. For further information on the recording and genetic evaluation for direct calving traits, see [19] and http://www.nordicebv.info.

Genotyping of bulls

All bulls were genotyped by using the BovineSNP50 BeadChip (Illumina, San Diego, CA) version 1 or 2 and genomic DNA extracted from whole blood or semen. Quality control (QC) to select single nucleotide polymorphisms (SNPs) were as follows: a minimum call rate of 85 % for samples, and exclusion of SNPs with a call rate less than 95 %, a minor allele frequency (MAF) less than 0.5 % and a significant deviation from Hardy–Weinberg proportions (HWP; p < 10⁻⁵). After QC, 43,415 SNPs were retained for analyses. Genomic positions of the SNPs were based on the UMD3.1 Bovine Genome Assembly [20].

Imputation to high-density (HD) and full-genome sequence

The 50 k SNPs were imputed to full sequence in two steps. First, 50 k genotypes for bulls were imputed to HD genotypes by using the IMPUTE2 software with default parameters, except for the effective population size (N_e) = 100 [21]. The reference population with HD genotypes consisted of 2036 bulls (902 Holstein, 735 Nordic Red and 399 Danish Jersey). The same QC parameters that were used for the 50 k chip were used for HD data. After QC, 648,219 SNPs remained for the HD chip.

In the second imputation step, using whole-genome sequencing (WGS) data from 242 dairy cattle as a reference, 12,322 bulls from three breeds (6032 Holstein, 1645 Jersey and 4645 Nordic Red) imputed to HD genotypes were further imputed to the full-sequence level by using Beagle software [22]. Sequences for the reference population used for imputed animals consisted of WGS data that were obtained by Aarhus University [23, 24] and by the 1000 Bull Genome Project run2 [14]. Chromosomes were divided into windows of about 20,000 consecutive SNPs, with an overlap of 250 SNPs at each end to minimize imputation errors at each end of the windows. For details on next-generation sequencing data and imputation steps, see [23].

Beagle v3.3.2 [22] was used to pre-phase reference data and to impute from imputed HD genotypes to the full-sequence variants. All SNPs with an imputation certainty (R²) value less than 0.9 were removed. A total of 8,938,927 SNPs remained that were distributed across the 29 bovine autosomes [25].

Statistical models for association analysis

First, a WGS scan was performed using a sire model without considering relationships between individuals, except for the sires and their sons. Based on the results from this scan, we selected a region on BTA23 for further analysis with an animal model that considered relationships between all individuals, using both single-marker and haplotype-based analyses.

Sire model for the whole-genome scan

A SNP-by-SNP analysis was carried out, in which each SNP was tested for association with the phenotype (de-regressed breeding value). The following linear mixed model (LMM) was used to estimate SNP effects:

$$y_{ij} = \mu + bx_{ij} + s_{i} + e_{ij} ,$$

where y _ij is the de-regressed breeding value of the jth son from the half-sib (sire) family i; μ is the general mean; b is the allelic substitution effect; and x _ij (ranging from 0 to 2) is the allelic dose of the jth individual for the SNP obtained by imputation. s _i is the random effect of the ith half-sib family, assumed to be normally distributed as $s_{i} \sim N\left( {0,\sigma_{s}^{2} } \right)$, where σ ² _s is the sire variance. e _ij is the random residual of son j from the half-sib family i and is assumed to be normally distributed as ${\mathbf{e}} \sim N({\mathbf{0}}, {\mathbf{W}}^{ - 1} \sigma_{e }^{2}$), where ${\mathbf{e}}$ is the vector of random residuals (e _i ); $\sigma_{e }^{2}$ is the error variance; and ${\mathbf{W}}$ is a diagonal matrix where the diagonal elements are weights of the DRP. The weight of the ith animal was estimated by w _i  = r ² _i /(1 − r ² _i ), where r ² _i is the reliability of the ijth animal’s DRP. Values of r ² _i that were higher than 0.98 were set to 0.98 to avoid having excessively large weights for sires with large numbers of progeny records.

All statistical analyses were conducted by using the software DMU [26]. The null hypothesis H ₀: $b = 0$ was tested with a two-sided t test. A SNP was considered to have a significant association with a trait if the $- \log_{10} ({\text{p}})$ was higher than 8.25 (after Bonferroni multiple-testing correction for 8,938,927 simultaneous tests), corresponding to a nominal genome-wide significance level of $p = 0.05$.

LMM analysis of the target region on BTA23

LMM analysis was used for association analyses of the target region on BTA23 for BI and its component traits (Table 1). Association between the SNP and phenotype was assessed by a single-locus regression analysis of the phenotype onto the allelic dosage for each SNP separately, using an LMM [27] as follows:

$$y_{j} = \mu + bx_{j} + u_{j} + e_{j} ,$$

where y _j is the phenotype (DRP) for the jth bull; μ is the overall mean; b is the allelic substitution effect; and x _j (ranging from 0 to 2) is the allelic dose of the jth individual for the SNP. u _j is the random polygenic effect with a joint multivariate normal distribution, ${\mathbf{u}} \sim N\left({{\mathbf{0}}, {\mathbf{A}} \sigma_{u }^{2} } \right)$, where u is the vector of polygenic effects (u _i ); ${\mathbf{A}}$ is the additive genetic relationship matrix; and $\sigma_{u }^{2}$ is the polygenic variance. e _j is the random residual for the jth animal distributed as described in the sire model.

Table 1 Most significantly associated SNPs with direct calving traits obtained from linear mixed model analyses for the target region on chromosome 23

The model was fitted by restricted maximum likelihood (REML) using the software DMU [26], which provided estimates of the fixed effects and their standard errors. Testing for the effect of a marker was done by using a two-sided t-test against a null hypothesis of H ₀: $b = 0$, with a Bonferroni correction for 8,938,927 multiple testing ($- \log_{10} \left( {\text{p}} \right) = 8.25$) similarly to the sire model.

Random haplotype model

Haplotypes for the 10 SNPs that were most significantly associated in the LMM analyses and located between 13,259,463 and 13,354,932 bp on BTA23 were extracted from Beagle’s WGS imputation output (Table 2). We used a LMM that included random polygenic effects and random effects of the haplotype, i.e. random haplotype model (RHM) following [28]. The RHM model was:

$$y_{j} = \mu + q_{h1j} + q_{h2j} + u_{j} + e_{j} ,$$

where q _h1j and q _h2j are random effects of the two haplotypes (one each from the sire and dam) carried by the jth individual, assumed to be normally distributed as $q_{hij} \sim N\left( {0,\sigma_{h}^{2} } \right)$, where σ ² _h is the variance of haplotype effects. Other terms in the RHM were the same as described for the LMM.

Table 2 Top 10 most significant SNPs selected to construct haplotypes for haplotype-based association mapping

The significance of the haplotype substitution effect was assessed by a likelihood ratio test that compared the RHM model to a null model that included mean, polygenic effect and random error terms, but no haplotype effects. Under H ₀, the test statistic has a χ ² distribution with one degree of freedom. The analysis was performed in the DMU software package [26]. A haplotype with a large negative predicted effect on BI, i.e. HAP_QTL was identified as the haplotype that carried the putative causal polymorphism for the QTL. The above LMM analysis was repeated with individual counts of HAP_QTL that were added as a fixed regression effect. Carrier status of each bull with respect to HAP_QTL was predicted for further analysis, as described below.

Analysis of calf survival as a function of haplotype carrier status of the sire and maternal grandsire

Recorded pregnancies were classified into types according to the carrier status of the sire and the maternal grandsire for HAP_QTL. For each pregnancy, we recorded whether the calf was alive or dead 24 h after birth. A recessive lethal mutation can be tested by comparing outcomes from mating types, because there is a 25 % probability that a pregnancy for which both the sire and dam are carriers will result in affected calves. A LMM was applied to test the effect of mating type on calf survival at birth.

Four classes of matings were defined according to the HAP_QTL carrier status, i.e. (I) noncarrier sire mated to the daughter of a noncarrier maternal grandsire; (II) noncarrier sire mated to the daughter of a carrier maternal grandsire; (III) carrier sire mated to the daughter of a noncarrier maternal grandsire; and (IV) carrier sire mated to the daughter of a carrier maternal grandsire. Stillbirths for heifers (first calving) and cows (later calvings) were analyzed separately. A total of 3,932,927 calving records from Denmark, Finland and Sweden were analyzed, with 2,785,085 records for mating type I, 536,614 for mating type II, 512,714 for mating type III and 98,514 for mating type IV. Numbers of records in each category are in Table 3.

Table 3 Number of observations (N) and calf survival rates (Survival) (%) for mating types according to genotypes of the sire (S) and maternal grandsire (MGS) for the 500-kb deletion

The fitted mixed model included parity and month of insemination (by year) as fixed effects and the maternal grandsire as a random effect:

$$y_{ijkl} = \mu + p_{j} + t_{k} + m_{l} + u_{i} + e_{ijkl} ,$$

where y _ijkl is an indicator of calf survival (0 for survival and 1 for stillbirth); p _j is the effect of parity; t _k is the effect of month and year of insemination; and m _l is the effect of mating type. Vector ${\mathbf{u}}$ of random grandsire effects u _i is assumed to be normally distributed ${\mathbf{u}} \sim N\left( {{\mathbf{0}},\sigma_{g}^{2} {\mathbf{A}}_{s} } \right)$, where ${\mathbf{A}}_{s}$ is the additive genetic relationship among sires of dams derived from the pedigree. ${\mathbf{e}}$ is a vector of random individual error terms $e_{ijkl}$ and is assumed to be normally distributed ${\mathbf{e}} \sim N\left( {{\mathbf{0}}, {\mathbf{I}}\sigma_{e}^{2} } \right)$. Assuming HWP, the expected proportion of conceptuses that are homozygous for HAP_QTL is equal to 1/4 of the probability that both parents are carriers, corresponding to 0, 0, $\frac{p}{{4\left( {1 + p} \right)}}$, and $\frac{{\left( {1 + p} \right)}}{{4\left( {2 + p} \right)}}$ for mating types I, II, III, and IV, respectively, where p is the frequency of the causative allele [16].

Search for chromosomal deletions within the QTL for stillbirth

A recessive lethal mutation can be caused by a chromosomal deletion [29, 30]. We searched for chromosomal deletions at the QTL using two approaches: (1) deviation from HWP, and (2) loss of genotype intensity based on SNP array data.

Deviation from HWP

If the recessive lethal effect is due to a chromosomal deletion, then an excess of homozygotes will be called for the SNPs that are located within the deleted region [16]. Imputed HD SNP genotypes in the QTL region were investigated for deviations from HWP. Among the imputed SNP genotypes, we tested for deviations of their frequencies from HWP by using a χ ² distribution with one degree of freedom. For SNPs that deviated from HWP, we investigated whether there was an excess or a deficit of heterozygotes.

Loss of genotype intensity

We had access to genotype intensity data for 243 RDCFIN bulls that had been genotyped with the Illumina BovineHD Genotyping BeadChip that included 725,293 autosomal SNPs [16]. We visually inspected the signal intensities around the QTL region using the SVS8 program (Golden Helix) for genotyped bulls. Reduced intensities were observed between SNP BovineHD2300003056 (23:12,289,281, rs109948445) and SNP BOVINEHD2300003186 (23: 12,816,660, rs132859742). An animal was considered to carry the deletion if it met two conditions: the average signal intensity (‘Log₂R ratio’) was less than 0, and 97 % of the SNPs were homozygous within the region defined by the signal intensities. Sires were categorized as carriers or noncarriers of the chromosomal deletion (CHR_DEL), and the average intensity per SNP was calculated for both groups. An individual that carries CHR_DEL will be erroneously genotyped as homozygous for all loci within the deleted region because only one allele is present. Thus, a small aberration (3 %) was allowed, to account for possible genotyping errors.

Search for recessive lethal haplotypes (absence of homozygotes) using 50 k genotype data

Absence of homozygotes for a common haplotype among live individuals is strong evidence for the presence of a recessive lethal allele [10]. Based on this, we searched for a recessive lethal allele on BTA23 segregating in Nordic Red dairy cattle. A total of 19,309 Nordic Red dairy cattle animals with 50 k genotypes, including 1119 SNPs on BTA23, were available. We discarded 136 SNPs that had a MAF less than 0.01 and 73 SNPs that had very strong deviations from HWP (p < 10⁻⁵). To avoid the influence of a deletion on phasing accuracy, we removed four SNPs that were located within the putative deleted region (between 12.2 and 12.9 Mb).

Next, we used Beagle 4.0 [31] to impute sporadic missing genotypes and to infer the haplotype phase. Finally, 901 SNPs remained on BTA23. To identify a haplotype that tagged the deleted region, we selected a 10-SNP region that spanned the putative deleted region (five SNPs upstream and five SNPs downstream) and covered a region between positions 11,818,761 and 13,122,318 bp. It was expected that one or more haplotypes in this window would be in tight linkage disequilibrium with the deleted region. If the deletion is recessive lethal, then the haplotype(s) that tag(s) it should only be present in the heterozygous condition; i.e. no live animal homozygous for the recessive lethal haplotype should be found.

Results

Genome scan for birth index

Figure 1 presents the Manhattan plot for the QTL scan for BI based on WGS data in Nordic Red dairy cattle. Additional file 1: Table S1 provides the top associated genome-wide significant SNPs for each chromosome. Strong association signals for BI were observed on BTA6, 14 and 23. Previously, we reported a QTL for BI on BTA6 which increased calf size at birth and adult stature [32] and we attributed the increases in calving difficulties and stillbirth to the increased calf size at birth. In the present study, we focused on the QTL for BI that is located on BTA23. The most strongly associated SNP was located at 13,313,896 bp (rs722178836) and had an allele substitution effect of -0.72 additive genetic standard deviations and $- \log_{10} \left( {\text{p}} \right) = 50.63$ (Fig. 2).

Fig. 1

Manhattan plot for single-variant association analysis for birth index in Nordic Red dairy cattle. Association analyses were carried out using a sire model for whole-genome sequence variants

Fig. 2

Association results for single-variant analysis, obtained using a sire model, for birth index on chromosome 23 in Nordic Red dairy cattle. The X-axis is chromosomal location in million base pairs and Y-axis is −log₁₀(p)

LMM analysis of the targeted region

BI and three calving traits (SB, CE and CS) were analyzed for association with SNPs in the target region on BTA23 using a LMM. The SNPs that were most strongly associated with each trait are in Table 1. The strongest association signal was observed for SBL (light gray dots in Fig. 3), followed by BI (light gray dots in Fig. 4) and SBF (light gray dots in Fig. 5), with a peak at 13,313,896 bp (rs722178836). Association peaks for both CE and CS were positioned at 17.63 Mbp (see Additional file 2: Figures S1, S2). The $- \log_{10} \left( {\text{p}} \right)$ values were between 15.9 and 24.4, i.e. much smaller than those observed for SB and BI.

Fig. 3

Single-variant association signal plots for calf survival at later calvings (SBL) (light gray dots) and single-variant analysis with the 10-SNP haplotype as a cofactor (blue dots). Haplotypes were constructed with the 10 top associated SNPs at the QTL peak and fitted as covariates in the model (blue dots)

Fig. 4

Single-variant association signal plots for birth index (BI) (light gray dots) and single-variant analysis with the 10-SNP haplotype as a cofactor (blue dots). Haplotypes were constructed with the 10 top associated SNPs at the QTL peak and fitted as covariates in the model (blue dots)

Fig. 5

Single-variant association signal plots for calf survival at first calving (SBF) (light gray dots) and single-variant analysis with the 10-SNP haplotype as a cofactor (blue dots). Haplotypes were constructed with the 10 top associated SNPs at the QTL peak and fitted as covariates in the model (blue dots)

Haplotype-based analysis of the targeted region

The strongest association signal among the analyzed traits was obtained for SBL (Table 1; Fig. 3). The 10 SNPs that were most strongly associated with SBL were used to construct a haplotype (Table 2). Diversity of these 10-SNP haplotypes was extremely low, with only five of the 1024 possible haplotypes being observed. Haplotypes and frequencies for the Nordic Red dairy cattle population were ATGTGAGTGA (0.9306), ATGTGAGTGC (0.0059), GCACAGCGAA (0.0019), GCACAGCGAC (0.0577) and GCACAGCGGA (0.0042).

When we analyzed haplotype as a random effect in the RHM, we observed that the haplotype GCACAGCGAC reduced the breeding values for SBL by 0.74 additive genetic standard deviations. The frequency of carriers of this haplotype (HAP_QTL) was highest in RDCFIN (15.5 %), followed by RDCSWE (10.1 %) and RDCDNK (3.9 %). No homozygous individuals for this haplotype were observed. When this haplotype was included in the model as a cofactor, no additional significant associations were observed (blue dots in Figs. 3, 4, 5), suggesting that this haplotype could fully explain the QTL variance.

Embryonic homozygosity for HAP_QTL is responsible for stillbirth

We analyzed the stillbirth rate for the four mating types (I, II, III and IV) based on the carrier status of the sire and maternal grandsire. The stillbirth rate was approximately 6 % higher for matings in which both the sire and maternal grandsire carried the putative causative haplotype HAP_QTL compared to matings in which neither was a carrier (Fig. 6). We consistently observed higher calf mortality rates from carrier-sire by carrier-maternal grandsire matings across the three Nordic countries, in both the first and later calvings (Table 3).

Fig. 6

Stillbirth rate (%) by mating type according to the carrier status of the sire (S) and maternal grandsire (MGS) for the lethal haplotype (HAP_QTL), based on the top 10 SNPs associated with stillbirth. I noncarrier S × MGS, II noncarrier S × carrier MGS, III carrier S × noncarrier MGS and IV carrier S × MGS matings, FIN Finland, SWE Sweden and DNK Denmark

Presence of a large chromosomal deletion in the QTL region

Deviation from HWP

Significant deviations from HWP were observed for SNPs on the HD chip in the interval between 12.1 and 12.9 Mbp (see Additional file 3: Figure S3), all of which were due to an excess of homozygotes for the SNPs. This result provides a strong indication of a chromosomal deletion.

HD genotype intensity confirms the structural variation

The plot of the average intensities from carriers and noncarriers of CHR_DEL (Fig. 7) showed a clear reduction in signal intensities for SNPs in the region between 12,289,281 and 12,816,660 bp, which confirmed the existence of a chromosomal deletion. Since the deleted region harbored several repetitive elements, it was difficult to identify the exact breaking points from the position of the split reads of the sequence data. Amplification across the deleted region was unsuccessful because unique primers for its extremities could not be designed. Sequence coverage increased in depth at the end region of the breaking point (see Additional file 4: Figure S4), which could be due to assembly problems or the presence of repetitive elements. These data, combined with the individual HD SNP intensity data and sequence alignments (results not shown), indicate that either the bovine genome assembly of this region is incorrect, or that the structural variation is more complex than just a simple deletion.

Fig. 7

Average genotype signal intensity per SNP for carriers (red) and noncarriers (blue) of the chromosomal deletion (CHR_DEL), plotted against their genomic positions. Carrier status of an individual was predicted based on loss of genotype intensity and homozygosity using SNPs from the Illumina BovineHD genotyping array

The deletion discovered here (CHR_DEL) is not the only structural variation found in the region. Bickhart et al. [33] published a copy number variation (CNV) in Angus cattle that originated from Great Britain within the region affected by CHR_DEL (BTA23: 12,500,344–12,515,633 bp), and Hou et al. [34] found a CNV (BTA23: 12,938,090–13,081,504 bp) in Nordic Red dairy cattle near CHR_DEL. Both CNV are gains of DNA sequence and are much smaller than the structural variation identified here.

In addition, we checked the genotype intensity (log₂R ratio) for eight SNPs on the 50 k chip located within the deleted region in a larger sample of 13,852 Nordic Red dairy cattle animals. The distribution of the intensities is in Fig. 8. Animals that were homozygous for all eight SNPs had low genotype intensities, which supports the existence of a chromosomal deletion. The length of this deletion was approximately 0.5 Mbp, as indicated by the HD intensity data and WGS depth.

Fig. 8

Density plot of genotype intensity for eight SNPs located within the chromosomal deletion in the Illumina BovineSNP50 BeadChip for an independent sample of Nordic Red dairy cattle

50 k genotype data confirm the presence of a recessive lethal haplotype spanning the deleted region

We identified a haplotype (HAP_50k) that spanned a region between 11,818,761 and 13,122,318 bp with a population frequency of 0.043, but we found no homozygous individuals for this haplotype. Under Hardy–Weinberg equilibrium, we would expect to have 36 homozygotes for this haplotype among the 19,309 50 k-genotyped animals. The probability of observing no homozygotes by chance is 2.31 × 10⁻¹⁶ (assuming a Poisson distribution), if this haplotype is not a recessive lethal variant. Absence of homozygotes for HAP_50k is strong evidence that Hap_50k either carries or tags a recessive lethal variant.

Concordance of carriers of HAP_QTL and CHR_DEL

WGS depth confirms the chromosomal deletion

We took advantage of the WGS information that was available for 83 Nordic Red dairy bulls. Six sequenced bulls carried the haplotype associated with stillbirth (HAP_QTL). All six also showed an approximately halved sequencing depth between 12.29 and 12.84 Mbp (one example shown in Additional file 4: Figure S4), whereas such a reduction was not found for any of the noncarrier bulls. None of the 59 Holstein or Jersey bulls for which WGS data was analyzed, exhibited this reduction in sequencing depth for this region on BTA23.

CHR_DEL carriers inferred from HD genotype intensity

Among the 243 RDCFIN bulls with HD genotypes, 47 animals were declared carriers based on the criteria of an average genotype signal intensity less than 0 and homozygosity for more than 134 of the 139 loci within the putative deleted region between 12.2 and 12.8 Mb. Of these 47 carriers, 25 were homozygous for all 139 loci, 13 were homozygous for 138 loci, and nine were homozygous for 134–137 loci. We examined the concordance between carriers of CHR_DEL and carriers of the recessive lethal haplotype based on 50 k data (HAP_50k). Among the 47 animals that carried the deletion, 42 carried the tagging haplotype. SNPs from the Illumina BovineSNP50 BeadChip located within the deleted region are in Additional file 5: Table S2.

Discussion

A large-effect QTL for BI was detected on BTA23 in a genome-wide association study in Nordic Red dairy cattle. A detailed study of its component traits showed that it had a large effect on stillbirth. We constructed haplotypes from the 10 SNPs that were most strongly associated with stillbirth based on WGS data and identified animals that carried the haplotype that had a negative effect on stillbirth (HAP_QTL). No homozygotes for HAP_QTL were observed. No other QTL effects were observed in the region when HAP_QTL was included in the model as a cofactor. Matings for which the sire and maternal grandsire carried HAP_QTL had significantly higher frequencies of stillbirth compared to matings for which the sire and/or maternal grandsire were a noncarrier.

We used the HD SNP genotypes to search for chromosomal deletions within the QTL region. Based on deviations from HWP and reduced genotype intensity, we observed a large chromosomal deletion (CHR_DEL) and identified animals that carried this deletion. Using a large number of animals with 50 k genotypes, we confirmed the existence of a recessive lethal haplotype (HAP_50k) that spanned the deleted region and identified carriers of HAP_50k.

After studying the concordance between animals carrying HAP_QTL, CHR_DEL and HAP_50k, an approximately 500-kb deletion on BTA23 emerged as a strong candidate QTL for stillbirth. Two lines of evidence support this conclusion. First, among the 83 sequenced Nordic Red dairy cattle animals, the six bulls that carried HAP_QTL showed a window where read depth was reduced by 50 % compared to the rest of the genome. In contrast, no noncarrier bulls showed this reduction in sequencing depth. However, uncertainty remains on the exact position of the breakpoints of the deleted region due to the presence of other structural variants in this genomic region (see Additional file 4: Figure S4). Second, analysis of 50 k genotypes for a large number of individuals confirmed a recessive lethal haplotype (HAP_50k) that spanned this deleted region.

Several observations support the conclusion that HAP_QTL is a recessive lethal haplotype that causes stillbirth. First, no homozygotes for HAP_QTL were observed. Second, the effect of HAP_QTL was observed in the case of carrier-sire and carrier-maternal grandsire matings for all three Nordic Red dairy cattle populations. There was a slight increase in stillbirths among type III matings in the Finnish population. This result may be due to a higher frequency of carriers of HAP_QTL, such that the dam inherits the deletion from its dam (i.e. maternal granddam of the calf) even when the maternal grandsire is a noncarrier. Frequencies of carriers of HAP_QTL were much lower in the Swedish and Danish RDC populations. Consequently, no effect on fertility was observed in type III matings in these populations.

The increase in the rate of stillbirths among type IV matings was smaller than the 12.5 % expected for a fully penetrant recessive lethal allele. Because stillbirths are recorded by farmers, the difference may be due to some of the stillbirths not being reported to the central registration system. It is also possible that some deaths of calves occurred during late gestation and, thus, were not recorded by farmers as stillbirths. The stillbirth is recorded as a binary trait but we analyzed it by assuming a normal distribution to test the effect of mating type on calf survival at birth. A logistic regression model might have been a better choice for the analysis of this phenotype.

Genes located within the deleted region

Three genes are located within the deleted region: (1) BTB (POZ) domain containing 9 (BTBD9, 12,405,340–12,462,065 bp), (2) glyoxalase I (GLO1, 12,483,468–12,509,232 bp) and (3) dynein, axonemal, heavy chain 8 (DNAH8, 12,550,367–12,883,939 bp) (Fig. 9). Which, if any, of these genes causes stillbirth if deleted is not obvious from their known gene functions. Phenotypes of BTBD9 knockout mice suggest that BTBD9 is involved in synaptic plasticity, learning, memory and protein alterations associated with vesicle recycling and endocytosis [35]. When glyoxalase I activity decreased in situ through aging and oxidative stress, levels of glycation and tissue damage increased. These effects may be associated with a risk of developing vascular complications due to diabetes and uremia [36]. Overexpression of glyoxalase 1 in bone marrow cells reversed defective neovascularization in streptozotocin-induced diabetic mice [37]. Diseases associated with DNAH8 include myosin storage myopathy and reduced body myopathy (http://www.genecards.org/cgi-bin/carddisp.pl?gene=DNAH8, accessed 23 December 2014). The most well-known disease associated with dynein malfunction in humans is polycystic kidney disease (PCD; also known as Kartagener syndrome), which is characterized by malfunctioning of the cilia and sperm, as well as respiratory problems. Around half of the PCD patients also have situs inversus. Extreme defects in laterality may be diagnosed as heterotaxy with asplenia or polysplenia [38].

Fig. 9

Genes located within the targeted region based on the Ensembl Bos taurus version 78 (UMD3.1, http://www.ensembl.org/)

The Arachnomelia mutation described by Buitkamp et al. [39], located near the QTL peak, was not the cause of the stillbirth cases described here. Arachnomelia syndrome in Simmental cattle is caused by a 2-bp deletion at the MOCS1 locus. This deletion was not observed among the sequenced RDC bulls. Haplotypes constructed around the MOCS1 locus (HAP_MOCS1) did not show any association with stillbirth in Nordic Red dairy cattle (see Additional file 6: Figure S5).

No other indexes in the Nordic breeding evaluation showed significant QTL near the deleted region

We carried out association studies for the other indexes included in the breeding goal for the Nordic Red dairy cattle breed (separate studies, details not presented here). The QTL for stillbirth on BTA23 was genome-wide-significant only for body conformation index, with a peak at 12,153,108 bp (rs110295087, $- \log_{10} \left( {\text{p}} \right) = 13.92$). The highest peak for body conformation on BTA23 was at 40,658,187 bp (rs133817421, $- \log_{10} \left( {\text{p}} \right) = 17.98$). No other indexes measured in Nordic Red dairy cattle showed genome-wide-significant QTL in this region. Thus, we found no evidence that this QTL for stillbirth exhibits pleiotropic effects on other traits in the breeding goal.

Kadri et al. [16] reported very high frequencies of carriers (13 to 32 %) of a large deletion on BTA12 that has a recessive lethal effect and segregates in three Nordic Red dairy cattle populations. They suggested that these high frequencies were due to the large positive effects on milk production traits. In contrast, low to moderate frequencies of carriers of HAP_QTL indicate that this haplotype has not been actively selected by hitchhiking such as the deletion on BTA12 in the Nordic Red cattle.

We used two imputation software packages, Impute2 and Beagle to impute 50 k genotypes to HD and thereafter to whole-genome sequence variants, respectively. These two imputations were performed in different projects. The performances of both software packages are very similar [40] and therefore, we did not expect to find any influence of using two software packages on our results. It is theoretically expected that a method that combines linkage information and LD at the population level is more efficient for genotype imputation and haplotype phasing in cattle for which good pedigree records are available; however, several studies have not confirmed this hypothesis (e.g. Ma et al. [41]).

We are surprised that in spite of its moderate frequency in some Nordic Red dairy cattle populations, this deletion was undetected, to date. This may be due to the lack of any obvious physical deformities of the dead calves. Examination of a calf that is homozygous for the deleted region may help to identify the gene for which loss of function constitutes the proximate cause of stillbirth. However, our efforts to obtain carcasses of dead calves from carrier–carrier matings have so far been unsuccessful.

Conclusions

In a whole-genome association scan for Nordic Red dairy cattle, we detected a QTL on BTA23 that affects BI. We dissected this QTL by analyzing the component traits of BI, observed that it had a very large effect on stillbirth and fine-mapped the QTL. We present evidence that the underlying causative variant acts as a recessive lethal variant. Based on analyses of WGS data of carrier and noncarrier individuals for HAP_QTL, an approximately 0.5-Mbp deletion/structural variation emerged as a strong candidate responsible for stillbirth.

Availability of supporting data

The sequence data of the animals are publicly available in the Sequence read Archive of NCBI (http://www.ncbi.nlm.nih.gov/sra) under accession numbers SRX527690-SRX527732 and are part of the 1000 bull genomes project (http://www.1000bullgenomes.com).