Origin and Evolution of Deleterious Mutations in Horses

Domestication has changed the natural evolutionary trajectory of horses by favoring the reproduction of a limited number of animals showing traits of interest. Reduced breeding stocks hampered the elimination of deleterious variants by means of negative selection, ultimately inflating mutational loads. However, ancient genomics revealed that mutational loads remained steady during most of the domestication history until a sudden burst took place some 250 years ago. To identify the factors underlying this trajectory, we gather an extensive dataset consisting of 175 modern and 153 ancient genomes previously published, and carry out the most comprehensive characterization of deleterious mutations in horses. We confirm that deleterious variants segregated at low frequencies during the last 3500 years, and only spread and incremented their occurrence in the homozygous state during modern times, owing to inbreeding. This independently happened in multiple breeds, following both the development of closed studs and purebred lines, and the deprecation of horsepower in the 20th century, which brought many draft breeds close to extinction. Our work illustrates the paradoxical effect of some conservation and improvement programs, which reduced the overall genomic fitness and viability.


Introduction
The domestication of the horse deeply impacted human history, enhancing the mobility of people, trade, and culture. For example, the diffusion of Indo-European languages has been associated with migration waves of horseback riders [1,2]. Following the incorporation of chariotry and cavalry into warfare, domestic horses have played a crucial role in the rise and fall of entire past civilizations. It was only with the onset of motor vehicles in the early 20th century that the horse remained consigned to farming and transportation in developing countries, and to recreation in Westernized societies. The equine industry remains instrumental today, with a census population size of 58.5 million horses, and a yearly market worth 300 billion US dollars [3].
Breeders have reshaped the horse evolutionary trajectory by controlling its reproduction for hundreds of generations [4]. Artificial selection has engendered several hundred breeds, showing striking phenotypic differences in a range of traits, including morphology, coat coloration, working capacities, and speed. However, developing such a broad array of breed-associated phenotypes has entailed profound genomic changes as demonstrated by Fages et al., who recently generated an extensive genome time-series spanning the last five millennia [5]. The authors found heterozygosity levels remaining relatively steady until they dropped by~16% some 200 years ago. This clearly

