Genetic variation is a key contributor to health and disease, and understanding the link between an individual’s genotype and the corresponding phenotype is a major goal of genetic research (Genomes Project et al., 2010). Whole genome sequencing (WGS) within and across populations enables highly efficient variant discovery and elucidation of the molecular nature of virtually all genetic variation, from single nucleotide polymorphisms (SNPs) to copy number variants (CNVs) and other large structural variants (Sudmant et al., 2015). Large-scale studies of genetic variation in humans have dramatically improved our understanding of genetic variation across a species and within populations (Genomes Project et al., 2012). Large-scale variant catalogs establish patterns of variation across the genome, including non-coding regions (Mu et al., 2011), permitting elucidation of regional variability in mutation and recombination rates. Knowledge of the background genetic “noise” helps to decipher the link between genotype and phenotype—allowing filtering of likely neutral variants from potentially deleterious variants in genomic regions of interest or within biologic candidate genes.

The equine reference genome (Wade et al., 2009) has provided a key basis for genetic investigations in horse populations (Rebolledo-Mendez et al., 2015; Raudsepp et al., 2019). However, despite consistent progress in our understanding of genetic disease in the horse, disease-causing variants have been identified for less than 20% of the currently recognized genetic diseases (Online Mendelian Inheritance in Animals, OMIA). Similar to humans, many of the significant GWAS regions of interest found for equine traits have been in non-coding regions of the genome. The unknown function of these regions has been a barrier to pinpointing and confirming the causal variant, which is necessary to develop effective genetic tests or targeted treatments. A more complete account of “normal” genetic variation in the horse is critical to establishing the link between genotype and phenotype (Yngvadottir et al., 2009; Genomes Project et al., 2010).

Here we provide the first large-scale catalog of genetic variation in the horse derived from WGS, with the intent of describing variant numbers, types, allele frequencies, and genomic location distribution across the general population and within and across major breeds.

This report comprises the first large-scale catalog of genetic variation developed for the horse, a species with potential as a translational model for many athletic phenotypes. We used WGS of 534 horses to determine overall genetic variation in the general equine population as well as 10 individual breeds, report variants with population and breed MAF discrepancies, identify variants with no homozygotes as well as variants that are present in all individuals, and in addition identify genomic regions with high or low genetic variation.

The significant association between the number of variants identified and the depth of coverage was unsurprising, as it has long been recognized that deeper coverage improves variant calling accuracy (Li R. et al., 2009). The nonlinear association between depth of coverage and number of variants identified is particularly relevant to future genetic studies, as there is minimal to no gain in the numbers of variants detected by increasing depth of coverage over 10X.

We elected to use the intersect of two commonly used callers based on evidence from previous work showing that specificity of identified variants could be improved (Field et al., 2015). This method identified 29,882,273 variants, which is higher than the 25,800,000 variants previously identified in a smaller cohort of 88 horses using only GATK haplotype caller (Jagannathan et al., 2019b). This difference is likely due to the increased number of horses and inclusion of more genetically distinct breeds in our study, increasing our ability to identify rare variants down to a MAF across breeds of 0.0009 compared with 0.0057 in the earlier study. Additionally, we use the intersect of two variant callers (GATK HaplotypeCaller and BCFtools) to improve specificity of variants identified in the previous publication (Jagannathan et al., 2019b). While this likely would have led to a reduced sensitivity, we were still able to identify an additional 4,082,273 variants.

The 1.54 heterozygous variants per kb of sequence is similar to that reported in cattle (1.44 per kb) (Daetwyler et al., 2014), but is higher than reported in the Yoruba (1.03 per kb) and European human populations (0.68 per kb) (Genomes Project et al., 2010). The cause of this higher variation in the horse is likely related to the heterogeneity of this population, which included 46 different breeds, compared to three cattle breeds (Daetwyler et al., 2014) and only a single population in the human studies (Genomes Project et al., 2010). Given the limited genetic diversity (Petersen et al., 2013) of the horse compared with human populations this is still somewhat somewhat surprising. However, a recent study of effective population size in the horse suggests that several breeds have larger effective population sizes (Beeson et al., 2019) than reported in human populations (Tenesa et al., 2007), and we would therefore, expect to see increased heterozygosity. Another reason for the increased number of heterozygous variants in horses is likely to be related to errors in the reference genome which, unlike the human reference genome that is based on multiple individuals, is based on a single horse (Kalbfleisch et al., 2018). In the original paper from the 1,000 human genomes consortium, it was concluded that a site where every individual was homozygous for an allele not present in the reference genome was a reference genome error. This accounted for ∼1 error per 30 kb of sequence. In this study, we identified 114,733 variants that were present in every individual in this population, which are presumed to be related to an error in the reference genome or a true rare variant that is present in the individual horse sequenced for the reference genome. The length of the RefSeq equine reference genome is 2,506.95 MB, therefore, we would expect to see ∼1.37 errors per 30 kb of sequence, which may partially explain the increased number of variants in horses compared with humans.

