Decline of genetic diversity in ancient domestic stallions in Europe.

Present-day domestic horses are immensely diverse in their maternally inherited mitochondrial DNA, yet they show very little variation on their paternally inherited Y chromosome. Although it has recently been shown that Y chromosomal diversity in domestic horses was higher at least until the Iron Age, when and why this diversity disappeared remain controversial questions. We genotyped 16 recently discovered Y chromosomal single-nucleotide polymorphisms in 96 ancient Eurasian stallions spanning the early domestication stages (Copper and Bronze Age) to the Middle Ages. Using this Y chromosomal time series, which covers nearly the entire history of horse domestication, we reveal how Y chromosomal diversity changed over time. Our results also show that the lack of multiple stallion lineages in the extant domestic population is caused by neither a founder effect nor random demographic effects but instead is the result of artificial selection-initially during the Iron Age by nomadic people from the Eurasian steppes and later during the Roman period. Moreover, the modern domestic haplotype probably derived from another, already advantageous, haplotype, most likely after the beginning of the domestication. In line with recent findings indicating that the Przewalski and domestic horse lineages remained connected by gene flow after they diverged about 45,000 years ago, we present evidence for Y chromosomal introgression of Przewalski horses into the gene pool of European domestic horses at least until medieval times.

INTRODUCTION

Present-day domestic horses are unique among livestock animals because of their almost complete lack of variability on the paternally inherited Y chromosome (–). In stark contrast to that, the maternally inherited mitochondrial DNA of horses shows huge variation (–). So far, only very few studies detected polymorphic Y chromosomal markers among modern horses, but these additional lineages seem to be derived from the ancestral domestic haplotype relatively recently () or were limited to a restricted geographical area ().

The pattern of high maternal diversity combined with extremely reduced paternal variation is generally explained by unequal sex ratios at the beginning of domestication. Likely, a large number of mares but only few stallions contributed to founding the domestic horse population (). However, further support for these discrepancies arose from continued sex-specific breeding enhancing the reproductive success of certain male lineages, particularly after the development of modern studs during the last few centuries (). Recent ancient DNA (aDNA) analyses of paternal variation have detected the coexistence of multiple Y chromosomal haplotypes in pre–Iron Age and Iron Age domestic stallions, which suggests that Y chromosomal diversity was much larger in past than in present-day horse populations (, ). However, because of limited sampling efforts, these earlier studies could not capture the diachronical changes in paternal diversity.

To clarify how the domestication process of horses shaped the diversity of the Y chromosome, we genotyped 16 loci on the Y chromosome in 96 (out of a total of 350) archaeological remains (table S1), spread across Eurasia and spanning the last ~5000 years. This data set provides the first full time series of Y chromosomal variation from early domestication stages until the Middle Ages.

RESULTS AND DISCUSSION

To clarify how the domestication process of horses and the potential input of the Przewalski horse shaped the diversity of the horse Y chromosome, we successfully genotyped 16 single-nucleotide polymorphism (SNP) loci on the Y chromosome in 96 (out of a total of 350) archaeological horse remains (table S1), spread across Eurasia and spanning the last ~5000 years. In agreement with previous studies, the Y chromosome showed low genetic diversity, with only 4 of 16 SNPs being polymorphic, defining four haplotypes (table S2). In addition to the haplotypes exhibited by present-day domestic horses (Y-HT-1), we detected two new haplotypes (Y-HT-3 and Y-HT-4) characterized by two differences and one difference from the domestic haplotype, respectively (Fig. 1). Somewhat surprisingly, we also detected in several of the successfully genotyped aDNA samples the Przewalski horse haplotype (Y-HT-2), which differs at all four polymorphic positions from the present-day domestic haplotype.

Fig. 1

Stallion lineages through time.

Temporal haplotype network of the four detected Y chromosome haplotypes. Age of the samples indicated by multiple layers separated by color; vertical lines connecting the haplotypes of consecutive layers/ages represent which haplotype was transferred into a later/younger period. Numbers constitute the respective number of individuals showing this particular haplotype for that period. Prz, Przewalski; Dom, domestic.

