Genotypes and animals

BovineSNP50 BeadChip (Illumina, San Diego, CA) [16] data for 3993 registered Angus animals born between 1955 and 2012 and representing 63 generations were available for analysis. Genotypes had been filtered to retain data with a call rate of ≥90% and minor allele frequency (MAF) ≥ 0.01 [2]. Pedigrees had also been validated by homozygous transmission incompatibility rates estimated as the frequency of loci for which the parent and offspring were alternate homozygotes. Pedigree relationships were expunged when the incompatibility rate exceeded 1.5%, and individuals who did not match their recorded parents or offspring were removed from the pedigree. After filtering, 52,545 loci and 3961 animals remained. These sites were next phased and missing genotypes were imputed using findhap (v2) based upon a combination of simple haplotype frequency sorting and the use of all available pedigree information [17].

Population and pedigree based haplotype analyses

Haplotypes were defined by 20-marker sliding windows genome-wide. Evidence for a haplotype harboring a lethal allele was evaluated using two statistics. The first was based on expected population frequency. For each haplotype in each window, frequency (p _i ) was estimated as the sample frequency (\( {\widehat{p}}_i \)) and the number of individuals in the population that were homozygous for the haplotype (H _i ) was tallied. When there was an absence of homozygotes, the likelihood of this occurrence was calculated using the Binomial distribution assuming Hardy Weinberg equilibrium and selective neutrality of the haplotype as P (H _i  = 0) = (1 - \( {\widehat{p}}_i^2 \))^N where N is the sample size.

In the second analysis, the actual matings within the pedigree were used to estimate the expected number of homozygotes for any given haplotype based upon the number of matings involving carriers and the assumption of selective neutrality. The pedigree for the genotyped animals was parsed to identify patrios, defined here as families for which genotypes were available for the offspring, sire and maternal grandsire [13]. Maternal genotypes were rarely available as DNA was primarily sourced from cryopreserved semen [2]. This family structure allowed the testing of segregation distortion when both the sire and maternal grandsire were heterozygous for the same haplotype. For each sliding window of 20 markers if a haplotype was never observed in the homozygous state and had a sample frequency of greater than 2%, the number of families for which the sire and maternal grandsire were both heterozygous was counted. The probability of observing at least one homozygote was then calculated based on the count of patrios and the assumption of selective neutrality. The probability that no homozygous hh haplotypes are observed in the progeny of C patrios when both the sire and maternal grand-sire are Hh heterozygotes and the h haplotype does not harbor a selected deleterious allele is (0.875–0.25q)^C where q is the frequency of the h haplotype in the population [13]. There were 2480 patrios represented in the sample. Regions of the genome with an identified deficit of homozygotes for a specific haplotype were examined for underlying genes using pybedtools [18] and the UMD3.1 genome assembly annotation run 104 [13, 14].

Generation of sequence data

To further analyze the variation within the genomic regions harboring haplotypes that were deficient for homozygotes, the whole-genome sequences of 109 animals from this population were examined [19]. These animals were selected for sequencing based upon their impact on the breed assessed by the expected numbers of genome equivalents (progeny have 0.5, grandprogeny have 0.25, etc) present within registered descendants in the population. They were sequenced with paired-end 2 × 100 bp sequence reads to an average of 27X depth of coverage of the UMD3.1 assembly with Illumina Genome Analyzer, GAII, HiSeq 2000 or 2500 instruments from two libraries with 350 bp and 550 bp average fragment sizes. FastQC was used to analyze the quality of the reads [20]. Exact duplicates were removed, and adapters were trimmed using a custom in-house Perl script. All remaining reads were error corrected using QuorUM [21]. Newly created duplicates (due to the trimming of low quality ends and correction of errors) and reads shorter than 35 bp were removed and the final data set was aligned to the UMD3.1 reference genome assembly using NextGENe 2.4.1 (SoftGenetics, LLC, State College, PA) alignment software. Reads were required to have a matching segment at least 35 bases long and 95% overall match, and a maximum of 2 bases of mismatch across the whole alignment. Up to 1000 alignments of equal likelihood were allowed genome-wide. NextGene 2.4.1 was also used for variant calling.