Materials and Methods
Sequencing data from modern horses were retrieved from ENA and NCBI repositories. Raw reads were processed and mapped against EquCab2 using PALEOMIX with default parameters [16]. Horse IDs, breed assignments, accession codes, and the resulting average depth-of-coverage per genome are reported in Table S1.
Evolutionary conservation is known to represent an excellent predictor of fitness effects, because mutations at sites that remained highly constrained during evolution are likely to be deleterious. We thus used phyloP conservation scores to estimate the potential impact of mutations. These scores summarize the evolutionary constraint position-wise, across a genome-wide alignment of 46 vertebrate species [17,18], including EquCab2 as reference assembly for the horse genome [19]. We considered as harmful all mutations at positions showing a minimum phyloP score of 1.5, which is a threshold that accurately discriminates fourfold and zerofold degenerate sites [20], i.e. synonymous (nearly-neutral) and non-synonymous (functional) variants, respectively. To calculate mutational loads per individual genome, we first estimated the genotype probabilities at each site, using ANGSD v0.917 [21] with the GATK likelihood model (- GL 2), and the following filtering parameters: -Uniqueonly 1 -remove_bads 1 -C 50 -baq 1 -minMapQ 30 -minQ 30 -only_proper_pairs 1. The most likely genotype was called and considered further, provided that its likelihood exceeded 0.99 and that it was homozygous (the phenotypic effect of deleterious variants at heterozygous sites depends on their unknown recessive, dominant, or co-dominant mode of inheritance). Sites were masked otherwise. Given both phyloP scores and genotype calls, we estimated the genetic load for each horse genome as: where i iterates over all the homozygous positions carrying a deleterious allele, and phyloP i is the phyloP score at the genomic position i. We assigned as deleterious the less frequent (or absent) allele in the 46-way alignment, provided that at most two variants were segregated.
Two approaches were applied to characterize the historical periods yielding inflated mutational loads. First, we retrieved all the previously published and radiocarbon-dated ancient genomes from ENA (PRJEB31613). After excluding mules, donkeys, and specimens belonging to other non-caballine lineages (i.e., not ancestral to modern horse breeds), we retained genomic data from 153 ancient horses (Table S2). These mostly lived during the last 3500 years. Then, we binned them within temporal windows of 1000 years, sliding every 250 years. For each time period, we estimated the frequency f of each deleterious allele from read count data using Maximum Likelihood (ML), where: where n_ind_t is the number of samples within each time interval, d i is the sequencing depth at that given position, and r di is the number of reads supporting the deleterious allele in individual i. We only considered those time bins showing at least 10 ancient horses genotyped, to minimize the variance associated with the estimation of f. By applying the same approach to the joint panel of 161 present-day domesticates (Table S1), we finally reconstructed the full temporal trajectory of each individual deleterious variant. Analyses were repeated considering transversions only, to limit the impact of post-mortem DNA misincorporations [22].
As a second approach to identify time periods of increasing loads, we calculated pairwise genetic distances with plink v2 [23] from a matrix that included the 161 modern domesticates, 14 Przewalski horses, and a donkey. We conditioned on 1,839,707,069 neutral positions (phyloP < 1.5), with one missing genotype at most. From these pairwise distances, we constructed a neighbor-joining (NJ) tree with subsequent topology refinement (-n option) [24]. Since low sequencing depths distort phylogenetic distances, all the genomes were pseudo-haploidized following a standard procedure in ancient DNA research. This consists in the random selection of one high-quality read at a given site as representative for the homozygous genotype. The tree branch lengths were used as proxies for neutral substitution rates, potentially revealing past episodes of elevated drift such as demographic collapses.
Pseudo-haploidized data were also used to characterize genetic purging, which involves selection against recessive mutations that are phenotypically exposed, such as those found at the homozygous state, both within and without Runs of Homozygosity (ROHs) resulting from inbreeding [25]. We approximated the strength selection by calculating the average genetic divergence of one given horse individual to the domestic donkey at constrained (dN; phyloP > 1.5) and neutral sites (dS; phyloP < 1.5). Strong negative selection purges out mutations at constrained sites, reducing dN and leading to negative dN-dS values. Conversely, deleterious variants are not efficiently removed under relaxed negative selection. Thus, dN behaves more neutrally, and effectively approaches dS, leading to dN-dS values closer to 0.
Finally, we estimated inbreeding coefficients, proceeding independently for Przewalski horses and modern domesticates. Specifically and for each group, we calculated the genotype posterior probabilities with ANGSD [21], and retained sites showing at most 10% missingness, provided that they segregated in approximate linkage equilibrium. The latter condition was satisfied through LD pruning and the calculation of r 2 [26] for Single Nucleotide Polymorphism (SNP) pairs located less than 50 Kb away in ngsLD [27]. We next clustered these SNP pairs into larger groups of linked variants using mcl [28], and selected the most central SNP as representative of each block. This yielded a total of 1,249,153 and 6,244,327 high-quality SNPs for Przewalski horses and modern domesticates, respectively. Inbreeding coefficients and IBD tracts were then co-estimated on these sites applying ngsF-HMM with a strict convergence criterion (min-epsilon = 1e −7 ) [29].

