We gathered horse remains encompassing all suspected domestication centres, including Iberia, Anatolia and the steppes of Western Eurasia and Central Asia (Fig 1a). The sampling targeted previously under-represented time periods, with 201 radiocarbon dates spanning 44426 to 202bc, and five beyond 50250 to 47950bc (Supplementary Table1).

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 toa 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.

Full size image

The DNA quality enabled shotgun sequencing of 264 ancient genomes at 0.10× to 25.76× average coverage (239 genomes above 1× coverage), including 16 genomes for which further sequencing added to previously reported data. Enzymatic¹³ and computational removal of post mortem DNA damage produced high-quality data with derived mutations decreasing with sample age, as expected if mutations accumulate through time (Extended Data Fig.1). We added ten published modern genomes, and nine ancient genomes characterized with consistent technology or covering relevant time periods and locations, to obtain the most extensive high-quality genome time series for horses.

Neighbour-joining phylogenomic inference revealed four geographically defined monophyletic groups (Fig 1b). These closely mirrored clusters identified using an extension of the Struct-f4 method⁵ (Fig 1d–f, Extended Data Fig.2,Supplementary Methods), except for the Neolithic Anatolia group (NEO-ANA), where the tree-to-data goodness of fit suggested phylogenetic misplacement (Fig 1c,Supplementary Methods).

The most basal cluster includedEquus lenensis (ELEN), a lineage identified in northeastern Siberia from the Late Pleistocene to the late fourth millenniumbc^5,14,15. A second group covered Europe, including Late Pleistocene Romania, Belgium, France and Britain, and the region from Spain to Scandinavia and Hungary, Czechia and Poland during the sixth-to-third millenniumbc. The third cluster comprised the earliest known domestic horses from Botai and Przewalski’s horses, as previously reported³, and extended to the Altai and Southern Urals during the fifth-to-third millenniumbc. Finally, modern domestic horses clustered within a group that became geographically widespread and prominent following about 2200bc and during the second millenniumbc (DOM2). This cluster appears genetically close to horses that lived in the Western Eurasia steppes (WE) but not further west than the Romanian lower Danube, south of the Carpathians, before and during the third millenniumbc. Significant correlation between genetic and geographic distances, and inference of limited long-distance connectivity with estimated effective migration surface¹⁶ (EEMS), confirmed the strong geographic differentiation of horse populations before about 3000bc (Fig 2a, Extended Data Fig.3a).