The sequence-derived variants for the 109 bulls were used for two purposes; identification of variants shared amongst animals identified as carriers for the putatively lethal BovineSNP50 20 SNP haplotypes, and for the imputation of the entire genotyped population to whole-genome sequence variation.

Examining carrier sequence data

Within the genomic regions identified from the marker data as containing candidate lethal haplotypes, all variants identified in our sequenced sample were analyzed for carrier concordance, implemented using python scripts. This involved examination of variants that were observed as heterozygous in all bulls predicted to be carriers of the homozygote deficient haplotype. The sample size for the sequenced animals was not sufficient to run the population allele frequency or patrio analyses.

Imputation of BovineSNP50 genotypes to whole genome sequence variation

We imputed the BovineSNP50 genotypes for this population to whole-genome sequence level variation in order to identify potential lethal variants, using the 109 sequenced animals as the reference population. We selected 24,974,785 SNPs from the full set that had been identified in these animals. We first included SNPs found in the 109 sequenced Angus bulls that were biallelic and located within coding sequence boundaries (455,585), UTRs (457,588) or that were splice site variants (20,695). These variants spanned the allele frequency spectrum but were identified at high sequence coverage. To enable imputation we also included 23,527,482 variants that had been independently identified in run 5 of the 1000 Bull Genomes project [22]. These variants represent filtered, high quality segregating sites identified from the whole genome sequences of 1578 animals from multiple taurine breeds including Angus. Genotypes called for the 24,974,785 variable sites genome-wide in the 109 registered Angus bulls were used as the reference set for whole-genome sequence imputation of the BovineSNP50 data using Fimpute [23].

The imputed genotypes were individually analyzed for the absence of homozygotes for alleles present within each of the candidate haplotyped loci using the frequency and pedigree approaches. We first identified high frequency variants for which no homozygous individuals were predicted. The pedigree analysis was also performed for all candidate variants identified in the frequency analysis. Variants identified by either of these processes were characterized using the variant effect predictor release 79 [24].

Identification of putative lethal haplotypes

Sensitivity to window size

For the purpose of identifying haplotypes that harbor candidate lethal mutations, two haplotypes that are identical by state (IBS) need only be identical for the alleles at L contiguous SNP loci, and not necessarily at the sequence level, which could include at a putative lethal locus. Conversely, two haplotypes that are identical by descent (IDB) will be identical at the sequence level except perhaps for occasional recent mutations. As L increases in size, we would expect that all haplotypes that are IBS are also IBD suggesting that they are essentially identical at the level of sequence. This means that the logic underlying the search strategy [13] is to identify haplotypes that all identical at the level of sequence and that all harbor the same lethal mutation such that the lack of observed homozygosity for a specific haplotype can be implicated as being due to the inevitability of homozygosity for the lethal mutation. Moreover, this also means that searching the sequences, within the haplotype boundaries, of animals that are all predicted to carry a candidate lethal haplotype to identify variants that are heterozygous in all of the carriers is a useful strategy to identify candidates for the lethal mutation. We evaluated the sensitivity of identification of these marker-based haplotypes to window size and concluded that a window size of 20 contiguous BovineSNP50 markers was appropriate for capturing the haplotypic diversity within the population (Fig. 1). This window size appears to discriminate between haplotypes that are identical by descent (IBD) and those that are identical by state (IBS) (See also Fig. 2). Our analysis shows that the number of common haplotypes detected rises as window size increases, but begins to plateau at 20 markers. Considering the moderate marker density of the BovineSNP50 (1 SNP per 50 kb), haplotypes that are defined by only two markers are assumed to be IBS for the purposes of analysis but likely actually represent a number of distinct haplotypes at the level of genome sequence. As the window size increases, the likelihood increases that two haplotypes found in different individuals that are IBS are also IBD and are thus concordant at the level of sequence variation. However, with large window sizes, recombination may lead to a lethal variant being present on more than one haplotype, thus decreasing the power of the analysis. Indeed, as window size increases, the overall rate of individuals homozygous for any haplotype declined. The window size selected for this study appears to achieve an appropriate balance of genome-wide homozygosity and rate of occurrence of high frequency haplotypes (Fig. 1). The 7 haplotypes reported in Table 1 were consistently detected for a range of haplotype window sizes from 10 to 40 markers.