To detect patterns of temporal variation, we divided the samples into four different time bins based on archaeological periods (before 2200 BCE, 2200–900 BCE, 900 BCE–400 CE, and younger than 400 CE; Fig. 2). On the basis of these, we constructed a temporal haplotype network (Fig. 1) and simulated the haplotype frequencies over time (Fig. 3). All analyses show the reduction of Y chromosome variation through time, which is especially pronounced for haplotypes Y-HT-1 and Y-HT-4. Wild and early domestic horses from the Neolithic and Bronze Age (>900 BCE) display comparatively high variation, despite their smaller sample size (n = 9 and 17 in time bins 1 and 2 versus n = 24 and 46 in time bins 3 and 4; Fig. 2), with all four detected haplotypes occurring at relatively even frequency in the first two defined periods (~0.10 to 0.36; Fig. 3 and table S3). Moreover, among European horses, the haplotype representing modern Przewalski horses (Y-HT-2) was very abundant during the two older periods and even constitutes the most frequent one before 2200 BCE (oldest time bin). Haplotype Y-HT-1, on the other hand, which dominates in present-day stallions, was only detected later than 2200 BCE. Although its estimated age roughly correlates with the onset of domestication (3500 BCE; Fig. 4), this haplotype only started to become more frequent during later periods. All analyses indicate low Y-HT-1 frequency before 2200 BCE (fig. S1). However, the first time bin (>2200 BCE) only includes samples from Europe, with the easternmost samples falling on the longitude of the Black Sea (Fig. 2). Accordingly, Y-HT-1 could have been present with higher frequencies at that time in populations from further East. In the Iron Age (time bin 3, 900 BCE–400 CE), Y-HT-1 already represented the most frequent haplotype.

Fig. 2

Decline of paternal diversity began in Asia.

Maps displaying age, locality, and haplotype (different colors) of each successfully genotyped sample.

Fig. 3

The frequencies of Y chromosome haplotypes started to change during the Late Bronze Age (1600–900 BCE).

Inferred temporal trajectories of haplotype frequencies. Each haplotype is displayed by a different color. The shaded area represents the 95% highest-density region. The trajectories were constructed taking the median values across frequencies from the simulations of the Bayesian posterior sample. The small chart represents the stacked frequencies; the amplitude of each colored area is proportional to the median haplotype frequencies (normalized) at a given time. The x and y axes of the small chart match those in the large one. Ka, thousands of years.

Fig. 4

Positive selection for modern domestic haplotype leads to the decline of diversity.

Violin plots of selection coefficients and initial haplotype frequencies, representing the posterior density functions of the selection coefficients of haplotypes Y-HT-1 and Y-HT-4, the age of Y-HT-1, and the initial haplotype frequencies of all haplotypes but Y-HT-1. The scale on the left is for the selection coefficients, whereas the scale on the right is for the age of Y-HT-1. The initial allele frequencies are on a scale from 0.0 (bottom) to 1.0 (top).

The first record of the modern domestic Y chromosome haplotype stems from two Bronze Age samples of similar age. Notably, both samples were found in two distantly located regions: present-day Slovakia (2000–1600 BCE, dated by archaeological context) and western Siberia (14C-dated: 1609–1436 cal. BCE). Although a very recent study proposes an oriental origin of this haplotype (), we cannot determine the geographical origin of Y-HT-1 with certainty, because this haplotype has not been found thus far in predomestic or wild stallions. There are two possible scenarios: (i) Y-HT-1 emerged within the domestic population by mutation and (ii) Y-HT-1 was already present in wild horses and entered the domestic population either at the beginning of domestication (but initially restricted to Asian horses) or later by introgression (from wild Y-HT-1 carrying studs during the Iron Age). Crosses between domestic animals and their wild counterparts have been observed in several domestic species (–); thus, the simplest explanation would be that we missed Y-HT-1 in older samples because of limited geographical sampling. However, the estimated haplotype age is contemporary (Fig. 4) with the assumed starting point of horse domestication ~4000–3500 BCE (), rendering it likely that Y-HT-1 originated within the domestic horse gene pool. Still, we cannot rule out definitively that it appeared before domestication.

