The genomic history of Iberian horses since the last Ice Age

Domain: earth_science
Published: 2025-08-02
Editorial outcome: minor_revision

Abstract
Horses have inhabited Iberia (present-day Spain and Portugal) since the Middle Pleistocene, shaping a complex history in the region. Iberia has been proposed as a potential domestication centre and is renowned for producing world-class bloodlines. Here, we generate genome-wide sequenced data from 87 ancient horse specimens (median coverage = 0.97X) from Iberia and the broader Mediterranean to reconstruct their genetic history over the last ~26,000 years. Here, we report that wild horses of the divergent IBE lineage inhabited Iberia from the Late Pleistocene, while domesticated DOM2 horses, native from the Pontic-Caspian steppes, already arrived ~1850 BCE. Admixture dating suggests breeding practices involving continued wild restocking until at least ~350 BCE, with IBE disappearing shortly after. Patterns of genetic affinity highlight the far-reaching influence of Iberian bloodlines across Europe and north Africa during the Iron Age and Antiquity, with continued impact extending thereafter, particularly during the colonization of the Americas.

Introduction
Horsepower revolutionized human history, arming ancient societies with fast mobility, new ways to make war and farm1. As a prominent symbol of prestige, the horse also deeply impacted past social systems2 by becoming a central element in elite rituals3. While once a matter of intense scientific debates, the domestication homeland has now been identified in the steppes of the lower Don-Volga in the second half of the third millennium Before Common Era (BCE)4. This process involved close reproductive control, shortened generation times, and a severe demographic bottleneck between ~2700 BCE and ~2200 BCE5. This was followed by a rapid expansion across Eurasia ~2200 BCE, initiating the era of large-scale horse-based mobility in human history, ultimately resulting in the near-total replacement of local horse lineages4,5. As a result, all modern domestic horses living on the planet belong to the DOM2 lineage.

Earlier archaeological evidence from Botai, Kazakhstan (~3500 BCE) has suggested horse milking, harnessing and corralling, indicative of horse husbandry6. This interpretation aligns well with the genetic signatures of reproductive control found on site5. The underlying process was, however, almost completely independent genetically from DOM2 domestication, and remained geographically isolated, with Przewalski’s horses as its only surviving descendants7. Today, the vast majority of horses belong to domestic bloodlines or their feral relatives. This stands in stark contrast to the rich diversity of wild horse populations that once roamed Eurasia at the time of domestication4,5,8. When, how, and why these wild lineages disappeared, and whether they contributed to the genetic makeup or characteristics of domesticated horses, remain open questions.

Iberia (present-day Spain and Portugal) provides a striking example of this lost diversity. Recent work has revealed the presence of a genetically-divergent, now-extinct lineage in the region—IBE—morphologically similar to other horses2. IBE carried genetic ancestry from a highly divergent, uncharacterized ghost lineage, particularly evident in its unique Y-chromosome haplogroup2. To date, only six IBE samples have been genetically characterized, spanning the Early Neolithic (Cova Fosca, ~5180 BCE)4, to the Bronze Age (El Acequión, ~1900 BCE)2. Although archaeological bone assemblages throughout the Upper Palaeolithic, Mesolithic and Neolithic periods have suggested favourable climatic conditions for horses in Iberia9–11, the question of whether the area was only inhabited by IBE horses remains unresolved, owing to the limited genetic data available. Key uncertainties also persist regarding (i) the timing of IBE extinction and potential local domestication, possibly as early as the Late Neolithic-Chalcolithic transition12–14, and; (ii) the introduction of DOM2 bloodlines in Iberia, with the earliest genetically characterized DOM2 horse dating to the Early Iron Age (Els Vilars, ~8th–7th centuries BCE)2,4. While the genetic legacy of IBE into modern domestic bloodlines appears minimal, the extent to which past breeding practices were permissive to admixture between DOM2 and local IBE horses remains largely unknown15. Additionally, Iberian horses were highly valued in Antiquity, as attested by Greek and Latin sources describing their military use until the first century BCE and their prominence in Roman chariot racing16,17. During the Iron Age, they served as riding and working animals, as a symbol of prestige with a prominent role in ritual practices, and, to a lesser extent, as a source of meat3. However, how Iberian bloodlines evolved in later periods remains unclear, particularly during the ~800-year-long Muslim rule of al-Andalus18, and the colonization of the New World19.

In this study, we generated genome-wide data from 87 archaeological horse remains to map the genetic diversity of horses in Iberia and across the Mediterranean during the last ~26,800 years. We identified the IBE lineage in Iberian horses from the Upper Palaeolithic, with the latest confirmed individual dating to ~571 BCE (734–408 BCE). The earliest DOM2 Iberian horse in our dataset dated to ~1213–1048 BCE, though the lineage arrived earlier, as evidenced by multiple admixture events with IBE beginning around ~1842 BCE (1973–1711 BCE), and continuing until ~358 BCE (557–159 BCE). This suggests that domesticated horses were not managed in complete isolation from wild populations and that IBE likely became extinct by the last quarter of the first millennium BCE. Finally, DOM2 horses across Europe exhibited strong genetic ties to Iberian populations between ~1200 BCE and ~500 CE, including during Roman times. This underscores Iberia as a major horse-breeding centre throughout the Iron Age and Antiquity, an influence that endured until the colonization of the Americas.

Results
Genomedataset
A total of 87 horse archaeological remains, including petrosal bones, teeth and postcranial elements, were extracted for DNA at the CAGT ancient DNA facilities (Toulouse, France) (Fig. 1; Supplementary Data 1). Low-depth sequencing of triple-indexed DNA libraries on the Illumina MiniSeq instrument revealed 26 females and 61 males showing endogenous DNA content compatible with relatively economical genome characterization by shotgun sequencing (median = 27.21%, range = 2.03–89.75%) (Supplementary Data 1). mapDamage220 base compositional profiles, and PMDtools21 nucleotide mis-incorporation profiles, were aligned with the molecular tools used for generating ancient DNA sequences, including USER enzymatic treatment and DNA library construction22, supporting data authenticity (Supplementary Figs. 1, 2 and 3). A total of 18.46 billion read pairs were produced on Illumina HiSeq4000 and NovaSeq 6000 instruments, of which 4.55 billion could be uniquely aligned with confidence against the horse reference genome, resulting in genome depth-of-coverage values between 0.10-fold and 4.12-fold (median = 0.97-fold; Supplementary Data 1).

The vast majority of the newly characterized genomes originated from Iberia (N = 63) with the remaining 24 originating from specimens excavated in France (N = 16), Italy (N = 3), Morocco (N = 4) and Israel (N = 1), to map possible Iberian genetic affinities and trade networks across the Mediterranean (Fig. 1b, c). The new genomes covered the Late Pleistocene (N = 4; ~25,265–14,179 BCE) to the Modern Period (N = 1; >1656 CE), with a majority, however, belonging to the Iron Age (N = 48), Roman Republic (N = 11), and the Roman Empire/Late Antiquity time periods (N = 5). The newly characterized genomes were combined with a total of 261 ancient horse genomes previously characterized encompassing the whole range of genetic lineages found across Eurasia (Fig. 1a; Supplementary Fig. 4; Supplementary Data 1). They were supplemented with the genomes of 83 modern domestic breeds from around the world, one modern and one historical Przewalski’s horse, and two donkey outgroups, to investigate the genetic landscape of Iberian horses from before and after the Last Glacial Maximum (LGM; 26.5-20/19 kyr scal. Before Present (BP))23, as well as their possible legacy to the present day.

IBE survived in Iberia throughout LGM to the Iron Age
Pseudo-haploidization of sequence data led to the recovery of 7.12 million autosomal transversion positions in at least 80% of specimens considered, including donkey outgroups (median missingness = 7.614%; range = 0.003–92.890%). Principal Component Analysis (PCA), excluding outgroups and divergent Late Pleistocene American horses, revealed that the majority of the newly sequenced specimens (N = 73) co-plotted with DOM2 (D2) (Fig. 2a, Supplementary Figs. 5, 6 and 7). Strikingly, the remaining 14 specimens occupied the most positive and negative scores along Principal Components 1 and 2 (PC1 and PC2), respectively, together with previously identified IBE horses2. PC2 scores discriminated mostly those lineages from the Ural Mountains (UR), Central Asia (BOTAI, BO; BORLY4, BL; NEOBOR, NB; NOVOLINKA, NV; KROKALEV, KR, and; TERSEK, TE), to northeastern Siberia (Equus lenensis, EL) (Supplementary Data 1).

Autosomal Neighbor-Joining phylogenetic reconstruction confirmed that the members of the IBE cluster formed a monophyletic group, receiving maximal (100%) bootstrap support (Fig. 2b; Supplementary Fig. 8). One Middle Bronze Age mare from Moncín (MNP20_Spa_m1558), however, slightly deviated from the other PCA cluster members. It was the first to diverge within the IBE phylogenetic clade, potentially reflecting a more complex genomic make-up.

The autosomal phylogeny also confirmed that several lineages of various genetic affinities once inhabited Europe. IBE showed closest affinities to lineages of Late Pleistocene from southern France (SF), as well as from northern France and Britain (NF). Another related cluster, more divergent from IBE, was also formed by horses from the fourth to third millennium BCE of Denmark (FB), Germany (CWC), Poland (POL), and the Carpathian basin (Czechia, CZ, and; Hungary, HU) (Fig. 2b; Supplementary Fig. 8).

The mitogenome phylogeny displayed, however, a different pattern compared to the autosomal phylogeny, in line with previous reports indicating incomplete lineage sorting and limited phylogeographic resolution for mitochondrial DNA in horses24. Interestingly, the majority of specimens showing autosomal affinity to IBE, including the distinctive MNP20_Spa_m1558 specimen, grouped together within three strongly supported mitochondrial clusters, showing SF or NF members as their closest relatives (Fig. 3a, Supplementary Fig. 9). In contrast, all of the eight newly characterized IBE male specimens grouped together with the three previously sequenced, forming a Y chromosomal clade supported by 100% bootstrap values (Fig. 3b, Supplementary Fig. 10). The latter cluster also included a specimen from Denmark (DJM130x6_Den_m3011; FB), dated to the transition of the fourth-to-third millennium BCE. This specimen was known for carrying one of the most divergent Y chromosomal haplotypes in Eurasia2,4, now linking to IBE.