We were also able to validate the IBD status of IBS haplotypes by an examination of the sequence data generated for animals that were predicted to be carriers of identical haplotypes. For each of the loci predicted to harbor recessive lethal haplotypes (Table 1), we identified the bulls among the 109 sequenced animals that were predicted to be carriers of each putative lethal haplotype and computed the pairwise rate of opposing homozygous sequenced sites between all pairs of carrier animals. For example, for the haplotype on chromosome 29, with 16 predicted carriers we made 120 pairwise comparisons of the average rates of alternate homozygote inconsistencies using all sequence called variant genotypes within the haplotype coordinates (Fig. 2). Amongst individuals that share a common haplotype, the alternate homozygote rate is related to the error rate for sequence-based genotype calls for heterozygotes. The rate observed in our predicted carriers within our 6 testable haplotype regions was similar to the rate for sire-son pairs. One region (Chr 1) had only one sequenced predicted carrier and could not be evaluated. The majority of the carrier pair comparisons with rates of opposing homozygous sequenced sites >0.01 were caused by the inclusion in the analysis of 3 individuals that had been sequenced to averages depths of <10X. Within the haplotype, opposing homozygous sequenced site rates were markedly lower for predicted carriers than for randomly selected individuals. When predicted carriers for a lethal haplotype were analyzed for a randomly selected 20 marker region elsewhere in the genome, they had opposing homozygote rates that were similar to those of the unrelated sire pairs. This indicates that the haplotypes generated from the BovineSNP50 data successfully identified genomic regions that were IBD at the level of the genome sequence. We also observed that this method is sensitive to the depth of sequence coverage due to the inaccurate identification of heterozygous loci as being homozygous when the alternate allele was never sequenced. In total, 16 animals had depths of sequence coverage of <10X and most of the remaining animals (70) had >20X. For a Sire-Son pair, low sequence coverage for at least one member of the pair led to a rate of homozygous inconsistency that was increased by an order of magnitude.

Analysis of sequence variation and candidate genes in bulls predicted to be lethal haplotype carriers

To identify the causal variants underlying these putatively lethal haplotypes, we first directly examined resequencing data from bulls that were predicted to be carriers of the lethal haplotypes. Among the 109 sequenced bulls up to 21 animals were predicted to be carriers of each of the BovineSNP50 putatively lethal haplotypes reported in Table 1. Within these seven genomic windows, we identified candidate variants that were never homozygous in the 109 sequenced bulls. However, none were observed to be exclusively heterozygous in the predicted carrier animals. That is, all of the alleles found to be heterozygous in all of the predicted carriers were also observed to be heterozygous in animals that did not carry the predicted lethal haplotype. A recessive lethal variant could be heterozygous in animals that were not carriers of the putatively lethal haplotype if the mutation is sufficiently old that recombination has occurred relative to the haplotype on which the mutation occurred. Table 1 reports, for each genomic region predicted to harbor a lethal haplotype, the number of variants that were heterozygous in all predicted haplotype carriers and that were never found to be homozygous in any of the 109 resequenced bulls. For instance, on chromosome 15, an intronic variant in SSRP1 was found in 9 of the 10 sequenced bulls that were predicted to carry the putatively lethal haplotype, and in 25 of all 109 sequenced Angus bulls, but was never found to be homozygous. Furthermore, this variant (Chr15: 81,825,933) was not observed to be homozygous in the 1000 bulls Run5 data. Although the variant is intronic, with no expected impact on splice site variation, it resides in an interesting candidate gene as SSRP1 is essential for mouse embryonic development [25].

