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

Sahana et al. Genet Sel Evol (2016) 48:35 DOI 10.1186/s12711-016-0215-z RESEARCH ARTICLE A 0.5-Mbp deletion onbovine chromosome 23 is a strong candidate forstillbirth inNordic Red cattle Goutam Sahana1* , Terhi Iso‑Touru2, Xiaoping Wu1, Ulrik Sander Nielsen3, Dirk‑Jan de Koning4, Mogens Sandø Lund1, Johanna Vilkki2 and Bernt Guldbrandtsen1 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 −log10(p) = 50.63 . Analyses of component traits showed that the QTL had a large effect on stillbirth. Based on the 10 most strongly asso‑ ciated SNPs with stillbirth, we constructed a haplotype. Among this haplotype’s alleles, HAPQTL had a large negative effect on stillbirth. No animals were found to be homozygous for HAPQTL. Analysis of stillbirth records that were cat‑ egorized by carrier status for HAPQTL 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. © 2016 Sahana et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/ publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. 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 highthroughput 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. Open Access G enetics S election Evolution *Correspondence: [email protected] 1 Department of Molecular Biology and Genetics, Center for Quantitative Genetics and Genomics, Aarhus University, 8830 Tjele, Denmark Full list of author information is available at the end of the article Page 2 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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−5). 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 (Ne)=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 (R2) 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 (deregressed breeding value). The following linear mixed model (LMM) was used to estimate SNP effects: where yij 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 xij (ranging from 0 to 2) is the allelic dose of the jth individual for the SNP obtained by imputation. si is the random effect of the ith half-sib family, assumed to be normally distributed as si ∼N  0, σ 2 s , where σs 2 is the sire variance. eij is the random residual of son j from the half-sib family i and is assumed to be normally distributed as e ∼ N (0 , W −1 σ 2 e ), yij =µ+ bxij + si + eij, Page 3 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 where e is the vector of random residuals (ei); σ2 e is the error variance; and W is a diagonal matrix where the diagonal elements are weights of the DRP. The weight of the ith animal was estimated by wi=ri 2/(1−ri 2), where ri 2 is the reliability of the ijth animal’s DRP. Values of ri 2 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 H0: b=0 was tested with a two-sided t test. A SNP was considered to have a significant association with a trait if the −log10(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: where yj is the phenotype (DRP) for the jth bull; μ is the overall mean; b is the allelic substitution effect; and xj (ranging from 0 to 2) is the allelic dose of the jth individual for the SNP. uj is the random polygenic effect with a joint multivariate normal distribution, u ∼N  0,Aσ 2 u , where u is the vector of polygenic effects (ui); A is the additive genetic relationship matrix; and σ2 u is the polygenic variance. ej is the random residual for the jth animal distributed as described in the sire model. yj = µ + bxj + uj + ej, 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 H0: b=0 , with a Bonferroni correction for 8,938,927 multiple testing ( −log10 ( p )= 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: where qh1j and qh2j 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 qhij ∼N  0, σ 2 h , where σh 2 is the variance of haplotype effects. Other terms in the RHM were the same as described for the LMM. 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 H0, the test statistic has a χ2 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. HAPQTL was identified as the haplotype that carried the putative causal polymorphism for the QTL. The above LMM analysis was repeated with individual counts of HAPQTL that were added as a fixed regression effect. Carrier status of each bull with respect to HAPQTL 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 HAPQTL. 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 HAPQTL carrier status, i.e. (I) noncarrier sire mated to yj=µ+qh1j+qh2j+uj+ej, Table 1 Most significantly associated SNPs with direct calving traits obtained fromlinear mixed model analyses forthe target region onchromosome 23 Suffix F and L are for first and later calvings Effects and SE are given in breeding index units BI birth index, SB stillbirth, CE calving ease, CS calf size, MAF minor allele frequency, SE standard error Trait Top associated SNP MAF Allele substitution effect SE p value BI Chr23:13,313,896 0.059 −8.76 0.45 4.79e−80 SBF Chr23:13,313,896 0.059 −7.35 0.45 8.09e−58 SBL Chr23:13,313,896 0.059 −10.53 0.44 1.00e−120 CSF Chr23:17,631,540 0.160 1.55 0.19 1.17e−16 CSL Chr23:17,631,540 0.160 1.60 0.18 6.27e−19 CEF Chr23:17,631,540 0.160 −1.86 0.21 3.34e−19 CEL Chr23:17,581,154 0.145 −2.52 0.24 1.26e−25 Page 4 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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. The fitted mixed model included parity and month of insemination (by year) as fixed effects and the maternal grandsire as a random effect: where yijkl is an indicator of calf survival (0 for survival and 1 for stillbirth); pj is the effect of parity; tk is the effect of month and year of insemination; and ml is the effect of mating type. Vector u of random grandsire effects ui is assumed to be normally distributed u ∼N  0,σ2 gAs  , where As is the additive genetic relationship among sires of dams derived from the pedigree. e is a vector of random individual error terms eijkl and is assumed to be normally distributed e ∼N  0,Iσ 2 e . Assuming HWP, the expected proportion of conceptuses that are homozygous for HAPQTL is equal to 1/4 of the probability that both parents are carriers, corresponding to 0, 0, p 4 ( 1 + p) , and (1+p) 4 ( 2 + p) 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 χ2 distribution with one degree of freedom. For SNPs that deviated from HWP, we investigated whether there was an excess or a deficit of heterozygotes. yijkl =µ+pj+tk+ml+ui+eijkl, 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 (‘Log2R 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 (CHRDEL), and the average intensity per SNP was calculated for both groups. An individual that carries CHRDEL 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−5). 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. Page 5 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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 −log10 (p) = 50.63 (Fig.2). 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 −log10 (p) values were between 15.9 and 24.4, i.e. much smaller than those observed for SB and BI. 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 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 −log10(p) Page 6 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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 (HAPQTL) 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 Table 2 Top 10 most significant SNPs selected toconstruct haplotypes forhaplotype-based association mapping Results are for calf survival at later lactations Effects and SE are given in breeding value index units MAF minor allele frequency, SE standard error SNP MAF Allele substitution effect SE p value Chr23:13,259,463 0.062 −10.15 0.43 4.11e−115 Chr23:13,260,528 0.062 −10.15 0.43 4.05e−115 Chr23:13,260,771 0.062 −10.15 0.43 4.02e−115 Chr23:13,261,303 0.062 −10.15 0.43 4.03e−115 Chr23:13,264,527 0.062 −10.15 0.43 4.14e−115 Chr23:13,269,718 0.062 −10.15 0.43 3.48e−115 Chr23:13,270,817 0.068 −10.22 0.44 2.39e−115 Chr23:13,271,582 0.062 −10.15 0.43 3.48e−115 Chr23:13,313,896 0.059 −10.56 0.44 1.00e−120 Chr23:13,354,932 0.066 −9.84 0.42 3.30e−113 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 con‑ structed 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 con‑ structed with the 10 top associated SNPs at the QTL peak and fitted as covariates in the model (blue dots) Page 7 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 observed (blue dots in Figs.3, 4, 5), suggesting that this haplotype could fully explain the QTL variance. Embryonic homozygosity forHAPQTL 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 HAPQTL 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). 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 CHRDEL (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. The deletion discovered here (CHRDEL) 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 CHRDEL (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 CHRDEL. Both CNV are gains of DNA sequence and are much smaller than the structural variation identified here. In addition, we checked the genotype intensity (log2R 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 Fig. 6 Stillbirth rate (%) by mating type according to the carrier status of the sire (S) and maternal grandsire (MGS) for the lethal haplotype (HAPQTL), 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 Page 8 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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. 50k genotype data confirm the presence ofa recessive lethal haplotype spanning the deleted region We identified a haplotype (HAP50k) 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−16 (assuming a Poisson distribution), if this haplotype is not a recessive lethal variant. Absence of homozygotes for HAP50k is strong evidence that Hap50k either carries or tags a recessive lethal variant. Concordance ofcarriers ofHAPQTL andCHRDEL 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 (HAPQTL). 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. CHRDEL 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 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 C carrier, NC non carrier, S sire, MGS maternal grand sire, FIN Finland, SWE Sweden, DNK Denmark Mother Mating type Total S‑NC×MGS‑NC (I) S‑NC×MGS‑C (II) S‑C×MGS‑NC (III) S‑C×MGS‑C (IV) N Survival N Survival N Survival N Survival DNK‑heifer 57,165 94.3 1971 94.9 2014 94.6 73 93.2 61,223 DNK‑cow 84,679 96.9 2551 96.7 3904 96.9 176 94.3 91,310 FIN‑heifer 535,406 94.1 120,843 93.6 105,221 93.0 21,459 88.3 782,929 FIN‑cow 805,528 96.3 178,770 96.0 178,395 95.3 36,778 90.4 1,199,471 SWE‑heifer 496,779 95.8 53,709 95.3 44,785 95.7 3250 90.7 598,523 SWE‑cow 828,052 96.7 97,536 96.5 81,253 96.4 8373 91.9 1,199,471 Total 2,785,085 95.9 536,614 95.6 512,714 95.3 98,514 90.5 3,932,927 Fig. 7 Average genotype signal intensity per SNP for carriers (red) and noncarriers (blue) of the chromosomal deletion (CHRDEL), 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 Page 9 of 12 Sahana et al. Genet Sel Evol (2016) 48:35 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 CHRDEL and carriers of the recessive lethal haplotype based on 50k data (HAP50k). 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 (HAPQTL). No homozygotes for HAPQTL were observed. No other QTL effects were observed in the region when HAPQTL was included in the model as a cofactor. Matings for which the sire and maternal grandsire carried HAPQTL 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 (CHRDEL) 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 (HAP50k) that spanned the deleted region and identified carriers of HAP50k. After studying the concordance between animals carrying HAPQTL, CHRDEL and HAP50k, 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 HAPQTL 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 (HAP50k) that spanned this deleted region. Several observations support the conclusion that HAPQTL is a recessive lethal haplotype that causes stillbirth. First, no homozygotes for HAPQTL were observed. Second, the effect of HAPQTL 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 HAPQTL, 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 HAPQTL 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, Fig. 8 Density plot of genotype intensity for eight SNPs located within the chromosomal deletion in the Illumina BovineSNP50 Bead‑ Chip for an independent sample of Nordic Red dairy cattle