Unsurprisingly, the degree of missingness per individual was negatively correlated with depth of coverage. This was not a linear correlation however, and beyond 10X coverage, there was minimal improvement in the degree of missingness, suggesting that 10X coverage is a reasonable target for population scale sequencing projects. The average missingness per individual varied greatly. Most of the horses with missingness >0.20 were horses that were not included in the breed analysis owing to being in the other breed category. Six of the Shetland ponies had missingness >0.20 and these were ponies with a targeted depth of coverage of ∼ 6X. The degree of missingness across the 32 autosomes and X chromosomes also varied, with the highest degree of missingness on chromosomes 12 and X. This may be related to larger uncharacterized regions on the X chromosome compared to most autosomes.

The TsTv and hetNRhom ratios are frequently utilized for quality control of sequencing data (Guo et al., 2013). The TsTv ratio here was similar (1.94) to reports of the expected TsTv ratio from genome sequencing data in humans of ∼2 (Genomes Project et al., 2010), which is a good indicator of SNP quality (Guo et al., 2013; Wang et al., 2015). The slight reduction in our study compared with human studies is likely related to decreased genetic diversity and smaller effective population sizes in horses (Petersen et al., 2013) compared with humans, as well as a lower quality reference genome (Wade et al., 2009). The TsTv ratio varies across the genome, but in humans does not vary based on ancestry (Wang et al., 2015). In the 10 target breeds, we did find that the TsTv ratio varied significantly but marginally by ancestry. This may be related to the varying depth of sequence coverage, which was not uniform across breeds. The hetNRhom ratio (2.24) was higher than the expected value of 2 based on Hardy-Weinberg equilibrium (Guo et al., 2013) and there were breed differences. This is consistent with human (Wang et al., 2015) and canine (Jagannathan et al., 2019a) reports that the hetNRhom ratio varies by ancestry. In humans, the highest median hetNRhom ratio was 2.0 in African populations with the lowest ratio of 1.4 in Asian populations; none of the populations investigated had a median hetNRhom ratio >2.0 (Wang et al., 2015). However, at least one dog breed had hetNRhom ratio of 3.3 in the catalog of canine genetic variation (Jagannathan et al., 2019a). This is likely related to the increased levels of inbreeding in horses compared with most human populations.

The variant totals differed by breed, which is consistent with reports in different cattle breeds (Daetwyler et al., 2014) and regional human populations (Genomes Project et al., 2010). Previously, this has been related to effective population size (Daetwyler et al., 2014), and we did see a significant association (p < 0.0001) between the number of variants identified in each of the 10 target breeds and a report of effective population size using 54K SNP array data (Petersen et al., 2013). However, this association was not seen with a more recent estimate of effective population size that used imputed genome-wide SNP data (Beeson et al., 2019). This may partially be related to different breeds studied, as the Beeson et al. (2019) paper did not include effective population size estimates for the Shetland and Clydesdale breeds. Additionally, we are not accounting for the degree of relatedness between the breeds studied and the reference genome. In this study, the Thoroughbred (which is the EquCab3 reference genome breed) has the fewest variants compared to the other breeds, consistent with both the reference genome being from a Thoroughbred and having the smallest effective population size (1,784) in Beeson et al. (Beeson et al., 2019). However, the Quarter Horse, which has the largest population size (6,516) (Beeson et al., 2019), but is more related to the Thoroughbred (Petersen et al., 2013) than other breeds, has a number of variants that is closer to the median, consistent with its close relatedness to the Thoroughbred (Petersen et al., 2013), but inconsistent with the large effective population size (Beeson et al., 2019). This would suggest, that while the number of variants does vary by breed, as seen in different human regional populations, the relationship between the horse breed and the reference genome appears to have had an effect in this population.