This analysis was also conducted after excluding the 16 animals sequenced to low average sequence depth of coverage (< 10X). Genotypes for these animals are likely to contain false positive homozygotes at many true heterozygous sites. This resulted in the detection of additional concordant loci but again none were exclusive to the predicted carriers. The chromosome 29 locus had 12 additional variants identified that were located in genes such as CAPN, NRXN2, PACS1 and PP25B. However, none of these mutations are in exons or had other interpretations in Variant Effect Predictor (VEP) and all 12 were heterozygous in 18 animals, including 4 that were not predicted by their marker data to be carriers. Thus, the region appears to comprise a large consistent haplotype with no one particular variant being simply implicated as causal for lethality. The locus on chromosome 4 had 118 variants identified as being heterozygous in haplotype carriers following the removal of the low sequence coverage animals, but none were predicted to alter protein amino acid sequences.

The locus on chromosome 1 for which we predicted a lethal haplotype contains only one gene, GBE1 that encodes a protein that catalyzes the branching of glycogen, the main form of energy storage in the body [23]. There was only one sequenced animal in our population predicted to carry this disease and further analysis of the unique heterozygotes within this gene were therefore not feasible.

Analysis of imputed sequence variation in 3961 registered Angus

When analyzed in the 3961 animals, 2504 of the 24,974,785 imputed variants had a MAF ≥ 2%, no predicted homozygotes, and were found in ≥30 double-carrier patrios. Of the 24,974,785 variants, 147,764 had a deleterious consequence predicted by their SIFT scores [26], but only 4 were among the 2504 candidate autosomal lethal alleles. Two of these sites were in LOC521645, an olfactory receptor gene (near, but not within the chromosome 15 locus, Table 1), and another was in LOC100336589 on chromosome 18. The fourth was on chromosome 21 in GEMIN2, for which a recessive lethal embryonic phenotype occurs in mouse knockouts and is associated with the survival motor neuron complex [27]. This allele had a frequency of 7.8% following imputation into the 3961 animals and was observed in 55 double-carrier patrios without producing a homozygous progeny.

Twelve variants within GBE1, none of which were predicted to be deleterious, were at a MAF ≥ 2%, had no homozygotes predicted, and were found in ≥ 30 double-carrier patrios, with one variant present in 133 double-carrier patrios. However, only one of the sequenced animals was predicted to carry this haplotype making it difficult to use the data for the sequenced animals to exclude any heterozygous sites within this animal’s diplotype from candidacy for lethality. The haplotype on chromosome 1 spanning GBE1 was the only predicted lethal to contain any of the 2504 candidate variants identified in the analysis of the imputed whole genome sequence data.

Management strategies

Angus breeders have historically selected against recessive lethal alleles which manifest as fatal calf defects in this population [4]. Testing for known genetic conditions is now required in order to register Angus animals in the U.S and animals that are identified as carriers are excluded from registration. As new deleterious alleles are discovered, including those causing early embryonic loss, this approach will become untenable. The seven new putative autosomal recessive lethal haplotypes may now be predicted in hundreds of thousands of genotyped animals [29] and many more deleterious loci may be discovered as the number of genotyped animals increases and new lineages rise to prominence that contain yet to be detected recessive alleles. Management of the U.S. registered Angus population must shift from registration exclusion to a means of incorporating marker diagnostics into genomic selection and mate selection protocols to enable the long-term sustainability of the population.

A mate selection procedure has been implemented in the MateSel software that applies a linear weight against either the number of recessive alleles present in the progeny generation (LethalA) or the number of homozygous progeny (LethalG) [30]. In a simulation with a higher count (N = 100) for the number of lethal alleles segregating in a population than identified in this study, it was not possible to sufficiently weight in either selection scheme to eliminate embryonic mortality. This scenario may approach reality as the number of identified deleterious recessive alleles causing calf defects as well as embryonic loss increases. These simulations also suggest that implementing the LethalG strategy is a more efficient means of achieving genetic gain while reducing the frequency of recessive alleles.

