The origins and spread of domestic horses from the Western Eurasian steppes

Abstract

Domestication of horses fundamentally transformed long-range mobility and warfare¹. However, modern domesticated breeds do not descend from the earliest domestic horse lineage associated with archaeological evidence of bridling, milking and corralling^2-4 at Botai, Central Asia around 3500 BC³. Other longstanding candidate regions for horse domestication, such as Iberia⁵ and Anatolia⁶, have also recently been challenged. Thus, the genetic, geographic and temporal origins of modern domestic horses have remained unknown. Here we pinpoint the Western Eurasian steppes, especially the lower Volga-Don region, as the homeland of modern domestic horses. Furthermore, we map the population changes accompanying domestication from 273 ancient horse genomes. This reveals that modern domestic horses ultimately replaced almost all other local populations as they expanded rapidly across Eurasia from about 2000 BC, synchronously with equestrian material culture, including Sintashta spoke-wheeled chariots. We find that equestrianism involved strong selection for critical locomotor and behavioural adaptations at the GSDMC and ZFPM1 genes. Our results reject the commonly held association⁷ between horseback riding and the massive expansion of Yamnaya steppe pastoralists into Europe around 3000 BC^8,9 driving the spread of Indo-European languages¹⁰. This contrasts with the scenario in Asia where Indo-Iranian languages, chariots and horses spread together, following the early second millennium BC Sintashta culture^11,12.

Fig. 1. Ancient horse remains and their genomic affinities.

a, Temporal and geographic sampling. The red star indicates the location of the two TURG horses (late Yamnaya context) showing genetic continuity with DOM2. The dashed line indicates the inferred homeland of DOM2 horses in the lower Volga-Don region. Colours refer to regions and/or time periods delineating genetically close horses. The radius of each cylinder is proportional to the number of samples analysed (for <10 specimens; radius constant above this), and the height refers to the time range covered. b, Neighbour-joining phylogenomic tree (100 bootstrap pseudo-replicates). Samples are coloured according to a and the main phylogenetic clusters are numbered from 1 to 4. c, Fold difference between neighbour-joining-based and raw pairwise genetic distances. d, Pairwise distance matrix of Struct-f4 genetic affinities between samples. Increasing genetic affinities are indicated by a yellow-to-red gradient. e, Struct-f4 ancestry component profiles. f, Ancestry profiles of selected key horse groups and samples. PRZE, Przewalski; UP-SFR, Upper Palaeolithic Southern France.

Fig. 2. Horse geographic and genetic affinities.

a–c, EEMS-predicted migration barriers and average ancestry components found in each archaeological site from before 3000 bc (a), during the third millennium bc (b) and after around 2000 bc (c). The size of the pie charts is proportional to the number of samples analysed in a given location (<10, constant above). Pie chart colours refer to K = 6 ancestry components, averaged per location. Regions inferred as geographic barriers are shown in shades of brown, and regions affected by migrations are shown in shades of blue. The base map was obtained from rworldmap.

Fig. 3. Population genetic affinities, evolutionary history and geographic origins.

a, Multi-dimensional scaling plot of f₄-based genetic affinities. The age of the samples is indicated along the vertical axis. CA, Central Asia. b, Horse evolutionary history inferred by OrientAGraph with three migration edges and nine lineages representing key genomic ancestries (coloured as in Fig 1a). The model explains 99.99% of the total variance. The triangular pairwise matrix provides model residuals. The external branch leading to donkey was set to zero to improve visualization. c, LOCATOR predictions of the geographic region where the ancestors of DOM2, tarpan and modern Przewalski’s horses lived. The tarpan and modern Przewalski’s horses do not descend from the same ancestral population as modern domestic horses. The map was drawn using the maps R package.

Extended Data Fig. 1. Proportion of missing derived mutations at sites representing nucleotide transversions.

Proportions are provided relative to the genome of a modern Icelandic (P5782) horse (Spearman correlation coefficient between total transversion errors and time, R=−0.77 p-value =0).

Extended Data Fig. 2. Struct-f4 validation.

a, Simulated demographic model. A single migration pulse is assumed to have occurred 150 generations ago from population E into B. The magnitude of the migration represents 5% to 25% of the effective size of population B. The model was also simulated in the absence of migration (i.e. m=0%). Five individuals are simulated per population considered, except for the outgroup where only one individual was considered. b, Correlation of the expected levels of gene-flow with the predicted E-ancestry component in individuals i belonging to population B, as well as with the average Z-scores of the f₄(A, B_i; E, Outgroup) configurations, which reflects the stochasticity resulting from the simulations, prior to any inference. Each point represents a simulated individual. Colors indicate the 10 independent simulation replicates carried out. c, Predicted ancestry profiles in the absence (m=0%) and with gene flow (m=25% and K=7, as per the number of internal nodes immediately ancestral to the 10 extant populations).

Extended Data Fig. 3. Mobility and demographic shifts.

a–c, Correlation between observed pairwise genetic distances between demes as inferred by EEMS and Haversine geographic distances prior to ~3,000 BCE (a), during the third millennium BCE (b) and after ~2,000 BCE (c). d, Isolation-by-distance patterns through time inferred from autosomal (red) and X-chromosomal (blue) variation. e–f, Bayesian Skyline plots reconstructed from mtDNA (e) and Y-chromosomal variation (f). The third millennium BCE is highlighted in blue. The red line indicates the median of the 95% confidence range, shown in grey.

Extended Data Fig. 4. Individual ancestry profiles.

a, NJ-tree shown in Fig 1b with sample labels as defined in Supplementary Table 1. b, Struct-f4 individual ancestry profiles. c, Model likelihood. A total of K=4 to K=9 ancestral populations are assumed. LnL = natural log-likelihood.

Extended Data Fig. 5. OrientAGraph population histories and genetic distances to the domestic donkey.

a–e, OrientAGraph models and residuals assuming M=0 to M=5 migration edges and considering nine lineages representing key genomic ancestries (colored as in Fig 1a). M=3 is shown in Fig 3b. f, Pairwise genetic distances between a given horse and the domestic donkey plotted as a function of the age of the horse specimen considered.

Extended Data Fig. 6. Inter-regional trade and chariot networks, marked by horse cheek pieces, connecting Bronze Age steppe societies, mineral rich Caucasian societies and the Old Assyrian trade network during the period 1,950-1,750 BCE.

Documented Near Eastern trade routes are marked with stippled lines (after, supplemented with data from^, and Pavel F. Kuznetsov).

Extended Data Fig. 7. DOM2 selection signatures.

a, Manhattan plot of F_ST-differentiation index between DOM2 and non-DOM2 horses along the 31 EquCab3 autosomes. F_ST outliers are highlighted using an empirical P-value threshold of 10⁻⁵ (red dashed line). The two outlier regions on chromosomes 3 and 9 are highlighted within red frames. b, F_ST-differentiation index and genomic tracks around the ZFPM1 gene. Depth represents the accumulated number of reads per position within DOM2 (blue) and non-DOM2 (magenta) genomes. c, Same as Panel b at GSDMC.

Extended Data Fig. 8. Normalized read coverage supporting the presence of causative alleles for coat coloration variation.

Each column represents a particular genome position where genetic polymorphisms associated or causative for coat coloration patterns have been described. The exact EquCab3 genome coordinates are indicated in the locus label. Specimens (rows) are ordered according to their phylogenetic relationships, as shown in Fig 1b. The color gradient is proportional to the fraction of reads carrying the causative variant. Loci that are not covered following trimming and rescaling of individual BAM sequence alignment files are indicated with a white cross.