Levels and Patterns of Mutational Loads, Across Site Categories and Breed Types
The Przewalski horse represents an excellent starting model for understanding the biological significance of mutational loads, owing to the population collapse experienced the 20th century, which led to their extinction in the wild in 1969. The now~2100 animals living on the planet descend from a foundational captive stock of only 12-15 animals [30]. We first estimated individual loads within protein-coding regions, for comparison with previous work [5,20,31,32]. Averaging over 13 Przewalski horse genomes, and one Przewalski × Domestic F1 hybrid, we identified an average number of 1703.43 deleterious mutations out of 6,201,743 protein-coding sites. This corresponded to a mean load of approximately 3.698 × 10 −4 . As expected, most of the investigated breeds showed lower protein-coding loads than Przewalski horses, except for Shetland and Welsh ponies, as well as Marwari, Noriker, and Akhal Teke horses. Many of these breeds are presently represented by only one or two horse genomes; hence, we caution that the full range of possible load values present in these breeds remains to be explored. Other breeds such as Haflingers show highly variable loads, with some horses reaching values similar to those found in Przewalski horses. It is noteworthy that Haflingers and Norikers are draft breeds that are traditionally used as farm and pack animals. While only five major sire lines are described for Norikers, all modern purebred Haflingers can trace their ancestry back to one sire, Folie 249. Outbreeding was strictly prohibited in both breeds until recently, which limited founding stocks over multiple generations, leading to high load values (4.292 × 10 −4 and 4.294 × 10 −4 , respectively; Figure 1A). Shetland ponies also ranked high, which was probably due to long periods of isolation and genetic drift in a small British island, and to the selective crossing policy developed since the creation of their breeding society in 1890 [33].
As they represent only a limited fraction of the genome, protein-coding regions might provide partial and potentially biased estimates for the mutational load. Thus, we expanded the calculation to also cover non-coding regions, including the 2 Kb flanking gene bodies, introns and intergenic regions. This increased the number of positions considered~14-fold, representing approximately 67.6 million homozygous sites per genome. This also helped recover slightly and mildly constrained sites, which remained under-represented in protein-coding regions. This was so because protein-coding sites show increased evolutionary constraint relative to other regions, even at positions with phyloP scores greater than 1.5 (average phyloP scores, protein-coding = 2.316, 2 Kb upstream = 2.096, 2 Kb downstream = 2.097, intergenic = 1.964).
In general, the load estimates showed moderate correlation between the different regions considered (Spearman correlation; ρ < 0.586; p-value < 0.019), except for the protein-coding and 2 Kb upstream regions, where the correlation was non-significant (p-value = 0.949). We also found substantial differences on absolute scales, with loads within non-coding regions one order of magnitude greater (intergenic = 4.405 × 10 −3 , 2 Kb upstream = 4.034 × 10 −3 and 2 Kb downstream = 3.972 × 10 −3 ) than in protein-coding regions (3.476 × 10 −4 ). Considering that most of the horse genome is non-coding, and that selection seems to be more efficient in purging strongly deleterious mutations within gene bodies, we concluded that, in horses, loads predominantly accumulate at non-coding sites, through multiple mutations of small fitness effect. We next estimated genome-wide mutational loads, aiming at obtaining a fully representative set of positions, and potentially providing finer resolution for assessing the genomic consequences of different breed management practices. Results indeed revealed two major groups of breeds, each reflecting a major determinant of the current mutational burdens.
The top half of the load distribution is clearly dominated by traditional working breeds, including coldblood draft horses, as well as other farm and pack breeds ( Figure 1B). We suggest that this owes to most working breeds being abandoned since the mechanization of locomotion. Their recent population collapse likely limited the efficacy of negative selection and inflated loads. Such is the case of South Korean Jeju horses, which collapsed after industrialization, and accumulated an excess of deleterious mutations (load = 4.301 × 10 −3 ), and other draft and farm horses, such Lipizzans (4.297 × 10 −3 ), Haflingers (4.292 × 10 −3 ), and Norikers (4.294 × 10 −3 ). Shetland ponies also show high genomic loads (4.310 × 10 −3 ), which are even larger than those of the closely related Icelandic horses (4.294 × 10 −3 ) ( Figure 2). Both consisted originally of draft and farm animals, but have now been reconverted for leisure activities, and their census population size is limited. Likewise, the Marwari horse was endangered during the first half of the 20th century, until a series of conservation initiatives were started. The only Marwari representative analyzed in this study showed a genomic load of 4.312 × 10 −3 . This estimate was greater than that of Sorraia horses (4.283 × 10 −3 ), which is a breed once thought to be extinct, until a relict population was discovered and recovered, albeit incorporating some farm specimens of uncertain genetic backgrounds [34]. It is noteworthy that the genomic burden was particularly pronounced for Friesian horses (4.428 × 10 −3 ), the only breed exceeding the Przewalski mutational load genome-wide (4.310 × 10 −3 ). However, two Friesian horses (SAMEA3951222 and SAMEA3951223) had much lower loads than the other three breed members. These two genomes were found to be more homozygous (37.1 versus 46.3 millions homozygous sites), despite being sequenced at comparable depths-of-coverage (~8-9 ×, Table S1). The bimodal pattern found for Friesian horses could be compatible with genetic purging [25], as further investigated and discussed below.
The bottom half of the full-genome load distribution is enriched in breeds that were originally engendered for sport and leisure. This mostly consisted of hotblood and warmblood lines, which show comparable loads despite being subject to different breeding strategies. On the one hand, hotblood lines such as Arabian (4.245 × 10 −3 ) and Akhal-Teke (4.249 × 10 −3 ) horses trace their origins deep in the past, possibly hundreds of years before the raise of modern breeding practices. The only exception pertains, precisely, to the more recently founded Thoroughbreds, for which studbook registration started in 1791, and loads are inflated (4.268 × 10 −3 ). On the other hand, warmblood lines are more recent, but follow open stud guidelines that tolerate introgression from exogenous alleles, hence minimizing the deleterious effects of inbreeding (4.211-4.277 × 10 −3 , Figure 1B). Trakehners are worth a special mention, because unlike most warmblood horses, they are managed through nearly closed studbook practices, and expectedly show more elevated loads (4.281 × 10 −3 ). The benefit of admixture is also evident in a series of breeds that incorporated ancient Arabian lines, such the Connemara, Welsh, Miniature, and Reit ponies, as well as the Percheron horse [35], which exhibit only moderated loads ( Figure 1B). Finally, the Yakutian horse had the lowest loads ( Figure 1B), with an average of 301,007 harmful alleles in the homozygous state, corresponding to a mean genomic load of 4.171 × 10 −3 . Their lowest genomic loads are in line with the incorporation of multiple reproductive stallions within the Yakutian genetic pool, as illustrated by their high Y-chromosomal diversity [36].
However, these groups showed no difference in their protein-coding loads ( Figure 1A; Wilcoxon test; p-value > 0.1428). In addition to relying on a more limited number of positions, protein-coding loads represent regions evolving under stronger functional constraint, as reflected by their greater However, these groups showed no difference in their protein-coding loads ( Figure 1A; Wilcoxon test; p-value > 0.1428). In addition to relying on a more limited number of positions, protein-coding loads represent regions evolving under stronger functional constraint, as reflected by their greater average phyloP scores. Computer simulations conducted by Fages et al. clearly indicated that, after a population collapse, load bursts are almost undetectable from strongly constrained sites as selection remains sufficiently effective, but can be detected at slightly deleterious variation ( Figure S2 in [5]). Given that protein-coding loads seem far less sensitive to population collapses, it is thus not surprising that they fail to recover significant differences caused by the recent history of working, leisure, and hotblood lines.