Cole (2015) suggested an alternative strategy in which parent average breeding indexes are adjusted for lost progeny by adding a cost associated with recessive haplotypes [31]. This method could be extended to account for the joint impact of multiple linked or unlinked segregating loci to reduce the frequency of deleterious alleles. Counterintuitively, this study also found a zero to inverse relationship between the embryo’s realized inbreeding coefficient and its probability of being homozygous for an allele responsible for a recessive disorder. This suggests that the goal of reducing the long-term rate of accumulation of inbreeding in breeding programs may not impact the rate of embryonic loss due to the action of recessive lethal alleles. This is consistent with the observation from studies of embryonic loss and coancestry in humans [7].

Managing false positives

Managing selection based on these results requires certainty about the lethal haplotype’s effect. In this study, a total of 12,020 haplotypes were found genome-wide (including partially overlapping haplotypes) that occurred at a MAF ≥ 2% but were never found as homozygotes. Assuming random mating, a haplotype at a frequency of 2% was expected to be observed in 1.6 homozygous individuals in our Angus sample. However, the existence of non-random mating within this population could substantially decrease the likelihood of observing homozygotes, and explain why so many of these haplotypes were detected. The patrio analysis directly incorporates the matings that created the population to detect deviations from selective neutrality in progeny genotypes. However, this approach is limited by the availability of genotyped patrios. In breeds where insemination and pregnancy records are more detailed, these have been crucial for validating putative lethal haplotypes [13].

In the absence of these data, both the MateSel approach and the index adjustment could be adapted to account for the uncertainty of the lethality of the predicted allele or haplotype. In addition to validated haplotypes, and the high confidence haplotypes that we observe here, this approach could enable the incorporation of haplotypes into the selection scheme that appear to be lethal but that are at low frequency in the population. There are now Angus pedigrees in which hundreds of thousands of animals have been genotyped world-wide and an analysis of these data would likely improve the resolution of lethal haplotype detection and could also identify many rarer variants [14]. If these lethal alleles are individually rare but each individual carries many of them, recessive lethals could affect a substantial portion of pregnancies.

False negatives

Even with adequate sample sizes to ensure statistical power, there are limitations to the methods that we have employed. The approaches employed will only detect recessive alleles that are perfectly concordant with a BovineSNP50 haplotype. If a recent autosomal recessive lethal mutation has occurred on a common haplotype, the population will comprise haplotypes harboring either the lethal mutation or the wild type allele and homozygotes that include the wild type allele will be observed. In this case, very large sample sizes are required to detect deviations from the number of homozygotes expected under Hardy-Weinberg equilibrium. This could also involve restricting the analysis to different pedigree lineages for an IBS haplotype. This may explain why we failed to identify any of the recent recessive genetic defects found in the Angus breed [5]. In recent years, the alleles responsible for these defects have frequencies that have ranged from 3 to 9% [4]. However, the regions of the genome that harbor these alleles were not detected in the marker-based haplotype analysis and the causal variants were not detected in the analysis of the imputed sequence data. The likely cause of this is that the haplotypes we examined did not always indicate the presence of the defective allele. For instance, the mutation causing Neuropathic Hydrocephalus, originated in bull G A R Precision 1680 born in 1990 [32]. Consequently, there are both wild-type and deleterious versions of the BovineSNP50 haplotype on which the mutation arose segregating in the genotyped population. When larger sample sizes become available, it would be useful to repeat these analyses and test for homozygote deficiency rather than complete absence.

We have also not attempted to model mutations in loci with parent of origin affects, such as imprinting associated defects [33]. SNP array genotypes identify large heterozygous deletions as being homozygous for the alleles present on the non-deletion chromosome, which may prevent the identification of carriers for the deletion. This type of mutation has previously been associated with lethal recessive diseases in cattle [6]. Additionally, lethal alleles with incomplete penetrance will not be captured by our analyses. Other errors in genotyping, phasing or imputation also likely contribute to reductions in power to detect lethal haplotypes.

Using sequence data

One potential solution to the limitations of array genotype data is to analyze sequence-derived and/or imputed genotypes. These data may help with the management of putative recessive lethal alleles. True causal variants can be tracked more effectively than marker haplotypes. When the causative alleles have been identified, gene editing may also present an efficient means of reducing the genetic load of elite sires in a manner that is complementary to the current breeding system [34].