Independent of its geographical origin, Y-HT-1 progressively replaced all other haplotypes—except for one additional lineage that is restricted to Yakutian horses (). Considering our data, this trend in paternal diversity toward dominance of the modern lineage appears to start in the Bronze Age and becomes even more pronounced during the Iron Age. The Bronze Age was a time of large-scale human migrations across Eurasia (–), movements that were undoubtedly facilitated by the spread of horses as a means of transport and warfare. At that time, the western Eurasian steppes were inhabited by highly mobile cultures that largely relied on horses (, , , ). The genetic admixture of northern and central European humans with Caucasians/eastern Europeans did correlate with the spread of the Yamnaya culture from the Pontic-Caspian steppe (), an area that has repeatedly been suggested as the center of horse domestication (, , ). Given the importance of domestic horses, it appears that deliberate selection/rejection of certain stallions by these people might have contributed to the loss of paternal diversity. The spread of humans out of this region might also have resulted in the spread of Y-HT-1 from Asia to Europe. This scenario also agrees with recent findings that the low male diversity of extant horses is not caused by recruiting only a limited number of stallions during early domestication ().

To evaluate whether the observed changes in Y haplotype frequencies can be explained by drift or whether selection needs to be invoked, we applied several analyses. The different methods assumed different demographic scenarios. All methods consistently provided substantial statistical support for positive selection on Y-HT-1 (100% positive coefficients; Fig. 4 and fig. S3). Thus, it appears that deliberate breeding contributed to the loss of paternal diversity. This inference was also supported by our model choice analysis, which provided overwhelming support for the selection model (Bayes factor = 9.4 × 1011) when compared to a pure drift model.

When we tested whether a temporally changing selection coefficient provided a better fit for the data, by assigning different coefficients to different periods, we found that one single coefficient for all periods obtained the best combined support of the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) (fig. S4), suggesting that the selective forces acting upon haplotype Y-HT-1 did not change substantially over time.

In addition, we detected strong selection acting also on other Y chromosome haplotypes. Although values were highly positive for Y-HT-4 (98.2% of selection coefficients were positive), selection for Y-HT-2 and Y-HT-3 was significantly negative (100% of coefficients negative). Our analyses indicate that Y-HT-2 declined throughout the sampled time frame, Y-HT-3 initially remained constant and then declined, Y-HT-4 rose before 1500 BCE and subsequently decreased in frequency, and Y-HT-1 almost aggressively invaded the population until fixation (Fig. 3). The whole pattern, including the turnover of Y-HT-4, is best explained by positive selection on Y-HT-4 and even stronger positive selection on Y-HT-1 (Fig. 4).

Y-HT-1 rapidly replaced the already advantageous haplotype Y-HT-4, its most closely related lineage. With Y-HT-4 probably being older than Y-HT-1, it seems possible that the selective advantage of Y-HT-4 was retained in Y-HT-1, which then acquired its own additional advantageous mutations and, consequently, an even higher selective value. Such a situation, where younger haplotypes are under increasingly higher selection, could be a case of selective escalation. In this scenario, the Y-HT-1–specific mutations further improved the advantageous phenotype of Y-HT-4, allowing Y-HT-1 to invade the gene pool of domestic horses. The speed with which Y-HT-1 reached fixation supports this escalation hypothesis. Future research could focus on identifying genes that display nonsynonymous mutations shared only by Y-HT-1 and Y-HT-4 and investigating their phenotypic expression. Possible traits could well be phenotypically visible features because increasing evidence indicates that, in addition to spermatogenesis and male fertility, several diseases and physiological and behavioral traits in males are regulated by Y chromosomal genes (); such Y-chromosome-associated traits could explain how early breeders chose the stallions they selected for breeding.

During the Iron Age, when Y-HT-1 started to obtain its dominating status, the Roman Empire was ruling Europe and its bordering regions. The Romans changed the mode of animal breeding in their territories by shifting the focus from female-based to male-based selective breeding (). This change in animal husbandry should have raised the reproductive success of the few stallions selected for reproduction, increasing the frequency of their Y chromosomal lineages. Although we cannot exclude the possibility that animal husbandry and breeding lost their organization with the end of the Roman Empire, the fact that horses were basically fixed for Y-HT-1 already at the beginning of the Middle Ages would have prevented other haplotypes from rising again during medieval times.