Phylogenetic Reconstruction Supports Recent Population Decays in Working Breeds
To investigate population collapses, we built an NJ tree, which recovered strong bootstrap support for known phylogenetic relationships (Figure 2). In particular, Przewalski horses formed a sister group to all modern domesticates, with the F1 hybrid occupying the most basal position in this clade ( Figure 2). Within domesticates, Mongolian, Jeju, and Yakutian horses split first, followed by a clade of Icelandic, Miniature, and Shetland ponies. Hotblooded Akhal-Tekke and Arabian horses clustered jointly, with a Reitpony specimen, which is a breed that is known to have been influenced by hotblood lines. However, Thoroughbreds formed their own cluster that was well separated from other hotblood lines. Coldblooded draft horses were also monophyletic, including Percherons, Friesians, Norikers, and Haflingers. Finally, warmblood horses were grouped per breed, but showed a more complex pattern of diversification, reflecting their more admixed nature and introgression from influential breeds, which were either cold or hotblooded.
Genes 2019, 10, x FOR PEER REVIEW 7 of 16 average phyloP scores. Computer simulations conducted by Fages et al. clearly indicated that, after a population collapse, load bursts are almost undetectable from strongly constrained sites as selection remains sufficiently effective, but can be detected at slightly deleterious variation ( Figure S2 in [5]). Given that protein-coding loads seem far less sensitive to population collapses, it is thus not surprising that they fail to recover significant differences caused by the recent history of working, leisure, and hotblood lines.