Combined, our autosomal, mitochondrial and Y-chromosomal phylogenetic analyses reveal that the divergent IBE lineage was genetically related to other wild lineages from western Europe, especially from France (Late Pleistocene). It was present in Iberia from at least ~26.8 kya (Z51_Spa_m24857)25 to ~2.5 kya (CBC11_Spa_m571) (Supplementary Data 1), which considerably extends its temporal timeline, previously documented between ~7.2–3.9 kya2,4. The genetic continuity observed between the Late Pleistocene to Early Iron Age also demonstrates that the IBE lineage survived during the LGM in an Ice Age refugium located south of the Pyrenean mountains, in line with climate niche modelling indicating palaeoenvironmental conditions viable for horses26.

Interestingly, IBE δ13C and δ15N isotopic profiles remained stable through time from the Late Pleistocene to ~2.5 kya, suggesting limited changes in diet and habitat after the LGM. While IBE δ13C values were within the range measured in other horse genetic lineages (Supplementary Fig. 11a), their δ15N values occupied the lower distribution range (Supplementary Fig. 11b). The δ13C and δ15N values in herbivores usually reflect those of the plants they consume, which vary with a range of natural and anthropogenic practices27. In a C3-dominated environment such as Iberia, plant δ13C and δ15N values will primarily vary depending on plant type (with cereals showing higher δ13C values, versus lower δ15N values for legumes), plant parts (with leaves and stems showing higher δ13C and δ15N values relative to grains), and habitat openness (with forested area associated with lower δ13C and δ15N values)28. Climate also plays a dominant role on δ15N values with plants growing in drier exhibiting higher measurements27. Last, domesticated animals consuming protein-rich diet or plants grown on soils fertilized with manure either intentionally, or as a result of repeated pasture use, show increased δ15N values27. The consistent δ15N differences between IBE and DOM2 horses across multiple sites and time periods, combined with the similar δ13C values between lineages, suggest systematic differences in habitat use, forage type, or physiology rather than purely regional effects. The low δ15N values of IBE horses align with those of other Pleistocene horses (such as LPSFR (SF)), supporting continuity in habitat and forage. In contrast, DOM2 horses exhibit increasingly higher δ15N values and a broader δ13C range, potentially due to foraging in open grasslands and increasingly human-controlled diet or habitats. We conclude that the observed differences could result from a combination of (1) distinct habitats (wetter forests for IBE, drier grasslands for DOM2), and; (2) differing diets (lower-protein forage in natural habitats such as legumes and grass for IBE, compared to increasingly anthropogenically-impacted habitats with higher-quality cereals and grasses for DOM2).

Horse genetic affinities in Europe pre-domestication
To further explore the relationship between IBE and other wild horse populations in Europe, we applied AdmixtureBayes29 graph modelling to a selection of seven major genetic horse lineages identified in Europe to date4,5. Lineages from the Carpathian basin and Poland, which cluster together in our autosomal phylogenetic tree (Fig. 2b, Supplementary Fig. 8), and are known to share genetic affinities with the Late Neolithic CWC5 lineage from Germany, were excluded to limit model complexity and computational time, while improving MCMC chain mixing efficiency.

The posterior graph distribution was dominated by a single model, representing 56.5% of the sampled graphs (Fig. 4a), with only two additional models occurring at frequencies above 5% (7.0% and 6.4%, respectively) (Fig. 4b, c). All three models depict a complex history of lineage divergence and admixture, in which Late Pleistocene populations from northern France and UK (LPNFR, NF) emerged as some of the earliest lineages investigated to expand and differentiate into Europe (divergence node n8). The descending populations further differentiated into distinct regional lineages, shaped by contributions from other, more divergent groups. For instance, the three models consistently support the IBE lineage tracking 1.7–6.0% of its ancestry back to a divergent population best modelled as the root and showing some of the largest branch drift estimates (admixture node A1; Fig. 4a–c). The deep branching of this contribution indicates the presence of an ancient horse lineage that had arrived in Europe much earlier, and survived in the Iberian climatic refugium. Furthermore, the Eneolithic population from Romania (ENEOROM, NR) is modelled to primarily derive from a node also forming the URAL (UR) population (n1 on the first two models, and n6 on the third). This node branches deep in the history of lineage divergence (n10), prior to the radiation leading to all other lineages from Europe, indicating yet another expansion into Europe than that forming NF and related descending populations.

In all three models, Late Pleistocene populations from southern France (LPSFR, SF) derived most of their ancestry from the same expansion that led to the formation of NF populations (admixture node A3), with 77.6–80.9% originating to node n9, immediately ancestral to n8. However, the remaining fraction of the SF genetic makeup originated from a population directly related to IBE (n2). As unadmixed members of the IBE lineage were only found in Iberia, thus likely differentiated in Iberia, this suggests that an IBE-related population expanded northward into southern France, where it formed SF by mixing with its local closest geographic neighbours.

The three models differ in the population history underlying the formation of the NR lineage and those occupying Germany ~2700 BCE (CWC) and Denmark ~3000 BCE (FBPWC, FB). The first two models indicate that CWC- and FB-related lineages participated to NR formation, located further east, likely as populations related to both lineages expanded eastward towards Romania and mixed with the divergent local lineage (admixture node A6; Fig. 4a, b). In contrast, the third model supports a reverse scenario, as an NR-related lineage (node n4) appears central in the formation of both CWC and FB lineages. This implies that a population related to the local lineage from Romania (NR) spread westward into Germany and Denmark, where it admixed with the descendants of the earlier expansion that established the NF lineage (Fig. 4c). While further genome-wide time-series are needed to fully resolve the complexity of these expansions and admixture events, the three population graphs reveal a history of significant genetic drift leading to FB and CWC populations. This aligns with early-to-mid Holocene bone assemblages, which indicate that horse populations were small and sparsely distributed across the continent during this period30.

The rise of DOM2 horses in Iberia
PCA and phylogenetic analyses identified specimen MNP20_Spa_m1558 as an outlier member of the IBE lineage. This was confirmed by qpADM31 modelling strongly rejecting MNP20_Spa_m1558 as entirely descending from the single IBE cluster (p-values=0; Supplementary Data 2). Furthermore, ADMIXTURE32 analyses based on the optimal K=2 number of ancestry components identified by cross-validation, indicated that the ancestry component maximized in IBE individuals (and other non-DOM2 lineages) accounted for only 79.7%-80.2% of its genetic ancestry (Supplementary Fig. 12). Considering the same number of components, the Struct-f433 profile of sample MNP20_Spa_m1558 also comprised a higher fraction of the ancestry component maximized in DOM2 horses than other IBE horses (Supplementary Fig. 13). We, thus, explored other possible contributions than IBE, considering all other genetic lineages identified in Europe, as well as DOM2.

qpADM modelling implementing the ‘rotating’ strategy from Harney and colleagues34 testing one possible source at a time while considering all others as outgroup (‘right’) populations, revealed DOM2 as the only compatible source (p-values > 0.01; Supplementary Data 2). Population graph modelling with OrientAGraph35, and including all known possible candidate sources from Europe, also supported the presence of gene flow from the DOM2 lineage into MNP20_Spa_m1558, indicating an admixture proportion of 30.6–34.4% (Supplementary Figs. 14 and 15). These analyses align with D-statistics of the form (EA, IBE; DOM2*, MNP20_Spa_m1558), explicitly testing whether MNP20_Spa_m1558 shares more derived mutations with members of the IBE lineage relative to the earliest DOM2 members known in Europe (DOM2*). The significantly positive D-statistics support an excess of IBE ancestry in MNP20_Spa_m1558 (Fig. 5a; Supplementary Fig. 16), reaching a fraction of 58.91% (range=57.8–60.0%), according to F4-ratios31 (Z-score=36.83; Fig. 5c).

The predominance of IBE horses in Iberia until the second millennium BCE contrasts with their almost disappearance from our dataset during the Iron Age. Only a single of the 20 IBE specimens investigated here (including MNP20_Spa_m1558, plus six previously published) was radiocarbon dated to the first millennium BCE (CBC11_Spa_m571, 734–408 BCE), in contrast to 48 DOM2 horses from Iberia (38 from this study plus 10 previously published, and associated with Iron Age and Roman Republic contexts) (Fig. 1bc; Supplementary Data 1). This indicates a history of genetic turnover following the expansion of DOM2 horses, as shown for other wild lineages across Eurasia4,5.