Sequence data also has potential advantages for the detection of recessive loci. If appropriately processed to capture SNPs, large and small indels and structural variants, they directly represent the pool of all recessive deleterious alleles. However, in practice, sample sizes have been small and the identification of large indels, particularly insertions, and complex structural variants has been challenging. Our analysis of sequence data failed to identify any candidate causative mutations in the marker-based haplotypes that were predicted to be lethal. The sample size for the sequenced animals was not sufficient to conduct the frequency or pedigree analysis with the genome-wide sequence variants. Furthermore, we did not attempt to analyze many types of complex variation, such as large indels or structural variants [35]. Large structural variants are enriched for deleterious variation but can be complex to analyze with short read data [36]. Alternative analyses of the sequence data to identify these variants or the use of methods which generate longer reads may be necessary to capture the causative variants.

In this study we did not detect a particular variant within a putative haplotype that was likely to cause a recessive lethal phenotype. However, the gene within the region on Chromosome 1, GBE1 appears quite promising. Mutations in this gene produce recessive phenotypes in mammals including horse, mouse and human [23, 24]. In the U.S. Quarter Horse population, phenotypes created by homozygotes for GBE1 mutations ranged from stillbirth to early failure to thrive, with death never occurring later than 18 weeks of age [37, 38]. Mouse knockout analysis revealed few visible or biochemical phenotypic effects in heterozygotes. Monitoring of embryonic development in homozygous knockouts revealed that deformities only occurred late in gestation and led to stillbirth or death shortly after birth. Mice with a construct with low GBE1 activity incorporated into their genome to replace the wild type allele demonstrated poor metabolic performance, and the accumulation of polyglucosan [39]. None of the mice with limited GBE1 function lived beyond 39 weeks, while all control mice survived the trial. While not all of the recessive GBE1 genotypes in other species have resulted in embryonic loss, the reduced growth associated with homozygosity for these mutations makes it unlikely that an affected animal would be selected as a sire or dam. They would therefore be highly unlikely to be included in our genotyped sample of Angus cattle. However, identifying homozygous calves from those produced by mating carriers and assaying their GBE1 functionality might be possible.

Mapping candidate variants without sequencing

Rather than generating expensive sequence data for identifying recessive lethals, two strategies might be useful: assay development and imputation. Novel variants of all classes that are detected by sequencing can readily be incorporated onto commercial genotyping platforms. This expedites fine-mapping within a known lethal haplotype in a commercial population. Variants that had predicted deleterious functional impacts based on bioinformatic analysis would also be excellent candidates for inclusion on commercial genotyping platforms.

Imputation accuracy, particularly for rare variants, may not be sufficient to identify candidate segregating recessive lethals. We previously analyzed the accuracy of imputation using variants from Run4 of the 1000 bull genomes project as a reference for our Angus BovineSNP50 genotyped population using the same imputation methods that were used in this study [19]. Our 109 sequenced animals had their BovineSNP50 genotypes imputed to the 1000 Bull Genomes Run4 sequence reference set as well as variants called directly from their whole-genome sequences. Comparing these two sets of genotypes revealed correlations between genotypes in the range of 80–90% for common variants, and 60–80% for variants with MAF ≤ 10% [19]. We would expect the causal variants underlying these lethal haplotypes to fall in the rare, more inaccurately imputed frequency class. This greatly complicates the utility of imputation for this application, and the results of imputation presented here are not appropriate for application to breeding decisions. Rather, the direct sequencing of larger samples of predicted carriers for each lethal haplotype will be required to identify the candidate causal variants. Had an interesting candidate locus with a plausible biological mechanism emerged, it would have been a good target for further confirmation and possibly immediate use. More sophisticated sequence imputation methods that provide higher imputation accuracies across the allele frequency spectrum are now becoming available. These have been used in human studies to identify individuals that are homozygous for rare variants with predicted deleterious effects [15]. Applying these methods to livestock populations with extensively described and genotyped pedigrees is feasible, and may prove useful for fine-mapping within candidate haplotypes.