Phylogenetic Reconstruction Supports Recent Population Decays in Working Breeds
To investigate population collapses, we built an NJ tree, which recovered strong bootstrap support for known phylogenetic relationships (Figure 2). In particular, Przewalski horses formed a sister group to all modern domesticates, with the F1 hybrid occupying the most basal position in this clade (Figure 2). Within domesticates, Mongolian, Jeju, and Yakutian horses split first, followed by a clade of Icelandic, Miniature, and Shetland ponies. Hotblooded Akhal-Tekke and Arabian horses clustered jointly, with a Reitpony specimen, which is a breed that is known to have been influenced by hotblood lines. However, Thoroughbreds formed their own cluster that was well separated from other hotblood lines. Coldblooded draft horses were also monophyletic, including Percherons, Friesians, Norikers, and Haflingers. Finally, warmblood horses were grouped per breed, but showed a more complex pattern of diversification, reflecting their more admixed nature and introgression from influential breeds, which were either cold or hotblooded. Neighbor-joining (NJ) cladogram depicting the horse phylogenetic relationships, mid-point rooted using the donkey as the outgroup. Color ranges highlight monophyletic lineages, including Przewalski horses (green), Mongolian-derived (brown) and Nordic (orange) breeds, as well as coldblooded (blue) and hotblooded (red) lines. Most of the non-colored taxa correspond to admixed warmblood lines, whose phylogenetic placement depends on the relative contribution of other influential breeds. Long internal branches are highlighted in red, as possibly reflecting episodes of elevated drift. Only bootstrap support values lower than 95% are displayed to avoid overloading the tree.
We further scrutinized the length of internal branch lengths to potentially reveal past episodes of increased genetic drift. The longest internal branches led to Sorraia horses (7 × 10 −5 substitutions per nearly-neutral site), Przewalski horses (6 × 10 −5 ), Haflingers (5 × 10 −5 ), Friesians (5 × 10 −5 ), and Lipizzans (2.6 × 10 −5 ) (Figure 2). The foundational branches of Icelandic (1.7 × 10 −5 ), Shetland (1.7 × 10 −5 ), and Jeju (1.5 × 10 −5 ) ponies were slightly shorter, on par with those leading to Akhal-Teke (1.9 × 10 −5 ) and Arabian (1.2 × 10 −5 ) horses. These long internal branches echoed the mutational loads carried by working horses, suggesting independent demographic bottlenecks reducing the efficacy of negative selection and inflating their mutational loads. It is noteworthy that the branch length leading to all domesticates was only 1.7 × 10 −5 substitutions per nearly-neutral site, despite encompassing the~45,000 years of divergence with Przewalski horses [5]. This suggests that the genomic signature left by the domestication bottleneck was mild, and relative to that observed in some modern breeds. This mild bottleneck is consistent with the large mitochondrial diversity found in horses, which was interpreted as a pervasive restocking of wild mares during the initial spread of horse husbandry [37,38].