The presence of the Y chromosome haplotype carried by present-day Przewalski horses (Y-HT-2) in early domestic stallions and a European wild horse (Pie05; table S2) could be the result of introgression of Przewalski stallions. Although the original distribution of the Przewalski horse is unknown, it was probably much larger than that of the relict population in Mongolia that produced modern Przewalski horses and might even have extended into Central Europe. However, it is also possible that either Przewalski horses were among the initially domesticated horses or that Y-HT-2 occurred both in Przewalski horses and in those wild horses that are the ancestors of domestic horses, based on autosomal DNA data (). Regardless of how Y-HT-2 entered the domestic gene pool, it was eventually lost, as were all haplotypes except Y-HT-1. In our sample set, Y-HT-2 was undetectable as early as the third time bin. However, it is possible that Y-HT-2 may have been present during this time period, but with a frequency below 0.11 (with 95% probability). The inferred time trajectories for Y-HT-2 frequencies suggest that it could nevertheless have persisted at very low frequencies until the Middle Ages (Fig. 3). On the basis of these simulations, this finding could be interpreted as a relic of this haplotype’s formerly higher frequency in the domestic horse gene pool. It is also possible that the presence of this haplotype could be the result of mating a wild stallion with a domestic mare, a frequently reported breeding practice when wild horses were still widely distributed. However, a significant contribution of the Przewalski horse to the gene pool of modern domestic horses has been almost ruled out by recent genomic studies (, , ).

MATERIALS AND METHODS

Sample preparation, aDNA extraction, amplification, and sequencing

We analyzed 350 samples from China; Mongolia; Siberia; eastern, middle, and northern Europe; and the Iberian Peninsula (table S1). The ages of the samples range from the Neolithic/Copper Age to medieval times and were determined either by 14C-dating of the samples or by their archaeological context. DNA was extracted from 100- to 200-mg samples of bone or tooth in laboratories dedicated to aDNA analyses at the Leibniz Institute for Zoo and Wildlife Research (Ludwig Laboratory) in Berlin, Germany, following the standard procedures to avoid contamination. Independent replications were carried out in a different laboratory at the University of Potsdam (Hofreiter Laboratory), Germany. The tissue surfaces of all samples were removed by abrasion to minimize contamination. Samples were ground with a cryomill and incubated in 1 ml of extraction buffer [0.45 M EDTA, proteinase K (0.25 mg/ml) (pH 8.0)] overnight under rotation at 37°C. An initial scan for the sex of the samples was done by polymerase chain reaction (PCR) using sex-specific primers [see the study of Lippold et al. () for details] and subsequent visualization on an agarose gel. The following steps were conducted according to an optimized extraction protocol for short DNA fragments (). All samples were analyzed for 16 previously published SNPs () amplified with 12 primer pairs (table S4). Multiplex PCRs were performed in 20-μl reactions using 4 μl of DNA extract and further containing 1× AmpliTaq Gold PCR buffer II (Applied Biosystems); 4 mM MgCl2; bovine serum albumin (1 mg/ml); 250 μM dATP (2′-deoxyadenosine 5′-triphosphate), dCTP (2′-deoxycytidine 5′-triphosphate), and dGTP (2′-deoxyguanosine 5′-triphosphate); 500 μM dUTP (deoxyuridine 5′-triphosphate); 150 nM of each primer; and 2 U of AmpliTaq Gold (Applied Biosystems). For each sample, one PCR was performed by adding 1 U of heat-labile uracil-DNA glycosylase and an initial incubation step of 15 min at 37°C to control for carryover contamination. Target regions were amplified with target-specific primer pairs that have been tagged with common sequence tags (CS1 and CS2, Fluidigm). After multiplex amplification, PCR products were sequenced using a combination of multiplexed Illumina MiSeq sequencing and traditional Sanger capillary electrophoresis. For Illumina sequencing, the amplicons were indexed using a maximum of 384 barcodes (containing Illumina adaptor sequences) in an additional amplification step and then pooled. To remove primer and adaptor dimers and other amplification products, we purified the multiplex PCR products using an AMPure PCR purification system (Agencourt) with a 1.8-fold ratio of SPRI (solid phase reversible immobilization) beads to reaction volume. The pooled and purified amplicon libraries were then sequenced on an Illumina MiSeq sequencer using CS1 and CS2 and their reverse complements as sequencing primers according to the manufacturer’s instructions. For a few samples, the PCR products were generated with the same target-specific primer pairs and then sequenced on the ABI 3130xl Genetic Analyzer by applying the BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific).

Processing and alignment of sequencing data