a–c, EEMS-predicted migration barriers¹⁶ and average ancestry components found in each archaeological site from before 3000bc (a), during the third millenniumbc (b) and after around 2000bc (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 toK = 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⁴⁶.

Full size image

Horse ancestry profiles in Neolithic Anatolia and Eneolithic Central Asia, including at Botai, maximized a genetic component (coloured green in Fig.1e, f) that was also substantial in Central and Eastern Europe during the Late Pleistocene (RONPC06_Rom_m34801) and the fourth or third millenniumbc (Figs.1e,3a, Extended Data Fig.4). It was, however, absent or moderately present in the Romanian lower Danube (ENEO-ROM), the Dnieper steppes (Ukr11_Ukr_m4185) and the western lower Volga-Don (C-PONT) populations during the sixth to third millenniabc. This indicates possible expansions of Anatolian horses into both Central and Eastern Europe and Central Asia regions, but not into the Western Eurasia steppes. The absence of typical NEO-ANA ancestry rules out expansion from Anatolia into Central Asia across the Caucasus mountains but supports connectivity south of the Caspian Sea prior to about 3500bc.

a, Multi-dimensional scaling plot off₄-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⁴⁷.

Full size image

The C-PONT group not only possessed moderate NEO-ANA ancestry, but also was the first region where the typical DOM2 ancestry component (coloured orange in Fig.1e, f) became dominant during the sixth millenniumbc. Multi-dimensional scaling further identified three horses from the western lower Volga-Don region as genetically closest to DOM2, associated with Steppe Maykop (Aygurskii), Yamnaya (Repin) and Poltavka (Sosnovka) contexts, dated to about 3500 to 2600bc (Figs.2a, b,3a). Additionally, genetic continuity with DOM2 was rejected for all horses predating about 2200bc, especially those from the NEO-ANA group (Supplementary Table2), except for two late Yamnaya specimens from approximately 2900 to 2600bc (Turganik (TURG)), located further east than the western lower Volga-Don region (Figs.2a, b,3a). These may therefore have provided some of the direct ancestors of DOM2 horses.

Modelling of the DOM2 population with qpADM¹⁷, rotating¹⁸ all combinations of 2, 3 or 4 population donors, eliminated the possibility of a contribution from the NEO-ANA population, but indicated possible formation within the WE population, including a genetic contribution of approximately 95% from C-PONT and TURG horses (Supplementary Table3). This was consistent with OrientAGraph¹⁹ modelling from nine lineages representing key ancestry combinations, which confirmed the absence of NEO-ANA genetic ancestry in DOM2 and confirmed DOM2 as a sister population to the C-PONT horses (Fig.3b).

Identifying discrete populations and modelling admixture as single unidirectional pulses, however, was highly challenging given the extent of spatial genetic connectivity. Indeed, the typical DOM2 ancestry component was maximized in the C-PONT group, but declined sharply eastwards (TURG and Central Asia) in the third millenniumbc as the proportion of NEO-ANA ancestry increased (Fig.2a). This suggests a cline of genetic connectivity east of the Western Eurasia steppes and Central Asia, ruling out DOM2 ancestors further east than the western lower Volga-Don and Turganik. A similar genetic cline characterized the region located west of C-PONT, where the typical DOM2 ancestry component declined steadily in the Dnieper steppes, Poland, Turkish Thrace and Hungary in the fifth to third millenniabc. This eliminates the possibility of DOM2 ancestors further west than C-PONT and the Dnieper steppes. Furthermore, patterns of spatial autocorrelations in the genetic data²⁰ indicated Western Eurasia steppes as the most likely geographic location of DOM2 ancestors (Fig.3c). Combined, our results demonstrate that DOM2 ancestors lived in the Western Eurasia steppes, especially the lower Volga-Don, but not in Anatolia, during the late fourth and early third millenniabc.

Analyses of ancient human genomes have revealed a massive expansion from the Western Eurasia steppes into Central and Eastern Europe during the third millenniumbc, associated with the Yamnaya culture^8,9,11,12,21. This expansion contributed at least two thirds of steppe-related ancestry to populations of the Corded Ware complex (CWC) around 2900 to 2300bc⁸. The role of horses in this expansion remained unclear, as oxen could have pulled Yamnaya heavy, solid-wheeled wagons^7,22. The genetic profile of horses from CWC contexts, however, almost completely lacked the ancestry maximized in DOM2 and Yamnaya horses (TURG and Repin) (Figs.1e, f,2a, b) and showed no direct connection with the WE group, including both C-PONT and TURG, in OrientAGraph modelling (Fig.3b, Extended Data Fig.5).

The typical DOM2 ancestry was also limited in pre-CWC horses from Denmark, Poland and Czechia, associated with the Funnel Beaker and early Pitted Ware cultures (FB/PWC, FB/POL and ENEO-CZE, respectively). DOM2 ancestry reached a maximum 12.5% in one Hungarian horse dated to the mid-third millenniumbc and associated with the Somogyvár-Vinkovci Culture (CAR05_Hun_m2458). qpAdm¹⁷ modelling indicated that its DOM2 ancestry was acquired following gene flow from southern Thrace (Kan22_Tur_m2386), but not from the Dnieper steppes (Ukr11_Ukr_m4185) (Supplementary Table3). Combined with the lack of increased horse dispersal during the early third millenniumbc (Fig.2b, Extended Data Fig.3b), these results suggest that DOM2 horses did not accompany the steppe pastoralist expansion north of the Carpathians.

By around 2200–2000bc, the typical DOM2 ancestry profile appeared outside the Western Eurasia steppes in Bohemia (Holubice), the lower Danube (Gordinesti II) and central Anatolia (Acemhöyük), spreading across Eurasia shortly afterwards, eventually replacing all pre-existing lineages (Fig 2c, Extended Data Fig.3c). Eurasia became characterized by high genetic connectivity, supporting massive horse dispersal by the late third millennium and early second millenniumbc. This process involved stallions and mares, indicated by autosomal and X-chromosomal variation (Extended Data Fig.3d), and was sustained by explosive demographics apparent in both mitochondrial and Y-chromosomal variation (Extended Data Fig.3e, f). Altogether, our genomic data uncover a high turnover of the horse population in which past breeders produced large stocks of DOM2 horses to supply increasing demands for horse-based mobility from around 2200bc.

Of note, the DOM2 genetic profile was ubiquitous among horses buried in Sintashta kurgans together with the earliest spoke-wheeled chariots around 2000–1800bc^7,9,23,24 (Extended Data Fig.6). A typical DOM2 profile was also found in Central Anatolia (AC9016_Tur_m1900), concurrent with two-wheeled vehicle iconography from about 1900bc^25,26. However, the rise of such profiles in Holubice, Gordinesti II and Acemhöyük before the earliest evidence for chariots supports horseback riding fuelling the initial dispersal of DOM2 horses outside their core region, in line with Mesopotamian iconography during the late third and early second millenniabc²⁷. Therefore, a combination of chariots and equestrianism is likely to have spread the DOM2 diaspora in a range of social contexts from urban states to dispersed decentralized societies²⁸.

Human-induced DOM2 dispersal conceivably involved selection of phenotypic characteristics linked to horseback riding and chariotry. We therefore screened our data for genetic variants that are over-represented in DOM2 horses from the late third millenniumbc (Extended Data Fig.7). The first outstanding locus peaked immediately upstream of theGSDMC gene, where sequence coverage dropped at two L1 transposable elements in all lineages except DOM2. The presence of additional exons in other mammals suggests that independent L1 insertions remodelled the DOM2 gene structure. In humans,GSDMC is a strong marker for chronic back pain²⁹ and lumbar spinal stenosis, a syndrome causing vertebral disk hardening and painful walking³⁰.

The second most differentiated locus extended over approximately 16 Mb on chromosome 3, with theZFPM1 gene being closest to the selection peak.ZFPM1 is essential for the development of dorsal raphe serotonergic neurons involved in mood regulation³¹ and aggressive behaviour³².ZFPM1 inactivation in mice causes anxiety disorders and contextual fear memory³¹. Combined, early selection atGSDMC andZFPM1 suggests shifting use toward horses that were more docile, more resilient to stress and involved in new locomotor exercise, including endurance running, weight bearing and/or warfare.

Our analyses elucidate the geographic, temporal and biological origins of DOM2 horses. This study features a diverse ancient horse genome dataset, revealing the presence of deep mitochondrial and/or Y-chromosomal haplotypes in non-DOM2 horses (Supplementary Fig 1). This suggests that yet-unsampled divergent populations contributed to forming several lineages excluding DOM2. This was especially true in the Iberian group (IBE), where the expected genetic distance to the donkey was reduced (Extended Data Fig.5f), but also in NEO-ANA according to OrientAGraph modelling (Fig 3b). Disentangling exact divergence and ancestry contributions of such unsampled lineages is difficult with the currently available data. It can, however, be stressed that Iberia and Anatolia represent two well-known refugia³³, where populations could have survived and mixed during Ice Ages.

Finally, our analyses have solved the mysterious origins of the tarpan horse, which became extinct in the early 20th century. The tarpan horse came about following admixture between horses native to Europe (modelled as having 28.8–34.2% and 32.2–33.2% CWC ancestry in OrientAGraph¹⁹ and qpAdm¹⁷, respectively) and horses closely related to DOM2. This is consistent with LOCATOR²⁰ predicting ancestors in western Ukraine (Fig 3c) and refutes previous hypotheses depicting tarpans as the wild ancestor or a feral version of DOM2, or a hybrid with Przewalski’s horses³⁴.

This work resolves longstanding debates about the origins and spread of domestic horses. Whereas horses living in the Western Eurasia steppes in the late fourth and early third millenniabc were the ancestors of DOM2 horses, there is no evidence that they facilitated the expansion of the human genetic steppe ancestry into Europe^8,9 as previously hypothesized⁷. Instead of horse-mounted warfare, declining populations during the European late Neolithic³⁵ may thus have opened up an opportunity for a westward expansion of steppe pastoralists. Yamnaya horses at Repin and Turganik carried more DOM2 genetic affinity than presumably wild horses from hunter-gatherer sites of the sixth millenniumbc (NEO-NCAS, from approximately 5500–5200bc), which may suggest early horse management and herding practices. Regardless, Yamnaya pastoralism did not spread horses far outside their native range, similar to the Botai horse domestication, which remained a localized practice within a sedentary settlement system^2,36. The globalization stage started later, when DOM2 horses dispersed outside their core region, first reaching Anatolia, the lower Danube, Bohemia and Central Asia by approximately 2200 to 2000bc, then Western Europe and Mongolia soon afterwards, ultimately replacing all local populations by around 1500 to 1000bc. This process first involved horseback riding, as spoke-wheeled chariots represent later technological innovations, emerging around 2000 to 1800bc in the Trans-Ural Sintashta culture⁷. The weaponry, warriors and fortified settlements associated with this culture may have arisen in response to increased aridity and competition for critical grazing lands, intensifying territoriality and hierarchy³⁷. This may have provided the basis for the conquests over the subsequent centuries that resulted in an almost complete human and horse genetic turnover in Central Asian steppes^11,21. The expansion to the Carpathian basin³⁸, and possibly Anatolia and the Levant, involved a different scenario in which specialized horse trainers and chariot builders spread with the horse trade and riding. In both cases, horses with reduced back pathologies and enhanced docility would have facilitated Bronze Age elite long-distance trade demands and become a highly valued commodity and status symbol, resulting in rapid diaspora. We, however, acknowledge substantial spatiotemporal variability and evidential bias towards elite activities, so we do not discount additional, harder to evidence, factors in equine dispersal.

Our results also have important implications for mechanisms underpinning two major language dispersals. The expansion of the Indo-European language family from the Western Eurasia steppes has traditionally been associated with mounted pastoralism, with the CWC serving as a major stepping stone in Europe^39,40,41. However, while there is overwhelming lexical evidence for horse domestication, horse-drawn chariots and derived mythologies in the Indo-Iranian branch of the Indo-European family, the linguistic indications of horse-keeping practices at the deeper Proto-Indo-European level are in fact ambiguous⁴² (Supplementary Discussion) . The limited presence of horses in CWC assemblages⁴³ and the local genetic makeup of CWC specimens reject scenarios in which horses were the primary driving force behind the initial spread of Indo-European languages in Europe⁴⁴. By contrast, DOM2 dispersal in Asia during the early-to-mid second millenniumbc was concurrent with the spread of chariotry and Indo-Iranian languages, whose earliest speakers are linked to populations that directly preceded the Sintashta culture^11,12,45. We thus conclude that the new package of chariotry and improved breed of horses, including chestnut coat colouration documented both linguistically (Supplementary Discussion) and genetically (Extended Data Fig.8), transformed Eurasian Bronze Age societies globally within a few centuries after about 2000bc. The adoption of this new institution, whether for warfare, prestige or both, probably varied between decentralized chiefdoms in Europe and urbanized states in Western Asia. The results thus open up new research avenues into the historical developments of these different societal trajectories.

Radiocarbon dating

A total of 170 new radiocarbon dates were obtained in this study. Dating was carried out at the Keck Carbon Cycle AMS Laboratory, UC Irvine following collagen extraction and ultra-filtration from approximately 1 g of osseous material. IntCal20 calibration⁴⁸ was performed using OxCalOnline⁴⁹.

Genome sequencing

All samples were collected with permission from the organizations holding the collections and documented through official authorization letters for partially destructive sampling from local authorities. Samples were processed for DNA extraction, library construction and shallow sequencing in the ancient DNA facilities of the Centre for Anthropobiology and Genomics of Toulouse (CAGT), France. The overall methodology followed the work from Seguin–Orlando and colleagues⁵⁰. It involved: (1) powdering a total of 100–590 mg of osseous material using the Mixel Mill MM200 (Retsch) Micro-dismembrator; (2) extracting DNA following the procedure Y2 from Gamba and colleagues⁵¹, tailored to facilitate the recovery of even the shortest DNA fragments; (3) treating DNA extracts with the USER (NEB) enzymatic cocktail to eliminate a fraction of post mortem DNA damage¹³; (4) constructing from double-stranded DNA templates DNA libraries in which two internal indexes are added during adapter ligation and one external index is added during PCR amplification; and (5) amplification, purification and quantification of DNA libraries before pooling 20–50 DNA libraries for low-depth sequencing on the Illumina MiniSeq instrument (paired-end mode, 2 × 80). All three indexes of each library were unique in a given sequencing pool.

Raw fastQ files were demultiplexed, trimmed and collapsed when individual read pairs showed significant overlap using AdapterRemoval2⁵² (version 2.3.0), disregarding reads shorter than 25 bp. Processed reads were then aligned against the nuclear and mitochondrial horse reference genomes^53,54, and appended with the Y-chromosome contigs from⁵⁵ using the Paleomix bam_pipeline (version 1.2.13.2) with the mapping parameters recommended by Poullet and Orlando⁵⁶. Sequencing reads representing PCR duplicates or showing a mapping quality below 25 were disregarded. DNA fragmentation and nucleotide misincorporation patterns were assessed on the basis of 100,000 random mapped reads using mapDamage2⁵⁷ (version 2.0.8). Paleomix returned provisional estimates of endogenous DNA content and clonality, as defined by the fraction of retained reads mapping uniquely against the horse reference genomes and those mapping at the same genomic coordinates, respectively. These numbers guided further experimental decisions, including (1) the sequencing effort to be performed per individual library; (2) the preparation of additional libraries from left-over aliquots of USER-treated DNA extracts, or following treatment of DNA extract aliquots with the USER enzymatic cocktail; and (3) the preparation of additional DNA extracts. After initial screening for library content, sequencing was carried out on the Illumina HiSeq4000 instruments from Genoscope (paired-end mode, 2 × 76; France Génomique), except for four samples (BPTDG1_Fra_m11800, Closeau3_Fra_m10400, Novoil1_Kaz_m1832 and Novoil2_Kaz_m1832), for which sequencing was done at Novogene Europe on an Illumina NovaSeq 6000 instrument (S4 lanes, paired-end mode, 2 × 150). Overall, we obtained sequence data for a total of 264 novel ancient horse specimens and 1,029 DNA libraries (980 new), summing up to 31.86 billion sequencing read pairs and 100.82 billion collapsed read pairs, which was sufficient to characterize 226 novel ancient genomes showing a genomic depth-of-coverage of at least 1× (median 2.80-fold, maximum 25.76-fold) (Supplementary Table1).

Allele sampling, sequencing error rates, genome rescaling and trimming

Following previous work^5,58, error rates are defined as the excess of mutations that are private to the ancient genome, relative to a modern genome considered as error-free. Mutations were polarized using an outgroup genome representing a consensus built from seven male specimens of diverse equine species (Equus africanus somaliensis,Equus asinus,Equus burchelli,Equus grevyi,Equus hartmannae,Equus hemionus onager andEquus kiang⁵⁹), according to a majority rule in which at maximum 2 of the 7 individuals showed an alternative allele. Minor and major alleles were identified using ANGSD⁶⁰ (version 0.933-86-g3fefdc4, htslib: 1.10.2-106-g9c35744) and the following parameters: -baq 0 -doMajorMinor 2 -uniqueOnly 1 -minMapQ 25 -minQ 30 -minind 7 -doCounts 1 -doMaf 1.

Error rate estimates ranged between 0.000337 and 0.003966 errors per site and revealed that nucleotide C→T and G→A misincorporation rates were still inflated relative to their reciprocal substitution types (T→C and A→G), despite USER treatment. Therefore, individual BAM alignment files were processed to further reduce nucleotide misincorporation rates. To achieve this, we used PMDtools⁶¹ (version 0.60) to bin apart reads likely containing post mortem DNA damage (--threshold 1; DAM) from those that did not (--upperthreshold 1; NODAM). NODAM-aligned reads were then directly trimmed by 5 bp at their ends, where individual base qualities generally drop. The base quality of aligned DAM reads was first rescaled using mapDamage2⁵⁷ (version 2.0.8), penalizing all instances of potential derivatives of post mortem cytosine deamination, then further trimmed by 10 bp at both ends. The resulting NODAM and DAM aligned reads were merged again to obtain final BAM sequence alignments. Final error rate estimates ranged between 0.000080 and 0.000933 errors per site (Supplementary Table1).

Uniparentally inherited markers and coat colouration

Mitochondrial genomes for the 264 newly sequenced samples were characterized from quality-filtered BAM alignment files (minMapQ=25, minQ=30), using a majority rule requiring at least five individual reads per position. Their resulting complete mitochondrial genome sequences were aligned together with a total of 193 sequences previously characterized^{3,5,14,15,58,62,63} using mafft⁶⁴ (version 7.407). Sequence alignments were split into six partitions, following previous work⁵, including the control region, all tRNAs, both rRNAs and each codon position considered separately. Maximum-likelihood phylogenetic reconstruction was performed using RAxML⁶⁵ (version 8.2.11) with default parameters, and assessing node support from a total of 100 bootstrap pseudo-replicates. The same partitions were provided as input for BEAST⁶⁶ (version 2.5.1), together with calibrated radiocarbon years (Supplementary Table1). Specimens lacking direct radiocarbon dates or identified as not belonging to the DOM2 cluster were disregarded (Supplementary Table1). While the former ensured precise tip-calibration for molecular clock estimation (assuming uncorrelated log-normal relaxed model), the latter prevented misinterpreting spatial variation in the population structure as changes in the effective population size⁶⁷. The best substitution model was selected from ModelGenerator⁶⁸ (version 0.85) and Bayesian Skyline plots⁶⁹ were retrieved following 1,000,000,000 generations, sampling 1 every 1,000 and disregarding the first 30% as burn-in. Convergence was visually checked in Tracer⁷⁰ (version 1.7.2).

The Y-chromosome maximum-likelihood tree was constructed calling individual haplotypes from trimmed and rescaled BAM sequence alignments against the contigs described by Felkel and colleagues⁵⁵, filtered for single copy MSY regions. The final multifasta sequence alignment included sites covered in at least 20% of the specimens, pseudo-haploidizing each position and filtering out transitions, as done with autosomal data. It was further restricted to specimens showing at least 20% of the final set of positions covered. This represented a total of 3,195 nucleotide transversions for 142 specimens. The final tree was computed using IQtree (version 1.6.12), following AICc selection of the best substitution model and 1,000 ultrafast bootstrap approximation for assessing node support^71,72. The Y-chromosome Bayesian skyline plot was obtained following the same procedure as above. Maximum-likelihood trees and Bayesian skyline plots are shown in Supplementary Fig 1 and Extended Data Fig.3e, f, respectively.

The presence of alleles associated with or causative for a diversity of coat colouration changes was investigated using individual BAM read alignments. For a total of 43 genomic locations representing biallelic SNPs, we simply counted the proportion of reads supporting the associated or causative allele. Results were summarized in the heat map shown in Extended Data Fig.8, with respect to the sample ordering displayed in the neighbour-joining phylogenetic reconstruction, and limited to those 13 loci that were polymorphic in our horse panel for clarity.

Neighbour-joining phylogeny, genetic continuity and population modelling

Phylogenetic affinities were first estimated by performing a BioNJ tree reconstruction with FastME⁷³ (version 2.1.4), based on the pairwise matrix of genetic distances inferred from the bed2diffs_v1 program¹⁶. Node supports were assessed using a total of 100 bootstrap pseudo-replicates. The ‘goodness-of-fit’ of the neighbour-joining tree to the data was evaluated by comparing the patristic distances and raw pairwise distances. Patristic distances were obtained from the ape⁷⁴ R package (version 5.5) and their ratios to raw pairwise distances were averaged for each given individual (Fig 1c). Averaged ratios equal to one support perfect phylogenetic placement for the specimen considered.

Genetic continuity between each individual specimen predating about 2200bc and DOM2 horses was tested following the methodology from Schraiber⁷⁵, which implements a likelihood-ratio test to compare the statistical support for placing DOM2 and the ancient specimen in a direct line of ancestry or as two sister groups. This methodology relies on exact allele frequency estimates within DOM2 and read counts for putatively ancestral ancient samples. To exclude residual sequencing errors within DOM2 horses, we, thus, conditioned these analyses on variants segregating at least as doubletons in positions covered in at least 75% of the DOM2 samples. Linked variation was pruned using Plink⁷⁶ (version v.1.9), with the following parameters, --indep-pairwise 50 10 0.2, which provided a panel of about 1.4 million transversions. Allele frequencies were polarized considering the outgroup genome used for measuring error rates. Results from direct ancestry tests are summarized in Supplementary Table2.

The complex genetic makeup of some individuals (CAR05_Hun_m2458 and Kan22_Tur_m2386) and/or group of individuals (DOM2) was investigated using thef₄-statistics-based ancestry decomposition approach implemented in qpAdm¹⁷ (version 7.0), in which one particular (group of) individual(s) is modelled as a linear, additive combination of candidate population sources (‘left’ populations). We followed the rotating strategy recommended by Harney and colleagues¹⁸ to assess all possible combinations of two, three and four donors (‘left’) selected from a total of 18 populations. The remaining 14, 15 and 16 populations were used as reference (‘right’) populations (Supplementary Table3).

We selected a total of nine horse lineages representing the main phylogenetic clusters, and carrying genetic ancestry profiles representative of the complete dataset, to model the population evolutionary history using OrientAGraph¹⁹ (version 1.0). By implementing a network orientation subroutine that enables throughout exploration of the graph space, OrientAGraph constitutes a marked advancement in the automated inference of admixture graphs. We considered scenarios from zero to five migration pulses (M = 0 to 5; Extended Data Fig.5a–e), and the population model assumingM = 3 is represented in Fig 3b. This analysis was conditioned on sites covered at least in one specimen of each population group. This filter yielded a set of 7,936,493 fully orthologous nucleotide transversions.

Struct-f4, ancestry components and multi-dimensional scaling

We extended the Struct-f4 package so as to assess individual genetic affinities within a panel of genomes, and to decompose them intoK genetic ancestries. Struct-f4, thus, achieves similar objectives to other clustering methods, such as ADMIXTURE⁷⁷ and Ohana⁷⁸, but does not assume Hardy–Weinberg equilibrium. The latter assumption is known to cause misinterpretation of highly drifted samples as ancestral homogeneous groups instead of highly derived mixtures from multiple populations, as thoroughly described elsewhere⁷⁹. To circumvent this, Struct-f4 relies on the calculation of the widely usedf₄ statistics, which were originally devised not only to test for admixture, but also to quantify the drift between the internal nodes of a population tree. The latter provides a direct representation of the true ancestral populations. Overall, Struct-f4 thus implements a more natural and robust (model-free) approach than other clustering alternatives.

Struct-f4 is based on a mixture model that parametrizes the drift that occurred between a given number ofK pre-defined ancestral populations, and the mixing coefficient of each individual. Model parameters are estimated using an adaptive Metropolis–Hastings Markov chain Monte Carlo integration, identifying optimal numerical solutions for parameters by means of likelihood maximization. Struct-f4 was validated following extensive coalescent simulations with fastsimcoal2⁸⁰ (version 2.6.0.3). An example of such simulation designed to mimic the complex horse evolutionary history is provided in Extended Data Fig.2, based on mutation and recombination rates of 2.3 × 10⁻⁸ and 10⁻⁸ events per generation and bp, respectively. Struct-f4 is implemented in Rcpp and only takes the full set off₄-statistics as input to automatically return individual ancestry coefficients, without requiring pre-defined, ad-hoc sets of reference and test populations.

Multi-dimensional scaling was carried out based on the co-ancestry semi-matrix summarizing the drift measured between each pair of individuals, as returned by Struct-f4, removing the domestic donkey outgroup prior to using the cmdscale R function.

Isolation by distance and spatial connectivity

Spatial barriers to gene flow prior to about 3000bc, between about 3000 and 2000bc and following about 2000bc were run using EEMS¹⁶ (built with Eigen version 3.2.2 and Boost version 1.57, and using rEEMSplots version 0.0.0.9000) for 50 million iterations and considering a burn-in of 15 million iterations. Convergence was ensured from visual inspection of likelihood trajectories as well as by the strong correlation obtained between the observed and fitted genetic dissimilarities. Pie-charts depicting the ancestry proportions inferred by Struct-f4 were overlaid on the migration surfaces to facilitate tracking the geographic position of each excavation site, averaging ancestry proportions or using individual ancestry profiles if only one sample was characterized genetically at that location. Spatial pie-chart projection was carried out using the draw.pie R function from the mapplots package⁸¹ (version 1.5.1). The size of each individual pie-chart was commensurate with the number of samples excavated at a given geographic location, provided that the number of samples was lower than 10, while set to a constant maximum radius otherwise.

Partial Mantel tests measuring the correlation between geographic and genomic distances over time were carried out using the ncf R package⁸² (version 1.2.9). This test corrected for the time variation present within each window, similar to the approach described by Loog and colleagues⁸³. Haversine geographic distances between pairs of ancient samples were computed using the geosphere package (version 1.5.10) in R⁸⁴, from the corresponding longitude and latitude coordinates, while radiocarbon date ages were considered as point estimates (Supplementary Table1). The matrix of pairwise genetic distances was obtained from the bed2diffs_v1 program provided together with the EEMS software¹⁶. The analysis was carried out for autosomes and the X chromosome separately, so as to investigate possible sex-bias in horse dispersal. Confidence intervals were calculated by sampling with replacement individuals within each time window.

Sliding time windows (step size = 250 years) were broadened forward in time until including at least ten specimens covering two-thirds of the total geographic area sampled in this study. The area delimited by a set or subset of GPS coordinates was calculated using the GeoRange R package⁸⁵ (version 0.1.0) and the age of the window was set to the average age amongst the samples included. Additionally, pairwise distances involving samples located less than 500 km away and separated by less than 500 years were masked in the corresponding matrices to estimate the patterns of isolation by distance between demes, instead of within demes. This whole scheme was designed to prevent regional effects, caused by the over-representation of particular regions in specific time intervals.

The LOCATOR²⁰ program (version 1.2) was run using a geolocated reference panel consisting of all non-DOM2 horses (n = 136), except the tarpan and the four Przewalski’s horses present in our dataset, and considering nucleotide transversions covered at least in 75% of the samples, for a total of 3,194,008 SNPs. The geographic origin of each DOM2 horse was then estimated from the geographic structure defined by the populations present in the reference panel. Default parameters were used, except that the width of each neural layer was 512 (instead of 256). The best run was selected as the one showing the lowest validation error from a total of 50 independent runs. The analysis was repeated for the tarpan as well as the four Przewalski’s horses present in our dataset.

Selection scans

To pinpoint genetic changes potentially underlying biological adaptation within DOM2 horses, we contrasted the frequency of each nucleotide transversion in our dataset (n = 10,205,277) in DOM2 (n = 141) and non-DOM2 horses (n = 142). The extensive number of samples represented provided unprecedented resolution into patterns of allele frequency differentiation, and encompassed the largest diversity of non-DOM2 horses characterized to date. Weir and CockerhamF_ST index values between both groups were calculated using Plink⁷⁶ (version 1.9) and visualized using the GViz R package⁸⁶ (version 1.36.2), together with external genomic tracks provided by the gene models annotated for EquCab3 (Ensembl v0.102) and the interrupted repeats precomputed for the same assembly and stored in the UCSC browser.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

All collapsed and paired-end sequence data for samples sequenced in this study are available in compressed fastq format through the European Nucleotide Archive under accession numberPRJEB44430, together with rescaled and trimmed bam sequence alignments against both the nuclear and mitochondrial horse reference genomes. Previously published ancient data used in this study are available under accession numbersPRJEB7537,PRJEB10098,PRJEB10854,PRJEB22390 andPRJEB31613, and detailed in Supplementary Table1. The genomes of ten modern horses, publicly available, were also accessed as indicated in their corresponding original publications^{59,63,87,88,89}.

The Struct-f4 software is available without restriction on Bitbucket (https://bitbucket.org/plibradosanz/structf4/src/master/).

We thank all members of the AGES group at CAGT. We are grateful for the Museum of the Institute of Plant and Animal Ecology (UB RAS, Ekaterinburg) for providing specimens. The work by G. Boeskorov is done on state assignment of DPMGI SB RAS. This project was supported by the University Paul Sabatier IDEX Chaire d’Excellence (OURASI); Villum Funden miGENEPI research programme; the CNRS ‘Programme de Recherche Conjoint’ (PRC); the CNRS International Research Project (IRP AMADEUS); the France Génomique Appel à Grand Projet (ANR-10-INBS-09-08, BUCEPHALE project); IB10131 and IB18060, both funded by Junta de Extremadura (Spain) and European Regional Development Fund; Czech Academy of Sciences (RVO:67985912); the Zoological Institute ZIN RAS (АААА-А19-119032590102-7); and King Saud University Researchers Supporting Project (NSRSP–2020/2). The research was carried out with the financial support of the Russian Foundation for Basic Research (19-59-15001 and 20-04-00213), the Russian Science Foundation (16-18-10265, 20-78-10151, and 21-18-00457), the Government of the Russian Federation (FENU-2020-0021), the Estonian Research Council (PRG29), the Estonian Ministry of Education and Research (PRG1209), the Hungarian Scientific Research Fund (Project NF 104792), the Hungarian Academy of Sciences (Momentum Mobility Research Project of the Institute of Archaeology, Research Centre for the Humanities); and the Polish National Science Centre (2013/11/B/HS3/03822). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie (grant agreement 797449). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements 681605, 716732 and 834616).

1. Naveed Khan

   Present address: Lundbeck Foundation GeoGenetics Centre, Copenhagen, Denmark

Authors and Affiliations

1. Centre d’Anthropobiologie et de Génomique de Toulouse, Université Paul Sabatier, Toulouse, France

   Pablo Librado, Naveed Khan, Antoine Fages, Mariya A. Kusliy, Tomasz Suchan, Laure Tonasso-Calvière, Stéphanie Schiavinato, Duha Alioglu, Aurore Fromentier, Charleen Gaunitz, Lorelei Chauvey, Andaine Seguin-Orlando, Clio Der Sarkissian & Ludovic Orlando

2. Department of the Diversity and Evolution of Genomes, Institute of Molecular and Cellular Biology SB RAS, Novosibirsk, Russia

   Mariya A. Kusliy & Alexander S. Graphodatsky

3. W. Szafer Institute of Botany, Polish Academy of Sciences, Kraków, Poland

   Tomasz Suchan & Magdalena Moskal-del Hoyo

4. Genoscope, Institut de biologie François-Jacob, Commissariat à l’Energie Atomique (CEA), Université Paris-Saclay, Evry, France

   Aude Perdereau

5. Génomique Métabolique, Genoscope, Institut de biologie François Jacob, CEA, CNRS, Université d’Evry, Université Paris-Saclay, Evry, France

   Jean-Marc Aury & Patrick Wincker

6. Earth System Science Department, University of California, Irvine, Irvine, CA, USA

   John Southon

7. Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, Santa Cruz, CA, USA

   Beth Shapiro

8. Howard Hughes Medical Institute, University of California, Santa Cruz, Santa Cruz, CA, USA

   Beth Shapiro

9. Department of Archaeology, Ethnography and Museology, Altai State University, Barnaul, Russia

   Alexey A. Tishkin, Kirill Yu. Kiryushin & Nikolai N. Seregin

10. Department of Archaeological Heritage Preservation, Institute of Archaeology of the Russian Academy of Sciences, Moscow, Russia

    Alexey A. Kovalev

11. Zoology Department, College of Science, King Saud University, Riyadh, Saudi Arabia

    Saleh Alquraishi, Ahmed H. Alfarhan & Khaled A. S. Al-Rasheid

12. Institute for Archaeology, Heritage Conservation Studies and Art History, University of Bamberg, Bamberg, Germany

    Timo Seregély

13. Museum Østjylland, Randers, Denmark

    Lutz Klassen

14. Saxo Institute, section of Archaeology, University of Copenhagen, Copenhagen, Denmark

    Rune Iversen

15. ArScAn-UMR 7041, Equipe Ethnologie préhistorique, CNRS, MSH-Mondes, Nanterre Cedex, France

    Olivier Bignon-Lau, Pierre Bodu & Monique Olive

16. Muséum d’histoire naturelle, Secteur des Vertébrés, Geneva, Switzerland

    Jean-Christophe Castel

17. UMR 5199 De la Préhistoire à l’Actuel : Culture, Environnement et Anthropologie (PACEA), CNRS, Université de Bordeaux, Pessac Cedex, France

    Myriam Boudadi-Maligne & Mélanie Pruvost

18. Geneva Natural History Museum, Geneva, Switzerland

    Nadir Alvarez

19. Department of Genetics and Evolution, University of Geneva, Geneva, Switzerland

    Nadir Alvarez

20. OD Earth & History of Life, Royal Belgian Institute of Natural Sciences, Brussels, Belgium

    Mietje Germonpré

21. Institute of Systematics and Evolution of Animals, Polish Academy of Sciences, Kraków, Poland

    Jarosław Wilczyński & Sylwia Pospuła

22. Institute of Archaeology and Ethnology Polish Academy of Sciences, Kraków, Poland

    Anna Lasota-Kuś & Krzysztof Tunia

23. Institute of Archaeology, Jagiellonian University, Kraków, Poland

    Marek Nowak

24. Department of Archaeology, Institute of History and Archaeology, Tartu, Estonia

    Eve Rannamäe

25. Department of Zoology, Institute of Ecology and Earth Sciences, University of Tartu, Tartu, Estonia

    Urmas Saarma

26. Diamond and Precious Metals Geology Institute, SB RAS, Yakutsk, Russia

    Gennady Boeskorov

27. Archaeological Research Collection, Tallinn University, Tallinn, Estonia

    Lembi Lōugas

28. Department of Natural Sciences and Archaeometry, Institute of Archaeology of the Czech Academy of Sciences, Prague, Czechia

    René Kyselý

29. Prague, Czechia

    Lubomír Peške

30. Vasile Pârvan Institute of Archaeology, Department of Bioarchaeology, Romanian Academy, Bucharest, Romania

    Adrian Bălășescu, Valentin Dumitrașcu & Roxana Dobrescu

31. Institute of Archaeogenomics, Research Centre for the Humanities, Eötvös Loránd Research Network, Budapest, Hungary

    Daniel Gerber, Anna Szécsényi-Nagy & Balázs G. Mende

32. Department of Genetics, Eötvös Loránd University, Budapest, Hungary

    Daniel Gerber

33. Institute of Archaeology, Research Centre for the Humanities, Eötvös Loránd Research Network, Budapest, Hungary

    Viktória Kiss, Gabriella Kulcsár & Erika Gál

34. Ásatárs Ltd., Kecskemét, Hungary

    Zsolt Gallina

35. Rippl-Rónai Municipal Museum with Country Scope, Kaposvár, Hungary

    Krisztina Somogyi

36. School of History, Classics and Archaeology, University of Edinburgh, Old Medical School, Edinburgh, UK

    Robin Bendrey

37. Trace and Environmental DNA (TrEnD) Lab, School of Molecular and Life Sciences, Curtin University, Perth, Western Australia, Australia

    Morten E. Allentoft

38. Lundbeck Foundation GeoGenetics Centre, GLOBE Institute, University of Copenhagen, Copenhagen, Denmark

    Morten E. Allentoft

39. Department of Academic Management, Academy of Science of Moldova, Chișinău, Republic of Moldova

    Ghenadie Sirbu

40. Center of Archaeology, Institute of Cultural Heritage, Academy of Science of Moldova, Chișinău, Republic of Moldova

    Valentin Dergachev

41. Archaeological Institute of America, Boston, MA, USA

    Henry Shephard

42. Centre National de Recherche Scientifique, Muséum national d’Histoire naturelle, Archéozoologie, Archéobotanique (AASPE), CP 56, Paris, France

    Noémie Tomadini, Sandrine Grouard, Benoit Clavel, Sébastien Lepetz & Marjan Mashkour

43. Institute for the History of Material Culture, Russian Academy of Sciences (IHMC RAS), St Petersburg, Russia

    Aleksei Kasparov, Vladimir Pitulko, Alexander Bessudnov & Nikolay A. Bokovenko

44. Geological Institute, Russian Academy of Sciences, Moscow, Russia

    Alexander E. Basilyan & Pavel A. Nikolskiy

45. Arctic and Antarctic Research Institute, St Petersburg, Russia

    Mikhail A. Anisimov & Elena Y. Pavlova

46. Institute of Animal Breeding and Genetics, University of Veterinary Medicine Vienna, Vienna, Austria

    Gottfried Brem & Barbara Wallner

47. Department of Prehistory and Western Asian/Northeast African Archaeology, Austrian Archaeological Institute, Austrian Academy of Sciences, Vienna, Austria

    Christoph Schwall

48. Estonian Biocentre, Institute of Genomics, University of Tartu, Tartu, Estonia

    Marcel Keller

49. Department of Archaeogenetics, Max Planck Institute for the Science of Human History, Jena, Germany

    Marcel Keller, Johannes Krause & Wolfgang Haak

50. SFB 1070 Resource Cultures, University of Tübingen, Tübingen, Germany

    Keiko Kitagawa

51. Department of Early Prehistory and Quaternary Ecology, University of Tübingen, Tübingen, Germany

    Keiko Kitagawa

52. UMR 7194 Muséum National d’Histoire Naturelle, CNRS, UPVD, Paris, France

    Keiko Kitagawa

53. Semenov-Tyan-Shanskii Lipetsk State Pedagogical University, Lipetsk, Russia

    Alexander N. Bessudnov

54. Museum of Natural History, University of Colorado-Boulder, Boulder, CO, USA

    William Taylor

55. Musée d’Anthropologie préhistorique de Monaco, Monaco, Monaco

    Jérome Magail

56. Institute of Archaeology, Mongolian Academy of Sciences, Ulaanbaatar, Mongolia

    Jamiyan-Ombo Gantulga & Turbat Tsagaan

57. Department of Archaeology, Max Planck Institute for the Science of Human History, Jena, Germany

    Jamsranjav Bayarsaikhan

58. Chinggis Khaan Museum, Ulaanbaatar, Mongolia

    Jamsranjav Bayarsaikhan

59. Department of Archaeology, Ulaanbaatar State University, Ulaanbaatar, Mongolia

    Diimaajav Erdenebaatar, Enkhbayar Mijiddorj & Tumur-Ochir Iderkhangai

60. Department of History, Kyrgyz-Turkish Manas University, Bishkek, Kyrgyzstan

    Kubatbeek Tabaldiev

61. Department of Biology, National University of Mongolia, Ulaanbaatar, Mongolia

    Bazartseren Boldgiv

62. Division of Archaeology, Biodiversity Institute, University of Kansas, Lawrence, KS, USA

    Sandra Olsen

63. Institute for Prehistoric and Protohistoric Archaeology, Kiel University, Kiel, Germany

    Cheryl A. Makarewicz

64. ROOTS Excellence Cluster, Kiel University, Kiel, Germany

    Cheryl A. Makarewicz

65. Archaeology of Social Dynamics, Institució Milà i Fontanals d’Humanitats, Consejo Superior de Investigaciones Científicas (IMF-CSIC), Barcelona, Spain

    Silvia Valenzuela Lamas

66. Departament d’Història i Arqueologia–SERP, Universitat de Barcelona, Barcelona, Spain

    Silvia Albizuri Canadell

67. Grup d’Investigació Prehistòrica, Universitat de Lleida, PID2019-110022GB-I00, Lleida, Spain

    Ariadna Nieto Espinet

68. Valencia, Spain

    Ma Pilar Iborra

69. Departamento de Medicina Animal, Facultad de Veterinaria, Universidad de Extremadura, Cáceres, Spain

    Jaime Lira Garrido

70. Centro Mixto UCM-ISCIII de Evolución y Comportamiento Humanos, Madrid, Spain

    Jaime Lira Garrido & Juan Luis Arsuaga

71. Instituto de Arqueología (CSIC–Junta de Extremadura), Mérida, Spain

    Esther Rodríguez González & Sebastián Celestino

72. Laboratori d’Arqueologia Prehistòrica, Universitat Jaume I, Castelló de la Plana, Spain

    Carmen Olària

73. Departamento de Geodinámica, Estratigrafía y Paleontología, Facultad de Ciencias Geológicas, Universidad Complutense de Madrid, Madrid, Spain

    Juan Luis Arsuaga

74. Department of Eneolithic and Bronze Age, Institute of Archaeology National Academy of Sciences of Ukraine, Kyiv, Ukraine

    Nadiia Kotova

75. Department of Archaeology, University of Exeter, Exeter, UK

    Alexander Pryor & Alan Outram

76. Center for the Study of Human Origins, Anthropology Department, New York University, New York, NY, USA

    Pam Crabtree

77. Department of Archaeology, Ethnology and Museology, Al Farabi Kazakh National University, Almaty, Kazakhstan

    Rinat Zhumatayev & Abdesh Toleubaev

78. Scientific Research Department, Orenburg State Pedagogical University, Orenburg, Russia

    Nina L. Morgunova

79. Department of paleontology, Faculty of Geology, Moscow State University, Moscow, Russia

    Tatiana Kuznetsova

80. Institute of Geology and Petroleum Technologies, Kazan Federal University, Kazan, Russia

    Tatiana Kuznetsova

81. Georgian National Museum, Tbilisi, Georgia

    David Lordkipanize

82. Tbilisi State University, Tbilisi, Georgia

    David Lordkipanize

83. Università degli Studi di Milano, Dipartimento di Beni Culturali e Ambientali, Milan, Italy

    Matilde Marzullo, Ornella Prato, Giovanna Bagnasco Gianni & Umberto Tecchiati

84. University of Tehran, Central Laboratory, Bioarchaeology Laboratory, Archaeozoology Section, Tehran, Iran

    Hossein Davoudi & Marjan Mashkour

85. Research Institute and Museum of Anthropology, Lomonosov Moscow State University, Moscow, Russia

    Natalia Ya. Berezina

86. Department of Archaeogenetics, Max Planck Institute for Evolutionary Anthropology, Leipzig, Germany

    Philipp W. Stockhammer, Johannes Krause & Wolfgang Haak

87. Institute for Pre- and Protohistoric Archaeology and Archaeology of the Roman Provinces, Ludwig Maximilian University, Munich, Munich, Germany

    Philipp W. Stockhammer

88. School of Biological Sciences, The University of Adelaide, Adelaide, South Australia, Australia

    Wolfgang Haak

89. Department of Biology, Universidad Autónoma de Madrid, Madrid, Spain

    Arturo Morales-Muñiz

90. Eurasia Department of the German Archaeological Institute, Berlin, Germany

    Norbert Benecke, Sabine Reinhold & Svend Hansen

91. Evolutionary Adaptive Genomics, Institute of Biochemistry and Biology, Faculty of Mathematics and Science, University of Potsdam, Potsdam, Germany

    Michael Hofreiter

92. Department of Evolutionary Genetics, Leibniz-Institute for Zoo and Wildlife Research, Berlin, Germany

    Arne Ludwig

93. Albrecht Daniel Thaer-Institute, Faculty of Life Sciences, Humboldt University Berlin, Berlin, Germany

    Arne Ludwig

94. ArchaeoBioCenter and Institute of Palaeoanatomy, Domestication Research and the History of Veterinary Medicine, LMU Munich, Munich, Germany

    Joris Peters

95. SNSB, State Collection of Anthropology and Palaeoanatomy, Munich, Germany

    Joris Peters

96. ArchaeoZOOlogy in Siberia and Central Asia—ZooSCAn International Research Laboratory, Institute of Archeology and Ethnography of the Siberian Branch of the RAS, Novosibirsk, Russia

    Sergey K. Vasiliev

97. Department of Eastern European and Siberian Archaeology, State Hermitage Museum, St Petersburg, Russia

    Konstantin V. Chugunov

98. Paleoecology Laboratory, Institute of Plant and Animal Ecology, Ural Branch of the Russian Academy of Sciences, Ekaterinburg, Russia

    Natalya A. Plasteeva & Pavel A. Kosintsev

99. Zoological Institute, Russian Academy of Sciences, St Petersburg, Russia

    Gennady F. Baryshnikov & Mikhail Sablin

100. Department of Archaeology, History Faculty, Vilnius University, Vilnius, Lithuania

     Ekaterina Petrova & Elina Ananyevskaya

101. Laboratory for Archaeological Research, Faculty of History and Law, Kostanay State University, Kostanay, Kazakhstan

     Andrey Logvin & Irina Shevnina

102. Department of History and Archaeology, Surgut Governmental University, Surgut, Russia

     Victor Logvin & Saule Kalieva

103. Saryarka Archaeological Institute, Buketov Karaganda University, Karaganda, Kazakhstan

     Valeriy Loman, Igor Kukushkin, Victor Varfolomeyev & Emma Usmanova

104. Toraighyrov University, Joint Research Center for Archeological Studies, Pavlodar, Kazakhstan

     Ilya Merz & Victor Merz

105. Faculty of History, L. N. Gumilev Eurasian National University, Nur-Sultan, Kazakhstan

     Sergazy Sakenov

106. Institute of Archaeology and Steppe Civilizations, Al-Farabi Kazakh National University, Almaty, Kazakhstan

     Viktor Zaibert

107. Department of Anthropology, Alumni Building, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA

     Benjamin Arbuckle

108. Nasledie Cultural Heritage Unit, Stavropol, Russia

     Andrey B. Belinskiy & Alexej Kalmykov

109. Research Center for the Preservation of Cultural Heritage, Saratov, Russia

     Aleksandr I. Yudin

110. Department of Russian History and Archaeology, Samara State University of Social Sciences and Education, Samara, Russia

     Alekandr A. Vybornov, Natalia Roslyakova & Pavel F. Kuznetsov

111. Russian and Foreign History Department, South Ural State University, Chelyabinsk, Russia

     Andrey Epimakhov

112. South Ural Department, Institute of History and Archaeology, Ural Branch of the Russian Academy of Sciences, Ekaterinburg, Russia

     Andrey Epimakhov

113. Archaeological School, Chuvash State Institute of Humanities, Cheboksary, Russia

     Natalia S. Berezina

114. Department of History of the Institute of Humanities, Ural Federal University, Ekaterinburg, Russia

     Pavel A. Kosintsev

115. Department of Human Evolutionary Biology, Harvard University, Cambridge, MA, USA

     David Anthony

116. Anthropology Faculty, Hartwick College, Oneonta, NY, USA

     David Anthony

117. Department of Nordic Studies and Linguistics, University of Copenhagen, Copenhagen, Denmark

     Guus J. Kroonen

118. Leiden University Center for Linguistics, Leiden University, Leiden, The Netherlands

     Guus J. Kroonen

119. Department of Historical Studies, University of Gothenburg, Gothenburg, Sweden

     Kristian Kristiansen

120. Lundbeck Foundation GeoGenetics Centre, Copenhagen, Denmark

     Kristian Kristiansen

121. Department of Biotechnology, Abdul Wali Khan University, Mardan, Pakistan

     Naveed Khan

Contributions

Designed, conceived and coordinated the study: L.O. Provided samples, reagents and material: A. Perdereau, J.-M.A., B.S., A.A.T., A.A.K., S.A., A.H.A., K.A.S.A.-R., T. Seregély, L.K., R.I., O.B.-L., P.B., M.O., J.-C.C., M.B.-M., N.A., M.G., M.M.-d.H., J.W., S.P., A.L.-K., K. Tunia, M.N., E.R., U.S., G. Boeskorov, L.L., R.K., L.P., A. Bălășescu, V. Dumitrașcu, R.D., D.G., V.K., A.S.-N., B.G.M., Z.G., K.S., G.K., E.G., R.B., M.E.A., G.S., V. Dergachev, H.S., N.T., S.G., A. Kasparov, A.E.B., M.A.A., P.A.N., E.Y.P., V.P., G. Brem, B.W., C.S., M.K., K. Kitagawa, A.N.B., A. Bessudnov, W.T., J.M., J.-O.G., J.B., D.E., K. Tabaldiev, E.M., B.B., T.T., M.P., S.O., C.A.M., S.V.L., S.A.C., A.N.E., M.P.I., J.L.G., E.R.G., S.C., C.O., J.L.A., N. Kotova, A. Pryor, P.C., R.Z., A.T., N.L.M., T.K., D.L., M. Marzullo, O.P., G.B.G., U.T., B.C., S.L., H.D., M. Mashkour, N.Y.B., P.W.S., J.K., W.H., A.M.-M., N.B., M.H., A. Ludwig, A.S.G., J.P., K.Y.K., T.-O.I., N.A.B., S.K.V., N.N.S., K.V.C., N.A.P., G.F.B., E.P., M.S., E.A., A. Logvin, I.S., V. Logvin, S.K., V. Loman, I.K., I.M., V.M., S. Sakenov, V.V., E.U., V.Z., B.A., A.B.B., A. Kalmykov, S.R., S.H., A.I.Y., A.A.V., A.E., N.S.B., N.R., P.A.K., P.F.K., D. Anthony, G.J.K., K. Kristiansen, P.W., A.O. and L.O. Performed radiocarbon dating: J.S. Performed wet-lab work: N. Kahn, A. Fages, M.A.K., T. Suchan, L.T.-C., S. Schiavinato, A.F., A. Perdereau, C.G., L.C., A.S.-O and C.D.S., with input from L.O. Analysed genomic data: P.L. and L.O. Analysed uniparental markers: D. Alioglu, with input from P.L. and L.O. Prepared the linguistic index: G.J.K. Interpreted data: P.L. and L.O., with input from B.A., S.R., S.H., D. Anthony, G.J.K., K. Kristiansen and A.O. Wrote the article: L.O., with input from P.L., B.A., S.R., S.H., D. Anthony, G.J.K., K. Kristiansen, A.O. and all co-authors. Wrote the Supplementary Information: P.L., A. Fages, G.J.K. and L.O., with input from all co-authors.

Corresponding author

Correspondence toLudovic Orlando.

Competing interests

The authors declare no competing interests.

Peer review informationNature thanks Samantha Brooks, Fiona Marshall and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visithttp://creativecommons.org/licenses/by/4.0/.

Reprints and permissions