Deleterious Mutations Segregated at Low Frequencies Until the Last~250 Years
We next aimed at reconstructing the past historical dynamics leading to present-day mutational loads. To achieve this, we first exploited genome-scale data from 153 ancient domestic horses and tracked the trajectories of all the deleterious alleles segregating in modern breeds over the last 3500 years (n = 1,313,308). Figure 3A-C provide illustrative examples of the temporal trajectories of a selection of alleles known to be associated with diseases [39][40][41], including increasing and decreasing trends as well as cases where variation does not follow simple temporal changes. Overall, we detected that~11.3% of the deleterious mutations were nearly fixed over time across all the periods, including in present-day domesticates (ML frequency ≥ 0.99). Thus, these harmful mutations spread prior to 3500 years ago, and probably prior to domestication. However, the vast majority of deleterious variants (~76.6%) remained nearly absent at all time periods (ML frequency < 0.01). This proportion increased to~84.2% and~86.9% when conditioning on more constrained positions (phyloP scores ≥ 2 and 2.5, respectively). This suggests that negative selection successfully maintained most of the deleterious variants at low frequencies during recent horse evolution.
We observed that the remaining fraction of deleterious variants (~12.0%, n = 158,448) followed a dynamic temporal trajectory, which is defined as detectable changes in frequency across successive time periods (∆). To further quantify ∆, we conditioned on non-overlapping time bins of 1000 years, centered at 250, 1250, 2250, and 3250 years ago. This was done to ensure no redundancy across adjacent intervals, and hence to avoid underestimating ∆, since overlapping windows comprised of almost the same horses would provide nearly identical allele frequencies. We found that the most recent time interval tested, representing the last 250 years, experienced the largest shift in allele frequency ( Figure 3D). Its median change over the 158,448 deleterious alleles was ∆ = 0.02399, while it was 0.02048 or less for older time periods (Wilcoxon test; p-value < 2.2 −16 ). These changes entailed increases or decreases in frequencies at equal proportions, except for the last 250 years, where 71% of deleterious variants became more common than in the previous time interval. This holds true when disregarding transitions, suggesting that the temporal trend was robust to the possible presence of post-mortem DNA misincorporations in the sequence data (although these were likely limited due to the treatment of most ancient DNA extracts with USER prior to DNA library preparation; Figure 3D).
As the results above supported those by Fages et al. [5], in which mutational loads were steady during millennia prior to the industrial revolution, we jointly considered all ancient specimens. This provided an approximately even number of 161 modern and 153 ancient horses to confidently estimate the frequencies of all the deleterious alleles (n = 1,313,308). On average, we detected that deleterious mutations are more common in modern horses, compared to their ancient relatives, by 0.3% according to transitions and by 0.8% to transversions ( Figure 3E). Note that transitions are more spread than transversions also in modern horses, suggesting that their greater frequency is not due to post-mortem damage, but to sequencing biases and/or biological processes, such as CpG hypermutability and selection against transversions [42,43]. Taken together, these findings confirm that current deleterious mutations segregated in the past, but that it was only recently that they significantly rose in frequency.
As the results above supported those by Fages et al. [5], in which mutational loads were steady during millennia prior to the industrial revolution, we jointly considered all ancient specimens. This provided an approximately even number of 161 modern and 153 ancient horses to confidently estimate the frequencies of all the deleterious alleles (n = 1,313,308). On average, we detected that deleterious mutations are more common in modern horses, compared to their ancient relatives, by 0.3% according to transitions and by 0.8% to transversions ( Figure 3E). Note that transitions are more spread than transversions also in modern horses, suggesting that their greater frequency is not due to post-mortem damage, but to sequencing biases and/or biological processes, such as CpG hypermutability and selection against transversions [42,43]. Taken together, these findings confirm that current deleterious mutations segregated in the past, but that it was only recently that they significantly rose in frequency.

Inbreeding Depression and Genetic Purging Shaped Loads in Modern Domesticates
Understanding the evolutionary mechanisms that forged current mutational loads is paramount to identifying (un)desirable breeding practices and designing more sustainable strategies. Given the recent time-scale delimited within the last 250 years, inbreeding depression represents the most likely mechanism underlying the recent increase of mutational loads in the horse genome [44,45]. Inbreeding depression is caused by recessive mutations that are phenotypically hidden in the heterozygous state, until they become effective once located within the runs of homozygosity (ROHs)

Inbreeding Depression and Genetic Purging Shaped Loads in Modern Domesticates
Understanding the evolutionary mechanisms that forged current mutational loads is paramount to identifying (un)desirable breeding practices and designing more sustainable strategies. Given the recent time-scale delimited within the last 250 years, inbreeding depression represents the most likely mechanism underlying the recent increase of mutational loads in the horse genome [44,45]. Inbreeding depression is caused by recessive mutations that are phenotypically hidden in the heterozygous state, until they become effective once located within the runs of homozygosity (ROHs) introduced by inbreeding. Assuming that inbreeding depression was the main driver of the mutational load patterns observed in this study, we should expect that: (1) recessive mutations segregated in the ancestral population, and (2) negative selection was not sufficiently strong to remove recessive mutations exposed within ROHs [25].
The first assertion was proved in the section above. In order to test the second, we characterized inbreeding and identified ROHs in each individual modern horse genome using ngsF-HMM [29].
Our inbreeding estimates replicated previous work for the 14 Przewalski horses [46], showing an average inbreeding coefficient of F = 18.5%, corresponding to~18.5% of the genome being identical-by-descendant (IBD) ( Figure 4A). The Przewalski x Domestic F1 hybrid analyzed (KB7903) showed no inbreeding, and also had the lowest mutational load in the group ( Figure 1B). This suggests that the inbreeding estimates that were recovered are genuine. We found that inbreeding coefficients and mutational loads were strongly correlated in Przewalski horses (Spearman correlation; ρ = 0.903; p-value = 9.740 × 10 −6 ), which was as expected if negative selection was not sufficiently strong to eliminate the deleterious variants exposed in ROHs.