Traditional Sanger sequences were edited and aligned in Geneious v8.1 (Biomatters). Illumina reads were automatically demultiplexed by the Illumina MiSeq software and further processed by applying the following steps. Illumina sequencing adaptors were cut from the reads using cutadapt (), and reads were then filtered for quality with Trimmomatic 0.33 (), discarding all reads that contain more than 8 base pairs (bp) with a Phred score below 20 and allowing only reads with a minimum length of 50 bp. The trimmed and filtered paired-end reads were merged to a single read with FLASH (), improving the overall quality. Using Bowtie 2 () and SAMtools 0.1.19 (), these merged reads were mapped to the reference sequences (G72335.1, G72336.1, G72338.1, G72339.1, and AB091794.1) for the respective target regions. For each sample and locus, a consensus sequence was generated in Geneious v8.1 (Biomatters) from which the genotype of the Y chromosomal SNPs could be deduced. With these, SNP haplotype frequencies were calculated using the R package “pegas” v0.8-1 (). To visualize the detected haplotypes during the different time periods, we constructed a heterochronous parsimony network using the R-script TempNet (). We also tested for temporal changes in haplotype frequencies among four time bins (before 2200 BCE, 2200–900 BCE, 900 BCE–400 CE, and younger than 400 CE) by means of a temporal test of allele frequencies ().

Statistical analysis

We applied the sequential Markov Chain Monte Carlo (S-MCMC) method reported by Wutke et al. (, ) to estimate coefficients of selection that historically acted upon the Y chromosome haplotypes in horses. This method performs statistical inference by using computer simulations in a Bayesian MCMC framework. The use of simulations allowed taking into account the uncertainties associated to radiocarbon or stratigraphic ages of the ancient samples, sampling error, generation time, population sizes, and the change of generation, which was modeled by a Wright-Fisher model (, ). Following previous work (, , ), the simulations also incorporated a growth in the horse population that started with the onset of domestication as in the studies of Wutke et al. (, ), but effective population sizes (Ne’s) and generation time were considered noise variables because of limited information in the data for estimating them. We tested alternative demographic histories with no population growth and an ample range of Ne values. Therefore, the estimated parameters were the selection coefficients, the age of the Y-HT-1 haplotype, and, alternatively, its initial haplotype frequency (see below). The simulations incorporated two alternative models addressing the initial state of Y-HT-1: (i) It was forced to be present since the beginning of the simulated time frame at a frequency that was a parameter or (ii) it was allowed to occur, emulating its appearance by mutation, at some point between the beginning of the simulated time frame and its first detection (the age became the parameter here). Both models were compared by means of Bayes factors. The MCMC was run under a Metropolis-Hasting algorithm incorporating a Gibbs sampler. The priors and details of the S-MCMC procedure followed those by Wutke et al. (, ). We ran four chains, 500,000 steps each, with a burn-in of 50%. We also replicated the estimation of selection coefficients by means of approximate Bayesian computation () and the path augmentation method (). The statistical significance of the selection coefficients was assessed by obtaining the proportion of positive or negative values in the Bayesian posterior sample, for example, >95% of coefficients being positive could be interpreted as “significant” in a one-tailed test, and the same for negative values. In addition, we obtained Bayes factors for comparing a model with no selection (changes provided only by gene drift) and the model with selection upon Y-HT-1. Because selective forces might have changed over time (), we performed a model choice analysis for selecting the number of time bins (with different selection coefficient each) that best fitted our data in the S-MCMC analysis. We tested from one to six time bins by means of both the AIC and the BIC. Time bins were defined considering main historical ages and sample sizes. In addition, considering that other haplotypes could be under selection, we replicated the inference by S-MCMC in a multiallelic mode, in which Wright-Fisher simulations were performed by multinomial, rather than binomial, samplings and several haplotypes were allowed to show selection. We assessed the number of haplotypes with signs of selection by means of AIC and BIC. Finally, we reported the temporal trajectories of haplotype frequencies of the Bayesian posterior sample to assess the changes in frequencies over time. The simulations were programmed in the language Fortran by means of the platform Microsoft Visual Studio and Intel Parallel Studio, using the code reported by Wutke et al. (, ), which is available at www.dropbox.com/sh/uns1gsprj7xbbzo/AAAEr4hfya16P5nejpXVVnUfa?dl=0.