A large number of variants were shared by at least one other breed, which is not surprising given the close relatedness of the breeds investigated (Petersen et al., 2013). However, there were also multiple variants with large allele frequency discrepancies between breeds. We defined a minor allele frequency <3% as rare and >10% as common due to limitations in the number of horses in each breed investigated, rather than values used in most human studies (<0.5% for rare variants and >5% for common (Bomba et al., 2017)). With the 534 horses here our power to detect variants present at a minor allele frequency of 3 and 0.5% in the population is 1.00 and 0.93, respectively. However, it is important to note that for the breed analysis with the least number of horses in the Clydesdales (19) our power to detect these allele frequencies is only 0.44 and 0.01, respectively. To have a power greater than 0.8 to detect all rare variants <3% or <0.5% allele frequency within a breed we would need to sequence 55 and 325 horses within that breed, respectively. We therefore, had 80% power to detect the rare variants (minor allele frequency 3%) in Quarter Horses, Shetlands and Thoroughbreds. The number of variants with marked population discrepancies was quite high, with ∼35% of variants considered rare in one breed and common in another. The reason for large allele frequency differences between populations is thought to be related to genetic drift (Hofer et al., 2009) in humans. However, given that different horse breeds have been selectively bred for different traits (Avila et al., 2018), it is also likely that selection at least partially accounts for some of the variants with large allele frequency discrepancies between breeds. Only ∼5% of variants were unique to a single breed, with about 10% of these being in coding regions. The ∼5% of unique variants in horses is also lower than seen in cattle where ∼31% of variants are unique to a single breed (Daetwyler et al., 2014).

Variants with no homozygotes were explored to determine if there were variants present in the general population that could be embryonic lethal in homozygous form. 2,888/2,889 of these variants were present at a minor allele frequency in the general population greater than would be expected for a homozygous lethal disease (MAF > 0.10). The one variant that was rare had a MAF of 0.03 and was an intergenic in frame deletion. Using the rules of Hardy-Weinberg equilibrium we would need to sequence 1,111 horses to identify just one homozygote, therefore using this dataset alone, we cannot determine if this variant is embryonic lethal. Given that 96% of the variants with no homozygotes have an allele frequency around 0.50, it is highly unlikely that these variants are embryonic lethal; rather it is possible that these regions are related to mapping errors due to the presence of paralogs or pseudogenes, or the presence of structural variants. This is supported by the fact that all but two of the genes containing these variants were uncharacterized or equine only transcripts or olfactory receptor genes.

Identifying regions of the genome with high or low variation in the general population is critically important for the investigation of possible disease-causing variants, as regions with high genetic variability are less likely to contain disease-causing variants for fully penetrant Mendelian diseases (Karczewski et al., 2019). However, our analysis unexpectantly found that genes in high variability regions had lower haploinsufficiency scores, suggesting that damaging variants are less well tolerated (Huang et al., 2010) than for genes present in the low variability regions. This is likely due to the large number of genes without haploinsufficiency scores (70% of high variation region genes and 23% of low variation region genes). 65% of genes in high and 18% of genes in low variation regions had the “LOC” designation that are genes unique to the horse or are only predicted to be equivalent to a human gene, and therefore would not be included in human databases of variant constraint. This would suggest that as expected, genes in the low variability regions are more similar to human genes than genes in the high variability regions.

Overall, this is the first large-scale catalog of genetic variation in the domestic horse and will be highly useful for evaluation of background genetic variation in any future genetic study. This catalog has paved the way for future investigation of the regions of the genome that are shared or have marked MAF discrepancies across breeds. Regions with marked discrepancies between breeds, or between a breed and the population, can then be interrogated to further look for signatures of selection both across breeds, as well as within breeds. Additionally, further investigation into which of the variants that are present in all individuals are related to poor genome annotation or instead are true variants in the reference genome is needed and will lead to improvement of the horse reference genome in the future. By improving knowledge of the poorly annotated regions of the horse genome we will be able to correct these for future versions of the reference genome. This will have benefits for future phenotype-causing variant identification studies. Given the relatively unique utility of the horse as a model for human athletic related traits (Hill et al., 2010; Rooney et al., 2018) and diseases (Ward et al., 2004; McCue et al., 2008; McIlwraith et al., 2012; McCoy et al., 2013; Norton et al., 2016), an improved ability to identify phenotype-causing variants in the horse may shed light on analogous genetic diseases in humans. This is especially true for complex phenotypes such as athleticism, osteoarthritis (McIlwraith et al., 2012), and exertional rhabdomyolysis (Norton et al., 2016) where the limited genetic diversity in the horse (Petersen et al., 2013) may further accelerate our ability to identify the true phenotype-causing variants. This improvement in the reference genome combined with an improved understanding of the background genetic variation (Yang et al., 2014; Farwell et al., 2015; Ellingford et al., 2016; Hartmannová et al., 2016; Känsäkoski et al., 2016; Kojima et al., 2016; Noll et al., 2016; Smedley et al., 2016; Dolzhenko et al., 2017; Eldomery et al., 2017; Schneider et al., 2017) in the horse should vastly increase the identification of phenotype-causing variants for important equine diseases.