Discussion
Recent work from Fages et al. revealed that the last few centuries have been accompanied by a ~16% drop in the horse heterozygosity genome-wide, and a ~4% raise in the mutational load within protein-coding regions [5]. However, the underlying drivers of these shifts remained unclear. To address this gap, we carried out an extensive characterization of mutational loads in horses, Modern domesticates returned slightly lower inbreeding coefficients (on average, F = 15.9%). This estimate is greater than the 8.8% recently estimated across nine breeds, based on 65,157 SNPs only [47], suggesting strong ascertainment bias within this SNP set. Ranking per breed revealed that Shetland, Sorraia, and Thoroughbreds were extremely inbred, with values approaching and even exceeding F = 30% ( Figure 4A). Their longest IBD tracts spanned 20 Mb, 10 Mb, and 11 Mb, respectively ( Figure 4B). Interestingly, and in contrast to Przewalski horses, inbreeding did not necessarily entail increased loads, as the correlation was weaker, albeit significant (ρ = 0.222; p-value = 4.708 × 10 −3 ). Limiting the calculations to those 53 domesticates sequenced above 15 × strengthened the correlation coefficient (ρ = 0.4754); however, it remained inferior to that inferred for Przewalski horses. A similar trend was found conditioning load estimates on protein-coding sites. This suggests that additional mechanisms, beyond inbreeding depression, have contributed to shape the mutational load present in modern breeds.
As genetic purging involves strong selection against recessive mutations exposed within ROHs, we propose that this mechanism could have reduced the mutational loads in the most inbred horses. This mechanism may have been inefficient in Przewalski horses due to the extremely limited reproductive stock available for the conservation program and/or to the favorable environmental conditions offered in captivity and reintroduction reserves. To further validate whether selection was stronger in modern domesticates than in Przewalski horses, we quantified the difference between non-synonymous and synonymous mutation rates, dN-dS, in each individual genome (see methods). All the horses had negative dN-dS values, as expected under negative selection ( Figure 4C). Yet, Przewalski horses appeared at the higher end of the dN-dS distribution (mean dN-dS = −3.328 × 10 −3 ), confirming more relaxed negative selection in this lineage relative to domestic horses. Amongst modern domesticates, the highly isolated Shetland ponies represent the only breed showing lower negative selection than Przewalski horses (−3.306 × 10 −3 ). Interestingly, Thoroughbreds were found to be at the tail of the dN-dS distribution (−3.366 × 10 −3 ). The correlation between load and dN-dS, which was calculated across the 19 Thoroughbreds investigated, was non-significant (ρ = 0.178; p-value = 0.467), supporting ongoing genetic purging in this breed.