Importantly, significantly positive D-statistics of the form (EA, IBE; DOM2*, X) indicated an excess of IBE ancestry in only 15 DOM2 specimens present in our dataset (plus MNP20_Spa_m1558; Fig. 5a; Supplementary Fig. 16). Four of these specimens were not from Iberia, revealing indirect genetic influence of IBE outside their native range, potentially through trade and/or exchange of admixed DOM2 horses originally bred in Iberia. These included one sample from the Roman settlement of Lixus, Morocco (LIX01_Mor_329; 255-402 CE), and three samples from Iron Age settlements in Sicily (Motya, MO2_Ita_m450; 475-425 BCE), southeast Britain (Chalk Hill, Chalk3_UK_m474; 541-406 BCE), and southern France (Pech501_Fra_m200; 200 BCE). The latter was excavated at the Pech Maho oppidum in southern Gaul, approximately 130 kilometres north of the Greek colony of Emporion (Empúries, northeastern Iberia). The site has yielded material evidence of Iberian contacts between ~500 and 200 BCE (Supplementary Text). Statistical power was sufficient to estimate that the IBE lineage contributed to 1.22% of the genomic makeup of the admixed specimen from Pech Maho (Z-score=134.72; Fig. 5c), and 1.65–1.90% (Z-scores = 74.04 and 102.61) in two other specimens from Spain (including MR7_Spa_m1131 from the final Bronze Age settlement of El Morredón, which represents the earliest Iberian DOM2 horse in our data set (1213–1048 BCE). While Z-scores were also largely significant in all other 12 admixed individuals tested (69.97 ≤ Z-scores ≤ 130.57), the confidence range estimated for the IBE ancestry could intercept zero, or even extend into negative values. This suggests limited sensitivity and that IBE genetic contribution to the genomes of the remaining horses with non-zero D-statistics was below the minimum successfully estimated value (i.e., 1.22%).

DATES36 analyses for estimating the admixture time between DOM2 and IBE lineages showed significant model fits for all admixed candidates (Fig. 5b), except LIX01_Mor_329 and Chalk3_UK_m474, which were, thus, conservatively considered false-positives (Supplementary Data 3, Fig. 5d). Admixture times were distributed across most of the mid-first millennium BCE, and extended into the first half of the second millennium BCE (Fig. 5b; Supplementary Data 3). The fact that multiple confidence intervals show no overlaps suggests several independent admixture events punctuating the Bronze Age and the Iron Age. These findings support husbandry practices allowing almost continuous contact with wild horses. Remarkably, the earliest admixture date was obtained on the MNP20_Spa_m1558 specimen from Moncín, Spain, and was estimated to ~1842 BCE (range = 1973–1711 BCE). This indicates that DOM2 horses have entered the region in the early second millennium BCE, only a few centuries after their initial domestication in the lower Don-Volga basin (Supplementary Data 3).

D-statistics calculations showed that none of the modern breeds investigated shared an excess of IBE ancestry (Supplementary Fig. 17). This supports previous reports of limited IBE genetic legacy into modern domestic horses2, likely due to a relatively early extinction of this lineage. The absence of measurable IBE-related genetic ancestry in the DOM2 ancient horses present in our dataset after ~358 BCE (midpoint; range = 557–159 BCE), suggests an extinction time by the end of the Iron Age, and likely prior to the Roman period.

The genetic influence of Iberian horses from the Iron Age onwards
Historical sources highlight the exceptional value of Iberian bloodlines from Antiquity to the Modern Period16,17,37–39. To investigate the potential influence of Iberian horse breeding over time, we analyzed ancient DOM2 horses lacking IBE ancestry, covering a period from ~2125 BCE to the modern era (Supplementary Data 1). We used f3-outgroup statistics with donkeys (EA) as an outgroup (Supplementary Fig. 18), to measure genetic affinities between each individual DOM2 horse from Europe and the Mediterranean, and the remaining DOM2 specimens from other archaeological sites. Pairs of samples from the same site or separated by more than 500 years were disregarded. This filter was especially implemented (1) to avoid only capturing within site affinities, (2) to avoid biasing the analysis due to uneven sampling across sites (i.e., those sites with multiple samples would show local affinities that sites represented with a single sample could not), and; (3) to only compare those horses living relatively close in time. This approach allowed us to detect genetic affinities spanning thousands of kilometres (Supplementary Fig. 19). We then mapped the pairs showing the highest detected affinity, visually representing their connection with colour-coded links proportional to the strength of the genetic relationship (Fig. 6).

Prior to ~1200 BCE, the number of unadmixed DOM2 specimens available was scarce across western Europe (Fig. 6a). However, Central European samples most often exhibited stronger genetic affinities with individuals from present-day western Russia, Kazakhstan, and Iran, suggesting that domestic horses in Europe were not yet locally bred but instead sourced from genetically similar stocks originating near the domestication centre.

In contrast, from ~1200 BCE to 500 CE, most horses analyzed—both within and outside Iberia—found their closest genetic relatives in Iberia or southeastern France (Fig. 6b). This likely reflects the emergence of locally bred bloodlines in the region, and their subsequent diffusion through trade networks across western Europe. Notably, Iberian genetic influence extended over vast distances, reaching as far north as present-day Estonia (Ridala) and Britain (Magor), as far east as Anatolia (Acemhöyük), and as far south as Morocco (Lixus) and Tunisia (Althiburos) (Fig. 6b).

Although our dataset limited was incompatible with extensive analyses beyond 500 CE, we detected a substantial expansion of Iberian genetic influence, extending to both Asia (southern Kazakhstan, Zhanaturmus, and; northeastern Iran, Belgheis), and the Americas (Kaw River and Chorillo Grande; Fig. 6c). The affinities with Asia likely reflect the shared genome heritage of Iberian and Oriental bloodlines, previously linked to the spread of horses following the Muslim conquest, which impacted both Europe and Asia2. Meanwhile, the genetic signal in the Americas corresponds to the early introduction of domestic horses by the Spanish and Portuguese colonizers (i.e., the Conquistadors)19.

Methods
Samples and archaeological contexts
All faunal remains in this study were authorized under the appropriate permits for destructive sampling. Sample acronyms (where available) can be found in the Supplementary Information. Furthermore, each sample has a unique laboratory code linking to associated metadata, including the archaeological or palaeontological site of origin. We processed 87 individuals from the Iberian Peninsula, Morocco, France, Italy and Israel spanning a temporal range from Late Pleistocene ~24,857 BCE to the Modern Period (~1800 CE). Our main aim was to clarify the diversity of horse lineages once populating Iberia, and their potential genetic influence in the global history of horse breeding. A detailed description of the archaeological contexts, sites, and individuals is reported in Supplementary Text, and summarized in Supplementary Data 1, together with their associated genetic groups and the various data partitions considered in the different genetic analyses.

Radiocarbon dating and isotopic analysis
A total of 29 new radiocarbon dates were obtained in this study (Supplementary Text). Dating and δ13C and δ15N isotopic measurements were carried out at the Keck Carbon Cycle AMS Laboratory, UC Irvine, following collagen extraction and ultra-filtration from approximately one gram of bone or dental material. Calibration was performed using OxCal Online54 and the IntCal20 calibration curve55.

Ancient DNA extraction and sequencing
The laboratory work was performed at the state-of-art ancient DNA facilities of the Centre for Anthropobiology and Genomics of Toulouse (CAGT, University of Toulouse, France), following the procedures from Fages and colleagues2 and Librado and colleagues4 (Supplementary Text). Petrosal bones were prioritized, then teeth and postcranial elements. The overall methodology involved the powdering of a total of 160–720 mg of bone or dental material using the Mixer Mill MM200 (Retsch) Micro-dismembrator. DNA extraction was carried out to facilitate the recovery of even the shortest DNA fragments56. All but one DNA extract was treated with the USER (NEB) enzymatic mix to eliminate a fraction of post-mortem DNA damage (Supplementary Data 1). Libraries were built from double-stranded DNA templates in which two internal indexes were added during adapter ligation, plus an additional external index during PCR amplification. Libraries were amplified, purified, quantified, and pooled at equimolar concentrations for low-depth sequencing on the Illumina MiniSeq instrument (paired-end mode, 2 × 80). After initial screening for library content, those libraries compatible with an economical characterization of the horse genome through shotgun sequencing were selected for deeper sequencing on the Illumina HiSeq4000 instruments from Genoscope (paired-end mode, 2 × 76; France Génomique), and the Illumina NovaSeq 6000 instruments from Novogene Europe (paired-end mode, 2 × 150).

Read processing and mapping
Demultiplexing, trimming and collapsing of the raw FASTQ read pairs were carried out using AdapterRemoval2 (v2.3.0)57 (Supplementary Text). Collapsed and un-collapsed reads were mapped against the horse reference genome EquCab358, appended with the Y-chromosome contigs from Felkel and colleagues59, and the mitochondrial reference genome (NC_001640)60. The alignments were performed using PALEOMIX bam_pipeline version 1.2.1361 and Bowtie2 (v2.3.4.1)62, with the parameters recommended by Poullet and Orlando63. For nucleotide mis-incorporation and DNA fragmentation profile analyses, we used PMDtools (v0.60)21 and mapDamage2 (v2.2.1)20, considering those alignments against chromosome 1 (Supplementary Figs. 2 and 3), or a random subset of 100,000 alignments (Supplementary Fig. 1), respectively. Overall, we generated sequence data for a total of 87 novel ancient horse specimens, with average depth-of-coverage, removing duplicates and filtering out alignment showing mapping qualities inferior to 25, ranging from 0.10-fold to 4.12-fold (median = 0.97-fold, with 40 above 1X average depth-of-coverage). Information about the number of independent DNA libraries, sequencing efforts, demultiplexing, endogenous content, (mito)genome coverage, and error rates (as calculated by Librado and colleagues4,5) are provided in Supplementary Data 1. The previously published sequence data of 262 ancient horses, 84 modern horses, plus two donkeys (one ancient, one modern), were processed following the same procedures to obtain a comparative genome panel including outgroups and encompassing all known horse genetic lineages (Supplementary Data 1).

SNP variation was identified using ANGSD (v0.940) (htslib 1.16, build on Dec 10 2022 14:56:47), following the pseudo-haploidization procedure from Taylor and colleagues19. This procedure disregarded low-confidence polymorphisms (-baq 0 -remove_bads 1 -uniqueOnly 1 -minMapQ 25 -minQ 30 -rmTrans 1 -rmTriallelic 1e-4 -SNP_pval 1e-16 -C 50), and sites above the 99.6% site-wise coverage threshold observed across all individuals, resulting in a total of 7,122,531 transversion sites, including outgroups.

Phylogenetic reconstructions
Horse mitochondrial DNA haplotypes were obtained from individual BAM alignments manipulation using BCFtools (v1.19)64 (Supplementary Text). Diploid genotype calling was employed to accommodate heteroplasmy, prevalent in horses60. A multifasta file was generated adding the donkey mitochondrial reference haplotype (Accession nb. CM027722.2)65 as an outgroup. We further excluded the highly repeated region located between positions 16,121–16,370 of the horse mitochondrial reference and those samples with over 50% missing data. Maximum Likelihood phylogenetic inference was performed using IQtree (v1.6.12)66 under GTR+F+R10 substitution model, selected via the Akaike Information Criterion (AIC), and a node support evaluated through 1,000 bootstrap pseudo-replicates and ultrafast bootstrap approximation (Fig. 3a, Supplementary Text; Supplementary Fig. 9).

Alignments against Y chromosomal contigs followed the same procedure as autosomal variation identification (with a threshold correspond to 99.5% of the site-wise coverage distribution across samples), however focusing on the single copy (“sc”) regions defined by Felkel and colleagues59, and excluding sites not covered in at least 10% of specimens. This provided data for a total of 14,992 transversion sites across 249 male individuals, including outgroups, which were converted into a multifasta alignment for phylogenetic reconstruction using IQtree (v1.6.12) with TIMe+R3 substitution model based on Akaike Information Criterion (AIC). Node support was assessed as for the mitogenome phylogenetic reconstruction (Fig. 3b, Supplementary Fig. 10; Supplementary Text).

Autosomal phylogeny was generated applying FastMe (v2.1.6.2)67 to the matrix of pairwise genetic distances estimated by Plink v1.9 (--distance square 1-ibs flat-missing)68 from the 7,122,531 transversion sites. Node support evaluated from 100 bootstrap pseudo-replicates using newick utils v1.669 (Fig. 2b, Supplementary Fig. 8; Supplementary Text).

Population structure and gene flow
Principal Component Analysis (PCA) was conducted on pruned autosomal transversions using smartPCA, EIGENSOFT (v6.1.4)70, and restricting to sites present at Minimum Allelic Frequency of 1% (Supplementary Text). Donkey outgroups, two divergent Late Pleistocene horses from north America, as well as sites under linkage disequilibrium (--indep-pairwise 500 10 0.2), were excluded via Plink (v1.9)68, resulting in a final matrix of 3,021,627 transversion sites in 431 individuals (Fig. 2a, Supplementary Figs. 5, 6 and 7).

Pairwise genetic proximity between samples was determined using f3-outgroup statistics31, computed with Calc-f333 using the full matrix of autosomal variation (Supplementary Data 1), plus two donkey outgroups. Standard errors were assessed jackknifing 5 Mb genomic blocks, following Librado and colleagues4. Samples defining individual genetic lineages were grouped together, resulting in N = 36 population groups (Supplementary Data 1), including two donkey outgroups, and N = 311 individual DOM2 horses. f3-outgroup statistics were visualized as heatmap, sorted to reflect the main genetic affinities (Supplementary Fig. 18).

The full matrix of autosomal variation was used for the calculation of D-statistics of the form (EA, IBE; DOM2*, X), where DOM2* represents the earliest DOM2 horses from Europe present in our dataset (N = 8), and X is a sample belonging to the IBE or the DOM2 genetic lineage (Fig. 5a, Supplementary Fig. 16, Supplementary Fig. 17a). We also calculated D-statistics of the form (EA, PR; DOM2*, X), where PR represents Przewalski’s horses (N = 2), and X represents the full list of modern DOM2 horses present in our genome panel (N = 83) (Supplementary Fig. 17b; Supplementary Data 1; Supplementary Text). D-statistics were calculated using qpDstat (v751)31 with default parameters, turning “f4mode” off. Statistical significance was evaluated using the Z-score threshold and Bonferroni corrections for multiple testing. Individuals showing admixed genomic profiles were analyzed using DATES v401036 to estimate the time of admixture between IBE and four pre-defined groups, representing different combinations of DOM2 horses as proxies for early DOM2 sources (Supplementary Data 1 and 3). In these analyses, we considered all autosomes (numchrom=31), and the following parameters: maxdis 1, binsize 0.001, qbin 10, and lovalfit 0.45, turning runfit and afffit options on.

DOM2 and IBE ancestry fractions were estimated calculating F4-ratios in qpF4ratio (version 310)31, assuming blocks of 5 Mb for standard error estimation. The topology tested was (A, B; X; IBE, Outgroup), where A and B represent a subset of N = 5 DOM2 horses from the 3rd and 2nd millennium BCE, and of N = 40 unadmixed Iberian horses from the 1st millennium BCE and early 1st millennium CE, respectively (Supplementary Data 1). The former provided a proxy for early DOM2 horses, while the latter was selected to approximate the Iberian DOM2 gene pool at the time of admixture with IBE.

Genetic ancestry profiling
Unsupervised ADMIXTURE (v1.3.0)32 analysis was carried out using the pruned SNP matrix (MAF=1%) to investigate the genetic profile of the horse specimens present in our dataset (Supplementary Fig. 12a, b), including or not the two divergent Late Pleistocene horses from north America (LPNAMR, NAM; Supplementary Fig. 12c, d; Supplementary Data 1). The optimal number of K=2 ancestry components was assessed through cross-validation. Confidence intervals for the proportions of genetic ancestries were determined through 100 bootstrap pseudo-replicates (-B 100 -c 0.0001).

Additionally, f4-statistics were computed using Calc-f433, using the same procedure as for the calculation of f3-Outgroup statistics, but considering more population groups to minimize permutations and running time (Supplementary Data 1). Unadmixed modern DOM2 horses were grouped together, as well as those ancient DOM2 individuals from the same archaeological site, when showing similar PCA placement and ADMIXTURE32 profiles. A total of K=2 ancestry components were assumed, in line with the optimal number of ancestry components identified by cross-validation in ADMIXTURE32 analyses (Supplementary Fig. 13; Supplementary Text).

Admixture modelling and population graphs
We used qpADM v81031 to model the genomic make-up of sample MNP20_Spa_m1558 (IBEm, identified as a genetic outlier to the IBE genetic cluster), considering various combinations of outgroup (‘right’) populations and potential source (‘left’) populations encompassing the entire range of horse genetic lineages present in Europe (Supplementary Data 1 and 2). To test for possible additional sources than the IBE lineage (Supplementary Data 2), we applied the procedure from Harney and colleagues34, ‘rotating’ one of the sources at a time. Analyses turned off the allsnps option, and significance was set at 0.01.

We also used OrientAGraph (version 1.0)35 to test for possible sources for the complex genomic makeup of sample MNP20_Spa_m1558 (IBEm) (Supplementary Figs. 14 and 15). This analysis included a total of N = 9 population groups from Europe, including CWC, IBE, FBPWC (FB), HUNG (HU), LPNFR (NF), LPSFR (SF), ENEOROM (NR), ENEOCZE (CZ), and Vert311, plus EA as outgroup, and a group of unadmixed ancient DOM2 horses from Iberia (N = 40). The analysis was repeated considering LPNAMR (NAM) as another outgroup, considering up to M=3 migration edges. A total of 100 SNPs per block (-k) were used for estimation of the covariance matrix, and a round of global rearrangements was carried out after adding all populations (-global).

Population graphs were reconstructed with AdmixtureBayes29 considering 32 parallel MCMC chains for mixing (--MCMC_chains), while limiting the number of possible admixture events to 10 (--max_admixes). The first 50% of 250 million iterations were used as burnin, and the remaining were thinned every 250. In this analysis, the full autosomal TPED matrix was restricted, using Plink v1.968, to seven population groups representing major genetic horse lineages in Europe (IBE, SF, NF, CWC, FB, UR and NR), along with two donkey specimens as outgroups (EA; Fig. 4; Supplementary Data 1; Supplementary Text).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, an ancient genomic analysis of 84 ancient horse samples demonstrates the genetic history of horses in\nIberia. The replacement of the IBE horse lineage with the DOM2 lineage and the global spread of Iberian horses is well-\nrepresented and very interesting. This manuscript analyzes many ancient genomes with very low coverage, but it needs\nto be longer to explain the methods that made this possible. Also, the authors should explain what percentage of SNPs they\nhave imputed if they imputed. Without these explanations in detail, it is difficult to know whether the analysis of an ancient\ngenome with an average of 0.95x is reliable or not. The figures and tables could be easier to read and need to be improved.\nMajor comment\nOn average, 0.95x coverage genomic data is used for the analysis, but the methods are not described in detail when I read\nthe methods in the Supplementary information. It should be more detailed, not just include citations. Since this analysis\nmethod is one of the main topics of this study, the extent of imputation and the percentage of imputed SNPs should be\nprovided. The authors should also state whether this analysis leads to any bias and, if so, how they addressed it.\nThe “Counts were further processed to randomly sample one nucleotide per position in every single sample” on line 850 of\nthe Supplementary information method, Does this randomly sample not affect the results?\nThis method yielded 4 million SNPs in transversion sites alone. Remarkably, such a large number of sites were obtained\nfrom an average of 0.95x ancient genomes.\nFor filtering SNPs it states “(-SNP_pval 1e-16 -remove_bads 1 -minMapQ 25 -minQ 30 -rmTrans 1 -rmTriallelic 1e-4 -C 50)”.\nI understand what these options mean, but what about depth?\nWithout these explanations in detail, it is difficult to know whether the analysis of an ancient genome with an average of\n0.95x is reliable or not.\nMinor comment\nLine 120: “84 ancient horse specimens”, please add the range of coverage and the median coverage of ancient genomes in\nabstract.\nLine: 122: add explanation of IBE.\nIntroduction: please explain a little more about DOM2.\nLine 167, 171", "reviewer_174": "do not show results in introduction section (i.e. Supplementary Table 1, Supplementary Text, Figs. 1a,\n1b and 1c).\nLine 190: please show the results of mapdamage2\nLine 215: mark “DOM2” in Fig. 2a. It was difficult to find.\nLine 218: which plot is IBE in Fig. 2a? Difficult to find.\nLine 222: mark IBE cluster in Fig. 2b. Difficult to find.\nLine 222: show the bootstrap values.\nLine 234 and 237: Fig. 3a and 3b are hard to see. Please revise to show what you want to show.\nLine 242-243: “divergent IBE lineage showed broad genetic affinities”. Why you did not use outgroup f3 statistics?\nLine 264-265: “the impossibility to model this specimen as belonging to a unique IBE source with qpADM”. What does it\nmean? MNP20_Spa_m1558 is not a source for IBE?\nLine 266: which program do you use for D-stat? What is the tree topology?\nLine 268: please mark MNP20_Spa_m1558 in Supp. Fig. 6.\nLine 280: why did you remove a possibility of DOM2 from “another of the wild horse lineages”?\nLine 283: please explain “five out of six secondary sources”.", "reviewer_2": "This is a nice study based on extensive work, and Nature Communications seems a good fit. As a general comment, I think\nthere is more that could be said to make the broader relevance of the study clear, for example about the process of\ndomestication and what the history of equid husbandry can tell us about human history, especially in Europe. These topics\nare mentioned briefly, but much of the text may be of greater interest to those who study this system in particular. I have a\nnumber of minor suggestions and comments (below).\nLines 118-119 – does this imply a horse exporting area (is that the implication of breeding)? Briefly say here why this would\nbe important to know.\nLine 122 – IBE is a genetic lineage – say so briefly here\nLine 135 – a map to clarify this part of the introduction (showing place names relevant to the study - supplemental perhaps)\nwould be useful.\nLines 136-137 – give reference for papers defining these genetic lineages.\nLine 147 – don’t need ‘the region of’\nLines 152-154 – should say more about the data that shows that the IBE lineage was from the wild populations, when first\ndetected, how mixing was determined, etc. In a general science journal, many readers won’t be familiar with this system and\nstudy specific terms – need to bring the general reader on board a bit.\nLine 157 – some further information on the DOM2 lineage would help at this point.\nLines 165-180 – seems to be quite a lot of results in the introduction.\nLines 184-195 – probably OK, but mixes methods and results.\nLines 212-213 – it is sensible to focus on transversions for these data, but the strategy should be briefly explained. Also,\nsaying that ‘at least 80%’ of specimens had 3.95 million identified positions begs the question of what the range of coverage\nwas for the other 20%.\nLines 215-216 – presumably the circle for published and triangle indicated in Figure 1 follows through for Figure 2, but\nshould be mentioned in the Figure 2 caption. Also say how many samples are collapsed for D2 (shown expanded in the\nsupplement).\nLine 221 – this is the first mention of this being an autosomal phylogeny – were sex chromosomes excluded from the\nreference? Probably it is to distinguish from mitochondrial, but if so, start Line 232 that way (‘For the mitogenome\nphylogeny…’).\nLines 232-234 – in the mtDNA phylogeny, IBE falls out in multiple lineages – worth mentioning here? What is the\nsignificance, since they cluster together in the nuclear tree?\nLine 234 – presumably the color-coding in Figure 3 matches that in Figure 2, but this should be indicated in the caption.\nLine 238 – ‘date’ should be dated\nLines 240-241 – DJM130x6_Den_ m3011 seems to be nearly the same as some IBE haplotypes in the Y phylogeny…\nLine 243 – presumably the wild lineages are identified by divergence from D2, and there are none that are extant – is that\nright? Anything else that identifies them as ‘wild’?\nLines 247-250 – why would this apply to the IBE lineage in particular – wouldn’t it be the same for any horse lineage?\nLines 254-257 – could this be associated with regional variation in suitable forage?\nLine 264- 265 – don’t follow – ‘impossible to model’, but assigned to IBE with qpADM with high significance?\nLine 268 – what are the ESSI and NAM (domestic from north America?) that stand out in these Admixture plots?\nLines 275-278 – not when K=2 (6b)\nLines 300-301 – seven out of how many, and how does that choice affect the graph?\nLines 303-304 – don’t quite follow this – in the consensus tree the IBE side of the lineage split seems to lead directly to all\nothers except UR, but this text seems to imply that only the other side of the split leads broadly to others.\nLines 322-323 – Figure 1 is a bit confusing with respect to the distribution of D2, which appears only in North America in\nfig1a, but then in Europe near Anatolia in fig1b. Neither seems a good fit with this text.\nLines 326-334 – this would be clearer if the map figure in fig 5d was clearer and referred to in the text. The caption notes 6\nD2 individuals in the map, but there are more than 6 individuals illustrated there. It isn’t easy to follow which are meant to be\nD2 and which not… A figure, perhaps as supplement, that tracked the discussion in the text would help.\nLines 360-361 – not sure I get ‘if any’, as there does seem to be evidence of admixture in supp fig 9. The point about\nconfounding from admixture PR could be much more clearly illustrated in the figure, and noted in the text.\nLine 371 – supplementary fig 10 is very difficult to interpret, certainly at normal size. At least make the key big enough to\nread, but also provide more information in the caption.\nLines 375-376 – figure 6 is very difficult to interpret. The x’s are essentially invisible, the color coding very hard to\ndiscriminate, and while open circles and x’s are defined, closed circles are not. Does the color indicate ancestry, while the\nposition where they were found is shown on the map? Not clear from the caption.\nLines 387-389 – not sure if this meant to mean domestic lineages from Iberia, or the IBE lineage descendants.\nLine 401 – strong rather than drastic?\nLines 402-404 – were there any other lines of wild horses identified in Europe, or just IBE so far?\nLine 405 – where might the origin of the IBE lineage be from? Say if that still isn’t know, other than ‘outside Iberia’.\nLines 419-423 – this is interesting, and I think could be one of the main contributions of this study (shifts in domestication\npractice), but I find the sentence hard to follow, and feel like it could be expanded on to consider that question more fully.\nLines 430-431 – but involving domestic lineages, after the IBE lineage had gone extinct – correct?\nLines 434-438 – why would this be important/ useful to know in detail?\nLine 438 – mobility in central…\nLine 441 – associated with\nLine 444 – ‘investigating’ instead on ‘entering into’?\nLine 500 – please give details on quality filters applied.\nLines 530-533 – need more details on the construction of the nuclear phylogeny especially (seems to be even less detail in\nthe supplement than here), such as partitioning, evolutionary model, strategy (transversions only), etc."}, "author_response": "Point-by-point answers to the reviewers’ requests:\nReferees' comments:\nReviewer #1 (Remarks to the Author):\nPOINT R1.1. In this manuscript, an ancient genomic analysis of 84 ancient horse samples\ndemonstrates the genetic history of horses in Iberia. The replacement of the IBE horse lineage\nwith the DOM2 lineage and the global spread of Iberian horses is well-represented and very\ninteresting. This manuscript analyzes many ancient genomes with very low coverage, but it needs\nto be longer to explain the methods that made this possible. Also, the authors should explain what\npercentage of SNPs they have imputed if they imputed. Without these explanations in detail, it is\ndifficult to know whether the analysis of an ancient genome with an average of 0.95x is reliable\nor not. The figures and tables could be easier to read and need to be improved.\nANSWER R1.1:\nWe appreciate the constructive suggestions provided by reviewer #1. We have now expanded the\ndescription of our methods, which follow gold standards in ancient DNA analyses (see Orlando\net al. Nat Rev Methods Primers 2021 for a recent review). Briefly, post-mortem DNA damage\nconsiderably reduces the quality and quantity of DNA molecules available for sequencing. As a\nresult, genome coverage-of-depth remains generally limited (median=0.97-fold, now that 3\nadditional genomes have been integrated for the revisions, leaving a final data set of 87 newly-\ncharacterized ancient genomes). As standard practice with such data, genetic variation is assessed\nfrom random sampling of one read at those genome locations covered multiple times. This creates\na so-called ‘pseudohaploid’ data set, in which individuals are homozygous for the sample\nnucleotidic base sampled. This approach has the advantage to mitigate differential drift across\npopulations (as one of the two possible alleles is drawn randomly at heterozygous sites). However,\nno genotypes are called, and analyses relying on the distribution of particular genotype classes\n(eg runs-of-homozygosity) cannot be carried out, unless genotypes are statistically imputed. Our\nprevious work demonstrated the accuracy of statistical imputation for SNPs with variants\nsegregating at minimum 5% frequency in modern domestic horses (Todd et al. iScience 2023).\nHowever, the performance drastically drops when applied to other genomic background than the\nonly two lineages represented in modern populations around the world (ie DOM2, and PRZW,\nthe Przewalski’s horse). This is due to the absence of divergent, now-extinct lineages in\nimputation panels. Since one of the objectives of our study was to track the extinction timing of\none divergent such lineage, imputation was not an option. Our analyses thus fully rely on the\nstandard pseudo-haploidisation procedure described above. To further assess the authenticity and\nquality of our data, we now present two new analyses, aimed at checking the presence of\nmolecular signatures of post-mortem damage (Supplementary Figs. 1-3), and assessing error rates\nin each individual genome (these are comparable, and even improved, relative to those of ancient\ngenomes; Supplementary Table 1). The expanded methodological section now reads as follows\n(page 19, lines 534-557):\n“Demultiplexing, trimming and collapsing of the raw FASTQ read pairs were carried out using\nAdapterRemoval2 (v2.3.0) (Schubert et al. 2016) (Supplementary Information). Collapsed and\nun-collapsed reads were mapped against the horse reference genome EquCab3 (Kalbfleisch et\nal., 2018), appended with the Y-chromosome contigs from Felkel et al. (2019), and the\nmitochondrial reference genome (NC_001640) (Xu & Arnason 1994). The alignments were\nperformed using PALEOMIX bam_pipeline version 1.2.13 (Schubert et al. 2014) and Bowtie2\n(v2.3.4.1) (Langmead et al. 2012), with the parameters recommended by Poullet and Orlando\n(2020). For nucleotide mis-incorporation and DNA fragmentation profile analyses, we used\nPMDtoolsv0.60 (Skoglund et al. 2014) and mapDamage2 (v2.2.1) (Jónsson et al., 2013),\nconsidering those alignments against chromosome 1 (Supplementary Figs. 2 and 3), or a\nrandom subset of 100,000 alignments (Supplementary Fig. 1), respectively. Overall, we\ngenerated sequence data for a total of 87 novel ancient horse specimens, with average depth-of-\ncoverage, removing duplicates and filtering out alignment showing mapping qualities inferior\n1\nto 25, ranging from 0.10-fold to 4.12-fold (median=0.97-fold, with 40 above 1X average depth-\nof-coverage). Information about the number of independent DNA libraries, sequencing efforts,\ndemultiplexing, endogenous content, (mito)genome coverage, and error rates (as calculated by\nLibrado and colleagues (2021); Librado and colleagues (2024)) are provided in Supplementary\nTable 1. The previously published sequence data of 262 ancient horses, 84 modern horses, plus\ntwo donkeys (one ancient, one modern), were processed following the same procedures to\nobtain a comparative genome panel including outgroups and encompassing all known horse\ngenetic lineages (Supplementary Table 1). SNP variation was identified using ANGSD (htslib\n1.16, build on Dec 10 2022 14:56:47), following the pseudo-haploidisation procedure from\nTaylor and colleagues (2023). This procedure disregarded low-confidence polymorphisms (-\nbaq 0 -remove_bads 1 -uniqueOnly 1 -minMapQ 25 -minQ 30 -rmTrans 1 -rmTriallelic 1e-4 -\nSNP_pval 1e-16 -C 50), and sites above the 99.6% site-wise coverage threshold observed\nacross all individuals, resulting in a total of 7,122,531 transversion sites, including outgroups.”\nWe naturally provide further detail in the Supplementary Material, and fully report the distribution\nof missing data across the newly-characterized data, and for the whole dataset (Supplementary\nMaterial, pages 30-31, lines 811-830):\n“Individual files of horses and donkey outgroups were processed in ANGSD (v0.940) (htslib:\n1.16, build Dec 10 2022 14:56:47 (Korneliussen et al. 2014) to retrieve counts of each\nnucleotide at those transversion positions that are covered in at least 80% of the specimens and\nare likely polymorphic in the overall genome sequence panel considered. Sequence filters\nincluded: (1) a minimum base quality equal to 30 (-minQ 30), (2) a minimum mapping quality\nequal to 25 (-minMapQ 25), (3) the removal of transitions (-rmTrans 1), to limit the impact of\nsequencing errors at sites that have been damaged post-mortem (following Green and\ncolleagues 2010), (4) the removal of reads flagged as bad (-remove_bads 1), (5) the down-\nscaling of the mapping quality of those reads showing excessive mismatches (-C 50), (6) the\nremoval of tri-allelic sites (i.e. with a p-value of being triallelic of 1e-4; -rmTriallelic 1e-4), (7)\nthe removal of sites unlikely to be polymorphic (i.e. showing a p-value for being a SNP lower\nthan 1e-16; -SNP_pval 1e-16), (8) -baq 0 (no Base Alignment Quality–BAQ– correction), and;\n(9) -uniqueOnly 1 (remove reads with multiple best hits). Counts were further processed to\nrandomly sample one nucleotide per position in every single sample, disregarding positions\nwithin the 99.6% quantile of the coverage distribution. This procedure returned a final TPED\nmatrix including 7,122,531 pseudo-haploid autosomal transversion positions in 433 horses, and\n2 donkeys. The fraction of missing data was between 1.340% (Pech1053_Fra_m200) and\n90.590% (GAR37_Spa_m14364) (median = 31.820%) for the newly characterized ancient\nspecimens, and between 0.003% (LR18x62_Rus_m1849) and 92.890% (NW22_USA_1662)\n(median = 7.614%) across the whole dataset (Missingness; Supplementary Table 1).”\nFinally, we have improved the readability of the different figures, as requested by both\nreviewers (see ANSWERS R1.10, R1.14 or R2.4 as examples).\nMajor comment\nPOINT R1.2. On average, 0.95x coverage genomic data is used for the analysis, but the methods\nare not described in detail when I read the methods in the Supplementary information. It should\nbe more detailed, not just include citations. Since this analysis method is one of the main topics\nof this study, the extent of imputation and the percentage of imputed SNPs should be provided.\nThe authors should also state whether this analysis leads to any bias and, if so, how they addressed\nit.\nANSWER R1.2:\nSee ANSWER R1.1. Pseudo-haploidizing data is a gold-standard procedure to handle low-\ncoverage ancient genomes, and is common to >99% of the ancient DNA literature. Previous work\nhas established that reliable ancestry information can be recovered from analyses of such data\n2\nwith coverages as low as 0.1-fold, such as through PCA, genetic clustering, and ancestry\ndeconvolution (Malaspinas et al. Bioinformatics 2014). These data are also suitable for testing\ngene flow and estimating admixture proportions and underlying timing through D-statistics\n(Patterson et al. Genetics 2012), F4-ratios (Patterson et al. Genetics 2012), and DATES\n(Chintalapati et al. 2022), representing the basic methodology underlying our work. It is\nnoteworthy that this is the very same methodology that has allowed our group to identify the\ngenetic homeland of domestic horses (Gaunitz et al. Science 2018, Librado et al. Nature 2021,\nLibrado et al. Nature 2024), and clarify the genetic affinities of major lineages throughout the last\n4,000 years (Fages et al. Cell 2019), including the colonial period (Taylor et al. Science 2023).\nTypically, these studies were based on data of similar quality (coverage) than the one we are\npresenting here. Supplementary Table 1 provides complete information pertaining to sequencing\nfor each of the individuals sequenced in this study (eg. the number of libraries, mitochondrial and\nnuclear coverage, endogenous content, or collapsed and uncollapsed reads), error rates,\nmissingness, and the analytical methods are fully described in the “Read Processing and mapping”\nsection from Supplementary Information, which reads as described in ANSWER R1.1.\nPOINT R1.3. The “Counts were further processed to randomly sample one nucleotide per\nposition in every single sample” on line 850 of the Supplementary information method, Does this\nrandomly sample not affect the results?\nThis method yielded 4 million SNPs in transversion sites alone. Remarkably, such a large number\nof sites were obtained from an average of 0.95x ancient genomes.\nANSWER R1.3:\nThe reason why approximately 4M transversion SNPs were identified is due to the presence of\ndonkey outgroups in the data set. Donkeys and horses show roughly 1% sequence divergence\n(Jonsson et al. PNAS 2014), versus 0.12% amongst modern horses (Schaefer et al. BMC\nGenomics 2017). Therefore, the vast majority of the variation is represented by sites\ndifferentiating the 4 million years-old divergence between and its closest relatives. Analyses such\nas phylogenetic reconstruction, F4-ratios, f3-outgroup and D statistics require the presence of an\noutgroup. Analyses such as PCA do not, and removing the donkey outgroups (plus the most\ndivergent lineage of Late Pleistocene American horses) left a total of 3,021,627 SNPs, as indicated\nin the Methods section (we now used the less restrictive Minimum Allelic Frequency filter of 1%,\nto identify potential pattern of genetic sharedness that are less common). The random sampling\nprocedure is a standard procedure in ancient DNA analysis, preventing bias that may arise from\ndifferential drift (and inbreeding) levels amongst individuals (Rizzi et al. Genet Sel Evol 2012).\nThe conditioning on transversions only is also a standard procedure in ancient DNA analysis,\nlimiting the introduction of sequencing errors resulting from the post-mortem DNA decay\nprocess, especially at cytosine residues. Such sites experience fast deamination post-mortem\n(Briggs et al. PNAS 2007), resulting in the introduction of C→T transitions, and their\ncomplementary G→A transitions. Therefore, filtering transitions, as proposed first by Green and\ncolleagues when reporting the first Neanderthal genome (Science 2010), is essential. Another\nimportant element of our procedure is to disregard sites showing exceptional coverage across all\nsamples. This was defined originally by the same authors, and implemented by our group when\ndeveloping the automated Paleomix pipeline used in this study (Schubert et al. Nat Protocols\n2014), to mitigate possible structural variation not accounted by the reference genome, resulting\nin the accumulation of reads from multicopy regions into the single copy region represented in\nthe genome reference assembly. It is noteworthy that the population structure identified in the\npresent study recapitulates that characterized over the last decade of ancient DNA research on\nhorses, building confidence in the variation identified in the present study.\nPOINT R1.4. For filtering SNPs it states “(-SNP_pval 1e-16 -remove_bads 1 -minMapQ 25 -\nminQ 30 -rmTrans 1 -rmTriallelic 1e-4 -C 50)”. I understand what these options mean, but what\nabout depth?\ngenome with an average of 0.95x is reliable or not.\nANSWER R1.4:\nThe rationale underlying filtering parameters is now explicated in the Supplementary Information,\nas follows (Supplementary Material, page 30, lines 814-825):\n“Sequence filters included: (1) a minimum base quality equal to 30 (-minQ 30), (2) a minimum\nmapping quality equal to 25 (-minMapQ 25), (3) the removal of transitions (-rmTrans 1), to limit\nthe impact of sequencing errors at sites that have been damaged post-mortem (following Green\nand colleagues 2010), (4) the removal of reads flagged as bad (-remove_bads 1), (5) the down-\nscaling of the mapping quality of those reads showing excessive mismatches (-C 50), (6) the\nremoval of tri-allelic sites (i.e. with a p-value of being triallelic of 1e-4; -rmTriallelic 1e-4), (7)\nthe removal of sites unlikely to be polymorphic (i.e. showing a p-value for being a SNP lower than\n1e-16; -SNP_pval 1e-16), (8) -baq 0 (no Base Alignment Quality–BAQ– correction), and; (9) -\nuniqueOnly 1 (remove reads with multiple best hits). Counts were further processed to randomly\nsample one nucleotide per position in every single sample, disregarding positions within the\n99.6% quantile of the coverage distribution.”\nGiven the random sampling procedure implemented generating pseudo-haploid data, no depth\nfilters are requested (ie we don’t attempt any genotyping, and no heterozygous sites are identified,\nsince the single allele drawn is conventionally considered homozygous). However, we typically\ndisregard those outlier sites in the site-wise coverage distribution across samples, as these may\nrepresent unidentified paralogs that are represented as single copies in the reference genome\nassembly.\nMinor comment\nPOINT R1.5. Line 120: “84 ancient horse specimens”, please add the range of coverage and the\nmedian coverage of ancient genomes in abstract.\nANSWER R1.5:\nWe have added the information suggested by reviewer #1 in the first Results section entitled\n‘Genome data set’ (page 6, lines 204-207):\n“A total of 18.46 billion read pairs were produced on Illumina HiSeq4000 and NovaSeq 6000\ninstruments, of which 4.55 billion could be uniquely aligned with confidence against the horse\nreference genome, resulting in genome depth-of-coverage values between 0.10-fold and 4.12-fold\n(median=0.97-fold; Supplementary Table 1).”\nPOINT R1.6. Line: 122: add explanation of IBE.\nANSWER R1.6:\nWe agree and have extended the Introduction to provide more contexts on the IBE genetic lineage.\nIt reads as follows (pages 4-5, lines 158-167):\n“Iberia (present-day Spain and Portugal) provides a striking example of this lost diversity. Recent\nwork has revealed the presence of a genetically-divergent, now-extinct lineage in the region —\nIBE— morphologically similar to other horses (Fages et al. 2019). IBE carried genetic ancestry\nfrom a highly divergent, uncharacterized ghost lineage, particularly evident in its unique Y-\nchromosome haplogroup (Fages et al. 2019). To date, only six IBE samples have been genetically\ncharacterized, spanning the Early Neolithic (Cova Fosca, ~5,180 BCE) (Librado et al. 2021), to\nthe Bronze Age (El Acequión, ~1900 BCE) (Fages et al. 2019). Although archaeological bone\nassemblages throughout the Upper Palaeolithic, Mesolithic and Neolithic periods have suggested\n4\nfavourable climatic conditions for horses in Iberia (Bernáldez-Sánchez et al. 2024; Sánchez-\nHernández et al. 2020; Sanz-Royo et al. 2020), the question of whether the area was only\ninhabited by IBE horses remains unresolved, owing to the limited genetic data available. Key\nuncertainties also persist regarding…”\nPOINT R1.7. Introduction: please explain a little more about DOM2.\nANSWER R1.7:\nWe agree and have developed the introduction to reflect current state of knowledge regarding the\nDOM2 lineage. It reads as follows (page 4, lines 137-146):\n“Horse power revolutionized human history, arming ancient societies with fast mobility, new\nways to make war and farm (Kelekna 2009). As a prominent symbol of prestige, the horse also\ndeeply impacted past social systems (Fages et al. 2019) by becoming a central element in elite\nrituals (Iborra Eres et al. 2023). While once a matter of intense scientific debates, the\ndomestication homeland has now been identified in the steppes of the lower Don-Volga in the\nsecond half of the third millennium Before Common Era (BCE) (Librado et al. 2021). This\nprocess involved close reproductive control, shortened generation times, and a severe\ndemographic bottleneck between ~2700 BCE and ~2,200 BCE (Librado et al. 2024). This was\nfollowed by a rapid expansion across Eurasia ~2200 BCE, initiating the era of large-scale\nhorse-based mobility in human history, ultimately resulting in the near-total replacement of\nlocal horse lineages (Librado et al. 2021; Librado et al 2024). As a result, all modern domestic\nhorses living on the planet belong to the DOM2 lineage.”\nPOINT R1.8. Line 167, 171"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Garrido et al. reported an ancient genomic analysis of 84 ancient horse samples, which demonstrates the genetic history of\nhorses in Iberia. The replacement of the IBE horse lineage with the DOM2 lineage and the global spread of Iberian horses\nare well-represented and very interesting. However, a paragraph titled “Horse genetic affinities in Europe pre-domestication”\nis ambiguous for me. There are many other points that should be revisited, and they are listed below. Please revise these\npoints.\nLine 129, “genetic evidence”: what kind of analysis showed this evidence?\nLine 242-245: too small to find southern France (SF), northern France and Britain (NF), Denmark (FB), Germany (CWC),\nPoland (POL), and the Carpathian basin\n(Czechia, CZ, and; Hungary, HU) in Fig. 2b and Supplementary Fig. 8. Also, the OUT names in Fig. 3 are too small to find.\nLine 258, which tree is “our phylogenetic analyses”?\nLine 259, can “genetic affinities” be mentioned from a phylogenetic tree?\nFig. 4 model to biogeographic DISCUSSION, but the connection between the model and biogeography is unclear.\nLine 305, I cannot find “19.1–22.4%” from Fig. 4.\nLine 306, I am not sure what the evidence of “IBE differentiated in Iberia” is.\nLine 309, what is “a deeply branching lineage”?\nLine 318-319, what does “a significant genetic contribution to the IBE lineage, with 1.7–6.0% of its ancestry tracking back to a highly drifted” mean? How can I find the result from Fig. 4?\nLine 346, it is better to explain “(EA, IBE; DOM2*, MNP20_Spa_m1558)” from general readers. Line 358: (EA, IBE; DOM2*,\nX).\nLine 369, “Statistical power was sufficient” is not sufficient. Please show Z scores.\nLine 375, why 1.22%? Why did you remove 1.65%-1.90%?\nWhat is “density” in Supplementary Fig. 19? F3 value?\nFig. 6: Why are close-distance samples not the highest in f3 value? The f3 value from close-distance samples, for example,\ntwo samples from the same site, should be the highest.", "reviewer_2": "This study is a continuation of and dependent on previous work on horse lineages, largely by the same group of\nresearchers. As such there was a tendency to assume knowledge of the previous work throughout. The revision does a\ngood job of addressing this by adding text that elaborates on previous knowledge and explains their interpretations more\nfully. This has extended the length of the paper, but in a worthwhile way in my opinion. The extent of the data, and the\nchallenges with ancient DNA, make the text, interpretation and conclusions complex, raising questions for the reader (some\nshared by myself and the other reviewer). I feel that the authors have made a thorough attempt to address these questions,\nand where appropriate provided useful caveats and cautionary statements. One aspect I feel they could take a bit further is\ntheir description and justification of their phylogenetic strategies. They note that they had settled on a strategy for their\nphylogenetic approach in earlier studies, but there needs to be some text referring both to the earlier studies (now provided)\nand a brief overview of their reasons for adopting the relatively unusual strategies they use. There are also still some minor\nissues with illustrations (e.g. figure 18 key size) though in general the illustrations are much improved. Overall, I think the\nauthors have addressed my questions and that their paper is now ready for publication."}, "author_response": "Point-by-point answers to the reviewers’ requests:\nReferees' comments:\nReviewer #1 (Remarks to the Author):\nGarrido et al. reported an ancient genomic analysis of 84 ancient horse samples, which\ndemonstrates the genetic history of horses in Iberia. The replacement of the IBE horse lineage\nwith the DOM2 lineage and the global spread of Iberian horses are well-represented and very\ninteresting. However, a paragraph titled “Horse genetic affinities in Europe pre-domestication” is\nambiguous for me. There are many other points that should be revisited, and they are listed below.\nPlease revise these points.\nPOINT 1.1:\nLine 129, “genetic evidence”: what kind of analysis showed this evidence?\nRESPONSE 1.1:\nWe have now clarified the analysis underlying this statement, as follows (line 129, page 3):\n“Admixture dating suggests breeding practices involving continued wild restocking until at\nleast ~350 BCE, with IBE disappearing shortly after.”\nPOINT 1.2:\nLine 242-245: too small to find southern France (SF), northern France and Britain (NF), Denmark\n(FB), Germany (CWC), Poland (POL), and the Carpathian basin\n(Czechia, CZ, and; Hungary, HU) in Fig. 2b and Supplementary Fig. 8. Also, the OUT names in\nFig. 3 are too small to find.\nRESPONSE 1.2:\nWe have now edited the figures as follows, to improve readability:\n-Fig. 2b: we increased the font size, going beyond the minimal size authorized in the instructions\nfor authors; we have also reported labels to highlight those population groups highlighted by the\nreviewer, including those used in AdmixtureGraph modelling;\n-Supplementary Fig. 8: we removed split the autosomal tree in three sections providing zoomed\nin parts of the overall tree; the captions have been updated accordingly;\n-Supplementary Fig. 9: we removed split the mitochondrial tree in four sections providing\nzoomed in parts of the overall tree; the captions have been updated accordingly;\n-Supplementary Fig. 10: we removed split the Y-chromosomal tree in two sections providing\nzoomed in parts of the overall tree; the captions have been updated accordingly;\n-Fig. 3: we increased the font size to match the minimal size authorized in the instructions for\nauthors; we also manually displaced those labels so as to prevent any overlap.\nPOINT 1.3:\nLine 258, which tree is “our phylogenetic analyses”?\nRESPONSE 1.3:\nThe corresponding sentence started with “Combined, our phylogenetic analyses” to indicate that\nall such analyses indicate the same. We have now made this explicit (lines 258-259, page 8):\n“Combined, our autosomal, mitochondrial and Y-chromosomal phylogenetic analyses reveal\nthat the divergent IBE lineage was genetically related to other wild lineages from western\nEurope, especially from France (Late Pleistocene).”\nLine 259, can “genetic affinities” be mentioned from a phylogenetic tree?\nRESPONSE 1.4:\nSee RESPONSE 1.3, as this is about the same sentence. We now write “genetically related to”\ninstead of show “broad genetic affinities” to reflect that the underlying analyses assess\nphylogenetic branching patterns.\nPOINT 1.5:\nFig. 4 model to biogeographic DISCUSSION, but the connection between the model and\nbiogeography is unclear.\nRESPONSE 1.5:\nWe rewrote entirely this to make the connection between the models and biogeography clearer,\nexplicating the geographic locations when necessary and pointing to the specific graph nodes that\nprovided sources in the different expansions discussed. We have also added labels for the\nadmixture nodes (A1 to A6) shown on Fig. 4abc to help readers track the various ancestry\nfractions discussed in the main text. This section now reads (lines 302-339, pages 10-11):\n“All three models depict a complex history of lineage divergence and admixture, in which Late\nPleistocene populations from northern France and UK (LPNFR, NF) emerged as some of the\nearliest lineages investigated to expand and differentiate into Europe (divergence node n8). The\ndescending populations further differentiated into distinct regional lineages, shaped by\ncontributions from other, more divergent groups. For instance, the three models consistently\nsupport the IBE lineage tracking 1.7-6.0% of its ancestry back to a divergent population best\nmodelled as the root and showing some of the largest branch drift estimates (admixture node\nA1; Fig. 4abc). The deep branching of this contribution indicates the presence of an ancient\nhorse lineage that had arrived in Europe much earlier, and survived in the Iberian climatic\nrefugium. Furthermore, the Eneolithic population from Romania (ENEOROM, NR) is modelled\nto primarily derive from a node also forming the URAL (UR) population (n1 on the first two\nmodels, and n6 on the third). This node branches deep in the history of lineage divergence\n(n10), prior to the radiation leading to all other lineages from Europe, indicating yet another\nexpansion into Europe than that forming NF and related descending populations.\nIn all three models, Late Pleistocene populations from southern France (LPSFR, SF)\nderived most of their ancestry from the same expansion that led to the formation of NF\npopulations (admixture node A3), with 77.6-80.9% originating to node n9, immediately\nancestral to n8. However, the remaining fraction of the SF genetic makeup originated from a\npopulation directly related to IBE (n2). As unadmixed members of the IBE lineage were only\nfound in Iberia, thus likely differentiated in Iberia, this suggests that an IBE-related population\nexpanded northward into southern France, where it formed SF by mixing with its local closest\ngeographic neighbours.\nThe three models differ in the population history underlying the formation of the NR\nlineage and those occupying Germany ~2,700 BCE (CWC) and Denmark~3,000 BCE (FBPWC, FB).\nThe first two models indicate that CWC- and FB-related lineages participated to NR\nformation, located further east, likely as populations related to both lineages expanded\neastward towards Romania and mixed with the divergent local lineage (admixture node A6;\nFigs. 4ab). In contrast, the third model supports a reverse scenario, as an NR-related lineage\n(node n4) appears central in the formation of both CWC and FB lineages. This implies that a\npopulation related to the local lineage from Romania (NR) spread westward into Germany and\nDenmark, where it admixed with the descendants of the earlier expansion that established the\nNF lineage (Fig. 4c). While further genome-wide time-series are needed to fully resolve the\ncomplexity of these expansions and admixture events, the three population graphs reveal a\nhistory of significant genetic drift leading to FB and CWC populations. This aligns with early-\nto-mid Holocene bone assemblages, which indicate that horse populations were small and\nsparsely distributed across the continent during this period (Sommer et al., 2011).”\nPOINT 1.6:\nLine 305, I cannot find “19.1–22.4%” from Fig. 4.\nRESPONSE 1.6:\nIn admixture graphs, square nodes represent admixture between two lineages (so-called admixture\nnodes), each contributing in respective proportions summing to 100%. By convention, only one\nof the two proportions is indicated to maximize readability, and the largest proportion is shown,\nwith reference to the contributing population; this helps tracking the main sources contributing to\nthe emergence of each lineage. In the case pointed out by the reviewer, Fig. 4a and Fig. 4b both\nshow the admixture node A3 (leading to LPSFR, SF), which receives 80.9% genetic ancestry\nfrom node n9, and, thus, 100-80.9%=19.1% genetic ancestry from a lineage related to IBE (n2).\nFig. 4c indicates, however, that 77.6% genetic ancestry is received from n9, leaving 100-\n77.6%=22.4% originating from n2. Deviating from this conventional reporting, which is standard\nin the field and common to many publications, will likely be extremely confusing for most\nreaders. We have, thus, decided to rewrite the full section, with explicit reference to specific\nadmixture nodes wherever necessary (see RESPONSE 1.5), and develop the corresponding\ncaptions to help readers understand where the fraction listed come from. The captions now read\n(lines 1029-1034, page 35):\n“Divergence nodes (n), exempt of admixture, are shown with ellipses coloured in dark grey,\nwhile admixture nodes (A) are shown with squares coloured in light grey, with percentages\nreflecting genetic contribution from the main contributing source. The fraction contributed by\nthe second admixture source is the difference between 100% and this percentage. Branches\nassociated with a proportionally larger admixture proportions are highlighted.”\nPOINT 1.7:\nLine 306, I am not sure what the evidence of “IBE differentiated in Iberia” is.\nRESPONSE 1.7:\nWe have developed the rationale for this (lines 321-324, pages 10-11):\n“As unadmixed members of the IBE lineage were only found in Iberia, thus likely differentiated\nin Iberia, this suggests that an IBE-related population expanded northward into southern\nFrance, where it formed SF by mixing with its local closest geographic neighbours.”\nPOINT 1.8:\nLine 309, what is “a deeply branching lineage”?\nRESPONSE 1.8:\nThis sentence has been rephrased, with specific reference to the underlying divergence node. It\nnow reads (lines 314-316, page 10):\n“This node branches deep in the history of lineage divergence (n10), prior to the radiation\nleading to all other lineages from Europe, indicating yet another expansion into Europe than\nthat forming NF and related descending populations.”\nLine 318-319, what does “a significant genetic contribution to the IBE lineage, with 1.7–6.0% of\nits ancestry tracking back to a highly drifted” mean?\nRESPONSE 1.9:\nPopulation graphs not only depict the history of lineage divergence and admixture; they also\nmeasure the genetic drift underlying each branch in the graph. Those estimates are provided on\neach individual branch of Fig. 4, multiplied by a factor 100 (as indicated in the figure captions).\nIn this sentence, we refer to admixture node A1 in all three models, which receive its minor\ngenetic contribution (1.7-6.0%) from a divergent node, best-modelled as the root. The\ncorresponding drift factor is amongst the largest estimated, with values reaching 422.638,\n410.354\nand 207.899 in the three models, respectively. This is worth indicating and may indicate a deep\nhistory of divergence combined or not with a significant history of demographic decline. We have\nrephrased the sentence as follows (lines 307-311, page 10):\n“For instance, the three models consistently support the IBE lineage tracking 1.7-6.0% of its\nancestry back to a divergent population best modelled as the root and showing some of the\nlargest branch drift estimates (admixture node A1; Fig. 4abc). The deep branching of this\ncontribution indicates the presence of an ancient horse lineage that had arrived in Europe much\nearlier, and survived in the Iberian climatic refugium.”\nPOINT 1.10:\nHow can I find the result from Fig. 4?\nRESPONSE 1.10:\nSee RESPONSE 1.6. The 1.7-6.0% ancestry commented can be found from the root to the\nadmixture node in which 98.3% is reported to come from node n7, leaving 100-98.3%=1.7% for\nthe contribution described. The same is true for panel b, excepting that the corresponding branch\nnow points to node n6. In panel c, the corresponding admixture node reports n5:94%, leaving\n100-94% for the contribution described.\nPOINT 1.11:\nLine 346, it is better to explain “(EA, IBE; DOM2*, MNP20_Spa_m1558)” from general readers.\nLine 358: (EA, IBE; DOM2*, X).\nRESPONSE 1.11:\nWe have edited the text to follow the reviewer’s suggestion. It now reads (lines 359-364, page\n12):\n“These analyses align with D-statistics of the form (EA, IBE; DOM2*, MNP20_Spa_m1558),\nexplicitly testing whether MNP20_Spa_m1558 shares more derived mutations with members of\nthe IBE lineage relative to the earliest DOM2 members known in Europe (DOM2*). The\nsignificantly positive D-statistics support an excess of IBE ancestry in MNP20_Spa_m1558\n(Fig. 5a; Supplementary Fig. 16), reaching a fraction of 58.91% (range=57.8-60.0%),\naccording to F4-ratios (Patterson et al., 2012) (Z-score-36,83; Fig. 5c).”\nPOINT 1.12:\nLine 369, “Statistical power was sufficient” is not sufficient. Please show Z scores.\nRESPONSE 1.12:\nWe have expanded the section, providing the Z-scores, as requested. It now reads (lines 384-392,\npage 13):\n4\n“Statistical power was sufficient to estimate that the IBE lineage contributed to 1.22% of the\ngenomic makeup of the admixed specimen from Pech Maho (Z-score=134.72; Fig. 5c), and\n1.65%-1.90% (Z-scores=74.04 and 102.61) in two other specimens from Spain (including\nMR7_Spa_m1131 from the final Bronze Age settlement of El Morredón, which represents the\nearliest Iberian DOM2 horse in our data set (1,213-1,048 cal. BCE)). While Z-scores were also\nlargely significant in all other 12 admixed individuals tested (69.97≤Z-scores≤130.57), the\nconfidence range estimated for the IBE ancestry could intercept zero, or even extend into\nnegative values. This suggests limited sensitivity and that IBE genetic contribution to the\ngenomes of the remaining horses”.\nAnother Z-score was provided for another F4-ratio described earlier in the text. See\nRESPONSE 1.11.\nPOINT 1.13:\nLine 375, why 1.22%? Why did you remove 1.65%-1.90%?\nRESPONSE 1.13:\nWe indicate that F4-ratios succeeded in estimating IBE fractions of genetic ancestries of 1.22%,\n1.65% and 1.90% in three specimens (lines 385-387, page 13). However, a total of 12 other\nspecimens showing evidence of admixture, as indicated by their significantly positive D-statistics\n(Fig. 5a), failed returning sensible F4-ratio estimates (i.e. the confidence range of the IBE-related\nfraction intercepted 0). This methodology, thus, did not have sufficient sensitivity to detect the\ncorresponding fraction of their genome received from IBE ancestors; we can, however, conclude\nthat this fraction is lower than the minimum fraction that was successfully detected amongst the\nthree successful individuals (i.e. 1.22%).\nPOINT 1.14:\nWhat is “density” in Supplementary Fig. 19? F3 value?\nRESPONSE 1.14:\nIn panels a, c, and e, we show the distribution of Haversine distances between pairs of samples,\nwith barplots breaking down the numbers present in each category, and the density showing the\nkernel density estimate smoothed according to Gaussian kernel (using the density() function in\nR). Panels b, d, and f show the same for the time difference between pairs of samples. We have\nrefrained using a density estimate as it seems to have been confusing (without bringing much to\nthe evidence), and now show instead blunt counts as a vertical axis. The edited Supplementary\nFig. 19 is shown below, together with its captions (lines 1684-1690, page 65):\nsamples considered in f3-Outgroup analyses.\nLeft: Distributions of Haversine distances between the pairs of horses considered in the\ncalculation of f3-Outgroup statistics before 1,200 BCE (a), 1,200 BCE–500 CE (c) and after\nthe calculation of f3-Outgroup statistics before 1,200 BCE (b), 1,200 BCE–500 CE (d) and after\nPOINT 1.15:\nFig. 6: Why are close-distance samples not the highest in f3 value? The f3 value from close-\ndistance samples, for example, two samples from the same site, should be the highest.\nRESPONSE 1.15:\nFor each ancient individual, we measure the f3-outgroup genetic affinity to any other ancient\nindividual present in our data set, aiming to reveal potential genetic connection at long geographic\nrange. As indicated in the figure caption: “Pairs involving individuals from the same site or\nseparated by >500 years were disregarded.” This filter is especially implemented (1) to avoid only\ncapturing within site affinities, (2) to avoid biasing the analysis due to uneven sampling across\nsites (i.e. those sites with multiple samples would show local affinities that sites represented with\na single sample could not), and; (3) to only compare those horses living relatively close in time.\nTherefore, pairs of individuals originating from the same site are not displayed. We, however,\nconfirm that individuals from the same site appear more closely related with each other than to\nsamples from other sites. This is visible in Supplementary Fig. 8, in which for example GE1, GE2\n6\na\nc\ne\nstnuoC\nstnuoC\nstnuoC\n6\n4\n2\n6\n4\n2\n0\n0\n0\n3\n2\n1\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n2 0\n2 5 0 0\n5 0 0 0\nv e r s\n0 0\nin e\n1 0\nd\n5\n0is\n0 0\n0 0\nta\n4\n0\n0\n0\n0\n1 5\n7 5 0\n0 0 0\n0\n6\n2\n0\n0\n0\n0\n0\n0\nb\nd\n0\nf\nstnuoC\nstnuoC\nstnuoC\n3\n2\n1\n6\n4\n2\n0\n0\n0\n1\n1\n0\n0\n0\n0\n0\n0\n0\n0\n5\n0\n5\n0\n−\n−\n5\n5\n0\n0\n0\n0\n− 3T\n− 2\n− 2\n0 0\nim\n5 0\n5 0\ne d iffe\n0\n0\n0r\ne n c e ( y\n2\n2\ne\n5\n5\n3a\n0\n0\n0r 0s\n)\n5\n5\n0\n0\n0\n0\nand GE6 samples form a highly-supported monophyletic cluster (the same holds true for CDL10,\nCDL13; CDL15; Z24, Z31, and Z37, and more).\nReviewer #2 (Remarks to the Author):\nPOINT 2.1:\nThis study is a continuation of and dependent on previous work on horse lineages, largely by the\nsame group of researchers. As such there was a tendency to assume knowledge of the previous\nwork throughout. The revision does a good job of addressing this by adding text that elaborates\non previous knowledge and explains their interpretations more fully. This has extended the length\nof the paper, but in a worthwhile way in my opinion.\nRESPONSE 2.1:\nWe thank the reviewer for their constructive suggestion and positive assessment of our work.\nPOINT 2.2:\nThe extent of the data, and the challenges with ancient DNA, make the text, interpretation and\nconclusions complex, raising questions for the reader (some shared by myself and the other\nreviewer). I feel that the authors have made a thorough attempt to address these questions,\nand where appropriate provided useful caveats and cautionary statements.\nRESPONSE 2.2:\nWe thank the reviewer to acknowledge that our conclusions provide sufficient nuance.\nPOINT 2.3:\nOne aspect I feel they could take a bit further is their description and justification of their\nphylogenetic strategies. They note that they had settled on a strategy for their phylogenetic\napproach in earlier studies, but there needs to be some text referring both to the earlier studies\n(now provided) and a brief overview of their reasons for adopting the relatively unusual strategies\nthey use.\nRESPONSE 2.3:\nThe addition of lineages of known phylogenetic affinities (as established here through autosomal\nphylogenetic reconstructions, and in previous population graph modelling; see Librado et al.\nNature 2024)) would significantly complexify the population model, increasing the number of\nnodes, potential admixture edges, and topological space, resulting in significant analytical\nchallenges. The description of a model including more populations would also become extremely\ndifficult to follow, considering the already important difficulty we had summarizing our present\nmodels. We have further expanded the rationale behind our choice, as follows (lines 293-299,\npages 9-10):\n“To further explore the relationship between IBE and other wild horse populations in Europe,\nwe applied AdmixtureBayes (Nielsen et al., 2023) graph modelling to a selection of seven major\ngenetic horse lineages identified in Europe to date (Librado et al., 2021; Librado et al., 2024).\nLineages from the Carpathian basin and Poland, which cluster together in our autosomal\nphylogenetic tree (Fig. 2b, Supplementary Fig. 8), and are known to share genetic affinities with\nthe Late Neolithic CWC (Librado et al., 2024) lineage from Germany, were excluded to limit\nmodel complexity and computational time, while improving M"}]

Decision letter