Discussion
Recent work from Fages et al. revealed that the last few centuries have been accompanied by a~16% drop in the horse heterozygosity genome-wide, and a~4% raise in the mutational load within protein-coding regions [5]. However, the underlying drivers of these shifts remained unclear. To address this gap, we carried out an extensive characterization of mutational loads in horses, leveraging previously published genome data from a total of 175 modern horses spanning 37 breeds and/or populations, and 153 ancient horses. We expanded the calculation of mutational loads outside protein-coding regions, which enhanced both resolution and accuracy.
Our findings support inbreeding depression as the main mechanism driving the load burst in domesticates and Przewalski horses. It is important to keep in mind that the last~250 years cannot generate sufficient amounts of de novo variants to explain the observed increment in load, given the low mutation rate inferred for horses [48]. Therefore, the excess of mutational load was almost entirely driven by standing variation; it was also likely located in non-coding regions, and associated with slightly deleterious and recessive inheritance (i.e., dominant deleterious mutations are less frequent, because they are phenotypically exposed to negative selection at the heterozygous state, and thus efficiently eliminated). We indeed confirm that deleterious variation segregated in the past, at very low frequencies, until rising recently. Nevertheless, their increments of only~0.3% (transitions) and 0.8% (transversions) is lower than the 4% raise in protein-coding loads [5]. Hence, the load burst cannot be only explained by a higher frequency of deleterious mutations, but requires that they increasingly became exposed at the homozygous state, owing to inbreeding. The significant correlation between inbreeding and the mutational load identified here further corroborates the role of inbreeding in causing a fitness depression.
Inbreeding was caused by two main historical shifts. The ubiquity of steam and combustion vehicles relegated breeds traditionally used for farming and transport to almost oblivion, resulting in fast population collapses (and even extinction in some cases) [49]. Although inbreeding exposes recessive deleterious variants to negative selection within ROHs, the reduced effective sizes considerably limited the efficacy of negative selection. For example, this increased mutational loads in Sorraia, Haflingers, Norikers, and especially Friesian horses, which show larger mutational loads than the endangered Przewalski horses ( Figure 1B). One exception to this pattern pertains to one single Percheron horse present in our genome dataset. This breed consisted of large, draft horses, and is known to have been extremely influenced by Arabian bloodlines while being founded in France, before becoming extremely popular worldwide, especially in the U.S., in the 19th century [50]. As a result, its geographic range was exceptionally large for a draft breed, which may have helped preserve sufficient genetic diversity, until conservation programs started in the 1960s.
Conservation programs often rely on closed stud practices to maintain un-admixed populations that are better adapted to native environments. For endangered breeds, this implies that foals can only be registered if descending from purebred studbook-registered parents. Similar rules govern selection programs to improve specific breeds, such as Thoroughbreds, for which a studbook was established well before Darwin formalized the concept of evolution through natural selection in 1859 [51]. Closed studs, and the preferential reproduction of influential stallions, increased inbreeding and the probability of exposing deleterious alleles in the homozygous state. In extremely large studs such as Thoroughbreds, which are intensively selected for performance, this may have provided sufficient strength to purge deleterious mutations. However, in other breeds, which are either less intensively selected or restricted to extremely low effective sizes, mutational loads could spread.
The abandonment of working breeds and the emergence of closed studbooks clearly post-date the onset of domestication by at least five millennia [32,52,53]. This has deep implications for the cost-of-domestication hypothesis [54], which posits that there was an increase in mutational load due to the repeated bottlenecks presumably experienced during the early domestication stages [31,[55][56][57]. In agreement with previous work in horses [20] and crops [58], we find that this is not necessarily the case, and highlight that two forces with opposing effects could have also contributed to shape mutational loads. On the one hand, the impact of recent events seems unprecedented in the horse evolutionary history, and appears to have eroded the horse genome more than the domestication bottleneck itself (Figures 2 and 3). On the other hand, restocking from the wild and cross-breeding during hundreds of generations could have counteracted the deleterious consequences of early population declines, e.g., through heterosis, which is a phenomenon involving greater vigor and fertility in hybrids than in their parental inbreed stocks [13].
A number of mathematical models for inbreeding depression, heterosis, and genetic purging [14,25,59,60] have predicted that populations with reproductive stocks that are comparable to what are found in many domestic breeds undergo a strong risk of extinction. For example, Caballero et al. recently used simulations to estimate that populations limited to approximately 70 reproductive individuals are under substantial risk of extinction, as defined by a >10% reduction in viability after 50 generations of evolution [61]. Note that a breeding stock of N m = 20 stallions and N f = 5000 mares corresponds to a population size of only N e = 4N m N f /(N m +N f ) ≈ 80 individuals. In line with this, and according to the Domestic Animal Diversity Information System (DAD-IS; accessed 19 July 2019 from the FAO website [49]), more than 200 horse breeds would thus be endangered or at the brink of extinction, while 88 are already extinct.
Encouragingly, the relation between inbreeding and mutational loads was found to be impacted by genetic purging in modern domesticates ( Figure 4). Thus, genetic purging adds an extra layer of complexity to the interplay of forces that forged genomic loads, not only helping to improve fitness, but also to optimize traits that are paramount to the equine industry. For example, in Thoroughbreds, genetic purging has been associated with improved racing performance [62]. This means that breeders can leverage genomic information to design mating strategies favoring the purging of deleterious mutations from the breeding stock, improving animal welfare and mitigating extinction risks.