E Andrew Bennett
email: eabennett@gmail.com

J Weber

W Bendhafer

S Champlot

J Peters

G Schwartz

T Grange

E.-M Geigl

The genetic identity of the earliest human-made hybrid animals, the kungas of Syro-Mesopotamia

Before the introduction of domestic horses in Mesopotamia in the late third millennium BCE, contemporary cuneiform tablets and seals document intentional breeding of highly valued equids called kungas, for use in diplomacy, ceremony and warfare. Their precise zoological classification, however, has never been conclusively determined. Morphometric analysis of equids uncovered in rich Early Bronze Age burials at Umm el-Marra, Syria, placed them beyond the ranges reported for other known equid species. We sequenced the genomes of one of these ~4,500-year-old equids, together with an ~11,000-year-old Syrian wild ass (hemippe) from Göbekli Tepe and two of the last surviving hemippes. We conclude that kungas were F1 hybrids between female domestic donkeys and male hemippes, thus documenting the earliest evidence of hybrid animal breeding.

INTRODUCTION

In the third millennium BCE, urbanized, socially-stratified and literate societies appeared for the first time in Syria and northern Mesopotamia [START_REF] Akkermans | The archaeology of Syria: from complex hunter-gatherers to early urban societies[END_REF][START_REF] Ur | Cycles of Civilization in Northern Mesopotamia, 4400-2000 BC[END_REF]. Part of this "second act" of the urban revolution was the breeding and employment of an equid of high status and prestige designated a "kunga." The precise taxonomical determination of the kunga and its identification in the archaeological record have been uncertain until now. Third millennium BCE cuneiform clay tablets from Syro-Mesopotamia describe several equids, using the generic term ANŠE associated with various logographs. Of these, the so-called "kunga", was represented by the cuneiform signs ANŠE.BARxAN [START_REF] Maekawa | The Donkey and the Persian Onager in Late Third Millennium B.C. Mesopotamia and Syria: A Rethinking[END_REF][START_REF] Maekawa | The Ass and the Onager in Sumer in the Late Third Millennium B[END_REF] (Fig. 1A). Texts from the Diyala region in Mesopotamia and the kingdom of Ebla in the Levant state that the prices for these equids were considerable, costing up to six times the price of a donkey [START_REF] Zarins | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF]. References for these valuable equids are found in multiple clay tablets [START_REF] Maekawa | The Donkey and the Persian Onager in Late Third Millennium B.C. Mesopotamia and Syria: A Rethinking[END_REF][START_REF] Maekawa | The Ass and the Onager in Sumer in the Late Third Millennium B[END_REF] (Fig. 1A) such as those detailing fodder expenses, e.g., barley for the equids of the god Shara and the deified king Shulgi from Umma [START_REF] Heimpel | Towards an understanding of the term SiKKum[END_REF], and dowries for royal marriages [START_REF] Biga | The marriage of Eblaite princess Tagris-Damu with a son of Nagar's king[END_REF]. Large-sized male kungas were used to pull the vehicles of "nobility and gods" [START_REF] Heimpel | Towards an understanding of the term SiKKum[END_REF], and their size and speed made them more desirable than asses for the towing of four-wheeled war wagons [START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF], which predate horse-pulled chariots. Smallersized males and female kungas were used in agriculture, where they were frequently reported pulling ploughs [START_REF] Maekawa | The Ass and the Onager in Sumer in the Late Third Millennium B[END_REF][START_REF] Zarins | The domestication of Equidae in third-millennium BCE Mesopotamia[END_REF]. Kunga foals were seldom born within the urban centers of Sumer and Syria, and Ebla purchased young kungas almost exclusively from what may have been the principal breeding-center at Nagar (modern Tell Brak), in northern Mesopotamia, whose rulers also provided them as gifts to the elites of allied territories [START_REF] Maekawa | The Donkey and the Persian Onager in Late Third Millennium B.C. Mesopotamia and Syria: A Rethinking[END_REF]. Presumed kungas featured prominently on royal seals throughout the region [START_REF] Dolce | Equids as Luxury Gifts at the Centre of Interregional Economic Dynamics in the Archaic Urban Cultures of the Ancient Near East[END_REF], and images of these hybrids likely appear on both the "war" and "peace" panels of the standard of Ur, a Sumerian artefact excavated from the royal cemetery in the ancient city of Ur (in modern-day Iraq). In one of the first depictions (2600 BCE) of a military expedition in human history, warriors stand on four-wheeled war wagons, each drawn by a team of unspecified equids (Fig. 1B). An example of the rein ring featured in this image has been found in a royal grave at Ur (Fig. 1C), decorated with a small statue of a non-caballine equid, either a kunga or hemione. Kunga use and traditions decreased and eventually vanished following the introduction of domestic horses in the region [START_REF] Zarins | The domestication of Equidae in third-millennium BCE Mesopotamia[END_REF][START_REF] Clutton-Brock | Horse Power. A History of the Horse and the Donkey in Human Societies[END_REF]. Early references to horses in cuneiform writing coincides with the Third Dynasty of Ur (late 3 rd millennium BCE), where they are referred to as anše-zi-zi and later anše-kur-ra (equids of the mountain) [START_REF] Heimpel | Towards an understanding of the term SiKKum[END_REF][START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF]. An introduction of domestic horses in Mesopotamia by the end of the third millennium is also supported by paleogenetic data illustrating their late arrival in Anatolia around 2000 BCE, presumably through the Caucasus [START_REF] Guimaraes | Ancient DNA shows domestic horses were introduced in the southern Caucasus and Anatolia during the Bronze Age[END_REF].

While the symbol for kunga (ANŠE.BARxAN) is used to describe a hybrid equid, the unambiguous assignment of this term to a species is difficult and controversially discussed. Some authors even argue that the kunga referred only to wild caught Persian onagers (also known as Iranian onagers; E. hemionus hemionus, a subspecies of the Asiatic wild ass) rather than hybrid animals [START_REF] Maekawa | The Donkey and the Persian Onager in Late Third Millennium B.C. Mesopotamia and Syria: A Rethinking[END_REF][START_REF] Heimpel | Towards an understanding of the term SiKKum[END_REF], although the difficulty in taming modern onagers, which are reportedly less tractable than zebras [START_REF] Mohr | Eine durch Hagenbeck importierte Herde des persischen Onagers, Equus hemionus Onager[END_REF], does not support this interpretation. One of the likely parents of the kunga is the donkey (E. africanus asinus), thought to be present in Sumer from at least the late 4 th millennium BCE [START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF]. The identity of the other parent, however, remains unclear. Another equid attested since the Early Dynastic period I (ca. 2800 BCE) is the anše-edin-na, literally translated as "equid of the desert". This animal was hunted for its meat and hide, but never used as a draught animal. The anše-edin-na is broadly considered to be a type of onager [START_REF] Maekawa | The Ass and the Onager in Sumer in the Late Third Millennium B[END_REF][START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF], although it is impossible to say whether it refers to the Persian onager (E. hemionus onager) or to the Syrian wild ass, or hemippe (E. hemionus hemippus) [START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF], sometimes also named "Syrian onager". Described as a light, swift animal [START_REF] Antonius | Beobachtungen an Einhufern in Schönbrunn 1: Der Syrische Halbesel (Equus hemionus hemippus J. Geoffr.)[END_REF], the hemippe was the smallest of all modern equids until the subspecies went extinct early in the 20 th century [START_REF] Clutton-Brock | Horse Power. A History of the Horse and the Donkey in Human Societies[END_REF]. It has been argued that in addition to its untamable, aggressive nature [START_REF] Antonius | Beobachtungen an Einhufern in Schönbrunn 1: Der Syrische Halbesel (Equus hemionus hemippus J. Geoffr.)[END_REF], its diminutive size made it an unlikely candidate for use in breeding kungas (3 and references therein). Some authors considered the non-donkey parent to be a horse (discussed in 4, 9).

In the elite burial complex of Tell Umm el-Marra (2600-2200 BCE), possibly belonging to the ancient city of Tuba, 55 km east of Aleppo in modern-day northern Syria (Fig. 2), men and women were interred with ceramics, bronze and silver vessels, bronze weapons and tools, and personal ornaments made of bronze, silver, gold and lapis lazuli [START_REF] Schwartz | A third millennium BC elite mortuary complex at Umm El-Marra, Syria: 2002 and 2004 excavations[END_REF]. Within this royal burial complex, complete skeletons of 25 male equids and bones from 6 additional animals were buried separately from humans, either in a sequence of pits or in their own mud-brick structures [START_REF] Schwartz | A third millennium BC elite mortuary complex at Umm El-Marra, Syria: 2002 and 2004 excavations[END_REF], akin to the 3000 BCE donkey burials at Abydos, Egypt [START_REF] Rossel | Domestication of the donkey: Timing, processes, and indicators[END_REF]. While some animals were interred after natural deaths, more than half appear to have been deliberately killed for burial in the complex. Morphometric values obtained from these bones indicate that these animals constitute a population outside of the typical ranges of horses, asses, and onagers, and it has been proposed that these skeletons represent hybrids, presumably kungas (17 and Supplementary Materials). In absolute size, the skeletons are closer to hemiones, but are more robust; commonly-used slenderness indices suggest greater affinities with asses than with hemiones. The leg characteristics of hemiones, responsible for a speed exceeding that of horses, is retained in these animals, suggesting that they were also fast [START_REF][END_REF]. Interestingly, discrepancies in wear between the incisors and cheek teeth of some of the equids indicate that the animals were foddered and not commonly grazed [START_REF] Weber | Archaeozoology of the Near East VIII[END_REF], features that would have been expected on the skeletons of the equids depicted on the standard of Ur, whose lip-or nose-rings would have made grazing difficult (Fig 1B). These animals would have been stronger and faster than donkeys, and must have been more tamable than hemiones [START_REF] Oates | Prehistoric Steppe Adaptation and the Horse[END_REF].

Taxonomic classification of equids uncovered in tombs across Mesopotamia (Ur, Kish, Lagashnow al-Hiba, Abu Salabikh, Tell Madhhur) is often controversial (for discussion see e.g., [START_REF] Zarins | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF][START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF][START_REF] Clutton-Brock | The animal remains from Abu Salabikh, a preliminary report[END_REF][START_REF] Clutton-Brock | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF]. The degree of variation within ancient populations is not fully known, and the degree of variation between individuals within a population -especially of domesticated animals -is large, making it difficult to differentiate between E. africanus and E. hemionus using solely bone morphological and metrical characteristics (for discussion see 22).

In order to clarify whether the burials of Tell Umm el-Marra indeed contained the remains of the politically and symbolically important hybrids referred to in numerous cuneiform tablets as kunga, and in order to determine the taxonomic status of those animals, we investigated the genomes in samples from the skeletons of the equid installations at Umm el-Marra, an equid sample from the Early Neolithic site of Göbekli Tepe (Turkey), as well as the last survivors of the Syrian wild ass conserved in the Natural History Museum of Vienna.

RESULTS AND DISCUSSION

Analysis of the maternal and paternal lineages of the Umm el-Marra equids

An initial PCR screening of equid samples from Umm el-Marra showed that DNA was extremely poorly preserved in these bones owing to the hot climate in Syria, detrimental to long-term DNA preservation, and the poor condition of the bones (phalanges, sternum) available for study (Fig. S1).

Therefore, we combined shotgun nuclear DNA sequencing with highly sensitive PCR, targeting taxonomically informative regions of both uniparental markers: mitochondrial DNA and the Ychromosome. To better pinpoint the genetic identity of the parental species, we increased the available Y-chromosome data by sequencing regions from additional populations of both modern and 19th and 20th century museum samples of hemiones and donkeys, for which the mitochondrial sequences were previously generated [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF]. Short, overlapping PCR products suited to the degraded DNA of the samples were designed to amplify a highly diagnostic mitochondrial control region fragment (324 base pairs long), including the site of a well-characterized 28 base pair-deletion exclusive to hemiones [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF], and three separate regions of the Y chromosomal DNA (in total 168 base pairs long) encompassing four single nucleotide polymorphisms (SNPs), which we show to be diagnostic between E. ferus (caballus and przewalskii), E. africanus, and E. hemionus.

The full targeted mitochondrial sequence was successfully amplified from two out of the six individual equids tested from Umm el-Marra. At every position divergent between E. ferus, E. africanus, and E. hemionus, both of these sequences contained the E. africanus specific bases and lacked the 28 base pair-deletion specific to E. hemionus [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF]. The maternal lineage of these equids thus unambiguously belongs to E. africanus as visible in a median joining network (24) (Fig. 3A).

All three Y-chromosome fragments were successfully amplified from these same two Umm el-Marra individuals. Within these three regions, four diagnostic positions differentiate E. africanus

(T/G/T/A) from E. hemionus (C/A/G/G), two of which also differentiate E. ferus (C/G/G/A) from either E.
hemionus or E. africanus. At each diagnostic position, the equids from Umm el-Marra were found to have the E. hemionus-specific base, and no diagnostic position of any product contained the E. africanus-specific base (table S1). The hemione-specific Y-SNPs were also confirmed previously in diverse hemiones from archaeological samples from the Caucasus, museum specimens from Tibet and Syria and present-day specimens from the Gobi in Mongolia (23) (Fig. 3B and table S1).

In addition to the E. hemionus diagnostic positions, both Umm el-Marra sequences contained two additional Y-chromosome SNPs observed only in the two hemippes from the 19th and 20th centuries analyzed herein (Fig. 3B), one of them being the last known member of the subspecies. This animal had been caught in the deserts north of Aleppo in 1911 and had been kept in the Schönbrunn Zoo in Vienna until its death in 1929 (see Fig. S2 for images of two of the hemippes used in this study).

Thus, the Umm-el-Marra equids harbor the maternal lineage of the domestic donkey and the paternal lineage of the Syrian wild ass, suggesting they could be F1 hybrids, since interspecific equid hybrids are generally sterile, or poorly fertile.

Analysis of the nuclear genomes of Umm el-Marra and Göbekli Tepe equids and the last

Syrian wild asses

To further establish the hybrid identity of these equids, we sequenced a subset of the nuclear genome of the best preserved Umm-El-Marra equid bone (table S5). In addition, we established the genome sequence of the extinct Syrian hemippe by sequencing a ca. 11,000-year-old wild ass from the early Neolithic site of Göbekli Tepe, present-day Turkey, representing the first temple [START_REF] Schmidt | Sie bauten die ersten Tempel[END_REF], and two 19 th century specimens from the Schönbrunn Zoo (table S5). These four newly generated genomes were compared to six modern horse [START_REF] Orlando | Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse[END_REF], six domestic donkey [START_REF] Wang | Donkey genomes provide new insights into domestication and selection for coat color[END_REF], three Mongolian khulans [START_REF] Wang | Donkey genomes provide new insights into domestication and selection for coat color[END_REF][START_REF] Huang | Donkey genome and insight into the imprinting of fast karyotype evolution[END_REF] (an E. hemionus subspecies from the Gobi), two kiang genomes [START_REF] Wang | Donkey genomes provide new insights into domestication and selection for coat color[END_REF][START_REF] Jonsson | Speciation with gene flow in equids despite extensive chromosomal plasticity[END_REF]

(E.
hemionus kiang or E. kiang), and one Persian onager genome (29) (an E. hemionus subspecies from Iran), (table S4). A set of 15.5 million SNPs residing outside of repeated sequences and being variable in the modern genome equid panel was used for calling the ancient genomes (see Supplementary Methods and table S5). Even though the best Umm-el-Marra extract contained only 0.18% endogenous DNA, we could obtain 45.6K SNPs, 40.4K of which were shared with either hemippe.

Of these, 15.2K SNPs (37%) were shared with both of the two best-covered hemippe genomes. First, we performed principal component analysis (PCA) to compare (a) the hemione and donkey genomes (Figs. 4A, S4A), or (b) hemione, donkey and horse genomes (Figs. S4B-E). Identical results were obtained whether we used the 15.2K SNPs shared between the Umm-el-Marra equid (UMM9), the Göbekli Tepe and the 1864 hemippe, or whether we used the 4.1 million SNPs shared between the Göbekli Tepe and the 1864 hemippe genome and projected the UMM9 equid onto the PCA (compare Figs.4A and S4A, S4B and S4D, S4C and S4E), showing that the 15.2K SNPs obtained allowed robust characterization of the status of the Umm-el-Marra equid. When only donkey and hemione genomes are used, PC1 separates donkeys from hemiones and PC2 separates the hemiones (Figs. 4A andS4A). The most differentiated in PC2 are the Persian onager and the kiang zoo specimens. The Mongolian khulan and the kiang from neighboring regions in China are very closely related and overlapping in the PCA, which is in accordance with the shared mitochondrial lineages we reported previously that led us to question the specific taxonomic status of the kiang as a separate species [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF] (see also the phylogenetic trees of the mitogenomes and the genomes in Figs. S6 andS7). The two modern hemippes as well as the Göbekli Tepe sample overlap as well (the lower coverage 1892 hemippe was projected) and are located at an intermediate position in PC2. The Göbekli Tepe sample, a mare, is clearly an ancient hemippe as observed from the phylogenetic trees constructed from both mitochondrial and nuclear genomes (Figs. S6-S7). The UMM9 equid falls exactly halfway between the donkeys and the hemippes. When the PCA also includes the horses, PC1 differentiates the horse from the non-caballine equids, PC2 separates the donkeys from the hemiones, and PC3 separates the hemiones in a similar way as PC2 does when horses are not included (Figs. S4B-E). In all analyses, the results for the UMM9 equid illustrate an intermediate position between the donkeys and hemippes.

The PCA analyses thus indicate that the UMM9 genome is a 50% mixture of donkey and hemippe.

We further explored this outcome through ADMIXTURE analysis (30) (Fig. 4B). A 4-population model separates horses, donkeys, Mongolian khulans and kiangs from onagers and the two bestcovered hemippes. The UMM9 equid is modelled as an admixture of equal proportion between donkey and hemippe/onager (Fig. 4B). Likewise, when considering the 4738 UMM9 equid SNPs for which all six donkeys differ from all hemippes that have the position covered, the UMM9 equid harbors SNPs corresponding to roughly half of those specific to each putative parent (Fig. 4C).

Finally, a bifurcating tree with gene flow analysis was performed using treemix (31) (Fig. 4D).

For the reference equids, genomic tree topology is similar in both the full mitogenome and genome tree topology obtained with different methods (Figs. S6-S7). The Persian onager and the hemippes are closely related, and the Mongolian khulans and the Tibetan kiangs are even more closely related.

In this respect, the genetic distances between the various hemiones correspond to the geographic distances between their native range (onager: Iran; hemippe: Syria; Mongolian khulan: Mongolia; kiang: Tibet). The UMM9 equid is represented on the tree as related to the donkey, but the residuals between the hemippe and the UMM9 equid are high and a gene flow event from the hemippe to the UMM9 equid best describes the phylogeny (see also Fig. S5). These results demonstrate the sufficiency of the Umm el-Marra SNPs to determine the phylogenetic relationships between the equids. The tree also placed the Umm el-Marra sample halfway between the asses and the hemippes (Fig. 4), which meets expectations when dealing with F1 hybrids of these species. Evidence from Ychromosome analysis indicates that the Syrian hemippe rather than the Persian onager was used to father the Umm el-Marra equids, whereas the mitochondrial DNA reveals that a donkey contributed the maternal genome. That the Umm el-Marra equids were in fact F1 hybrids and not back-crossed hybrids is also supported by the relative hybrid sterility between donkeys and horses, as well as experiments in the 1940s crossing female donkeys with male hemiones, the Turkmenian kulans (E. hemionus kulan), which produced sterile offspring [START_REF] Ščekin | Rezultaty gibridizacii lošadei i oslov s kulanami[END_REF].

Expectedly, the 19 th -20 th century hemippes, representing some of the last survivors of the subspecies are genetically similar, whereas the ~11,000-year-old Göbekli Tepe hemippe is more divergent (Fig. S7). We also noted that the divergence between the three sequenced hemippes is much larger than that observed between the six domestic donkeys (Fig. S7). The higher diversity between the sequenced hemippe genomes versus between the donkey genomes suggests that the donkey mother of the UMM9 equid is more closely related to present-day donkeys than the Syrian hemippe father of the UMM9 equid is to the other sequenced hemippes. This difference, albeit small, may account for the slightly higher affinity of the UMM9 equid to present-day domestic donkeys that is visible on Figs. 5 B &D. Both the genomic phylogenetic tree (Fig. S7) and the PCA analyses (Figs.

4B & S4) indicate that differentiation between the donkey genomes is low, far less than between the various present-day hemiones, presumably because donkeys went through a major bottleneck, possibly upon domestication and translocation to southwest Asia outside the range of the ancestors of donkeys. The observation that the ~4500-year-old UMM9 equid appears more closely related to present-day donkeys than to the last hemippe that disappeared a century ago suggests that the bottleneck of the donkey population had already taken place by the third millennium BCE.

It has been noted that the Syrian wild ass (hemippe), whose range once extended across the Levant, was the smallest form of modern equids [START_REF][END_REF]. Both historical specimens analyzed in this study stood ca. 100 cm at the shoulder (14) (Fig. S2). In contrast, the hybrids of Umm el-Marra were estimated to average 130 cm at the shoulder [START_REF] Weber | Archaeozoology of the Near East VIII[END_REF]. Regarding this difference in size, previous work had recovered mitochondrial haplotypes from larger-sized Bronze Age equids recovered from Tell Munbaqa, situated in northern Syria east of Umm el-Marra, as well as from three historical hemippe samples dating from the mid-19th to early 20th century. Both the larger ancient and smaller more recent animals were shown to cluster together in a single, separate mitochondrial clade [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF]. The Göbekli Tepe wild asses were on average even slightly larger than those of roughly contemporaneous Tell Mureybet (10 th -9 th millennia BCE) and 3 rd -2 nd millennia BCE Tell Munbaqa, two sites located in the direct vicinity of Umm el-Marra ( 22), (Fig. 2). It was concluded, therefore, that the small Syrian wild ass was likely to have been a dwarfed descendant of a genetically continuous population of larger, more robust animals populating Syria in the third millennium BCE and earlier [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF]. The genomic analyses of both ancient and historical hemippes in the present study support this earlier finding. Indeed, no dwarf form has ever been reported from Late Pleistocene and Holocene sites in Mesopotamia [START_REF] Bökönyi | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF] or Anatolia [START_REF] Meadow | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF]. Nearly 2,000 years after the equid burials of Umm el-Marra, sixth century BCE palace reliefs featuring hunted hemiones from Nineveh (in modern-day northern Iraq) show already relatively small animals (Figs. 1D &S8).

CONCLUSION

The genomic results from the rare equid burials at the elite mortuary complex of Umm el-Marra confirm earlier hypotheses based on morphological data that these animals are indeed hybrids [START_REF] Postgate | Equids in the Ancient World (Beihefte zum Tübinger Atlas des Vorderen Orients, Reihe A[END_REF][START_REF] Zarins | The domestication of Equidae in third-millennium BCE Mesopotamia[END_REF][START_REF] Champlot | An efficient multistrategy DNA decontamination procedure of PCR reagents for hypersensitive PCR applications[END_REF] and, given their interment in high-status tombs, are most likely identical with the valuable kungas frequently mentioned in cuneiform texts and depicted in images and royal seals throughout Mesopotamia. This study now offers a firm zoological classification of the historical kunga as an F1 cross between a female donkey and a male Syrian wild ass, or hemippe, putting to rest past speculations regarding the taxonomic identification of the BAR x AN. We further show that the 3 rd millennium BCE ancestors of the hemippe were likely larger than those first described by European travelers visiting Syria in the 19 th century. Our study also presents the earliest known case in human history of inter-species hybridization practiced by Early Bronze Age breeders at sites such as Nagar (Tell Brak) (Fig. 2), to generate animals famous for their power, both physical and symbolic, in ancient warfare and diplomacy. This result also deepens our insight into the economic and political relationships between contemporary royal households of Greater Mesopotamia, and the dynamics by which these social elites fostered distant alliances. It also increases our understanding of the ways in which the earliest stratified urban societies of the Middle East developed and maintained their positions of authority. In this respect, genomic characterization of additional equids from comparable contexts, particularly from Nagar, may help clarify the scale of hybrid breeding in 3 rd millennium BCE Mesopotamian societies prior to the introduction of domestic horses.

Materials & Methods

Sample description

Bone remains from equid skeletons dated between ca. 2550 and 2300 BCE and excavated in 2006 at Tell Umm el-Marra, a Bronze Age elite cemetery in northern Syria [START_REF] Schwartz | A third millennium BC elite mortuary complex at Umm El-Marra, Syria: 2002 and 2004 excavations[END_REF][START_REF] Weber | Archaeozoology of the Near East VIII[END_REF][START_REF] Schwartz | From Urban Origins to Imperial Integration in Western Syria: Umm el-Marra 2006[END_REF] were sampled for ancient DNA analysis. Further descriptions of the samples are given in the Supplementary Materials and photos of the samples from the two individuals from which sufficient DNA was recovered appear in Fig. S1. A petrous bone excavated from a layer dated between 9500-8300 BCE from the site of Göbekli Tepe in southeast Turkey [START_REF] Schmidt | Sie bauten die ersten Tempel[END_REF] and Supplementary Material) was also analyzed in the present study. Furthermore, we analyzed two samples of the extinct Equus hemionus hemippus originating from the desert of Aleppo in Syria and kept in the zoo of Schönbrunn, Vienna, Austria, a tooth from the NMW 6048/ST345 specimen and a hair and skin sample from the NMW1308/B4690 specimen, corresponding to animals who died in the Schönbrunn zoo in 1864 and 1892, respectively (36) (Fig. S2). Finally, hair of a male Somalian ass (E.africanus somaliensis) from the "Réserve Africaine de Sigean" (Sigean, France) was provided for the analysis of the Y chromosome by Dr. Elodie Trunet (Sample "As.Somalie").

aDNA Extraction, Amplification, and Sequencing

Hair samples were added to 1.5 mL hair digestion buffer (100 mM Tris-HCl pH 8.0, 100 mM NaCl2, 40 mM DTT, 3 mM CaCl2, 2% N-lauryl sarcosyl, 250 µg/mL proteinase K) and incubated 4-24 hours at 50°C, shaken 300 RPM. Solutions were then pelleted and the supernatant was purified using a QIAquick Gel extraction kit (Qiagen, Hilden, Germany) according to instructions. DNA from archaeological bone samples was extracted, purified and prepared for either qPCR or sequencing in the ancient DNA lab described previously [START_REF] Champlot | An efficient multistrategy DNA decontamination procedure of PCR reagents for hypersensitive PCR applications[END_REF][START_REF] Bennett | Library construction for ancient genomics: single strand or double strand?[END_REF][START_REF] Pruvost | Freshly excavated fossil bones are best for amplification of ancient DNA[END_REF]. Bone cleaning and treatment protocols were as described previously [START_REF] Guimaraes | Ancient DNA shows domestic horses were introduced in the southern Caucasus and Anatolia during the Bronze Age[END_REF][START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF]. Briefly, after removal of the surface with a razor's blade or surface-cleaning with bleach the bones were either sawed using a flame-sterilized diamonddisc of a Dremel Fortiflex (Dremel Europe, The Netherlands) and ground to fine powder in a 6775

Freezer/Mill®Spex SamplePrep in liquid nitrogen or drilled at low speed with a flame-sterilized bit.

The dense pyramidal part of the petrous bone GT64 was isolated using a flame-sterilized diamonddisc of a Dremel and then ground to fine powder in a 6775 Freezer/Mill®Spex SamplePrep in liquid nitrogen. Half of the GT64 powder was treated with diluted hypochlorite (1/20) and both halves washed with phosphate buffer according to Korlević et al., 2015 (40). DNA extraction was performed by incubating the bone powder at 37°C for 48 to 90 hours in either 1 to 10 mL extraction buffer A (0.5 M EDTA, 0.25 M PO4 3-pH 8.0, 0.14 M ß-mercaptoethanol) or in twice 1 mL extraction buffer B (0.5 M EDTA, 0.05% Tween20, 250 μg/mL Proteinase K, 0.14 M ß-mercaptoethanol) that were pooled prior to purification. Samples were purified using silica membrane spin-columns (QIAquick Gel Extraction kit) with a vacuum manifold (Qiagen) and 25 mL extenders (Qiagen) as described [START_REF] Bennett | Library construction for ancient genomics: single strand or double strand?[END_REF] [START_REF] Gorgé | Analysis of Ancient DNA in Microbial Ecology[END_REF], as well as with either the 5M guanidine HCl, 40% isopropanol (5M40) buffer as described in Dabney et al [START_REF] Dabney | Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments[END_REF], or the 2M guanidine HCl, 70% isopropanol (2M70) buffer as described in Glocke and Meyer [START_REF] Glocke | Extending the spectrum of DNA sequences retrieved from ancient bones and teeth[END_REF]. The elution was performed twice in 25 μl 10 mM Tris-HCl pH 8.0, 0.05% Tween-20 (EBT) made from gamma-irradiated water (8 kGy).

Purified DNA was amplified by qPCR, the extract making up 5-20% total volume (10-20µl per reaction). Inhibition characteristics were determined for failed samples indicating possible inhibition, and once optimal dilutions were determined, qPCR was attempted again. To protect against crosscontamination, the UQPCR method was used [START_REF] Champlot | An efficient multistrategy DNA decontamination procedure of PCR reagents for hypersensitive PCR applications[END_REF][START_REF] Pruvost | Freshly excavated fossil bones are best for amplification of ancient DNA[END_REF][START_REF] Pruvost | Minimizing DNA contamination by using UNG-coupled quantitative realtime PCR on degraded DNA samples: application to ancient DNA studies[END_REF], in which uridine was substituted for thymidine in all PCRs, and incubation with uracil N-glycosylase (UNG, extracted from G. morhua; Biotec Marine Biochemicals, Norway) was performed prior to each reaction. Mock extracts were included with each extraction and amplified to control for contamination. qPCR reactions varied slightly depending on the sample, but a typical reaction included 1.77µL of LC FastStart DNA MasterPLUS mix1b and 0.23uL of either FastStart DNA MasterPLUS mix1b, mix1a, or FastStart Taq (Roche Applied Science, Mannheim, Germany), a final concentration of 1µM of each primer and 1U per reaction of UNG in 10µL total volume. Primers were obtained from Sigma-Aldrich (St. Louis, USA). Mitochondrial primers were designed to amplify 357 base pairs of the hyper variable region of E. africanus and E. hemionus mitochondria using short, overlapping fragments (table S3). Ychromosome primers were designed to amplify three short sections of Y-chromosome DNA containing the target SNPs (table S1-S3). Several modifications of these primers were designed to increase sensitivity of qPCRs by minimizing the likelihood of primer-dimers and artifacts, and increasing primer efficiency. A list of primers used and product sizes is given in table S3. qPCR was performed using Lightcycler 1.5 or Lightcycler 2 (Roche Applied Science, Mannheim, Germany). qPCR programs varied depending on primer requirements and product length, but a typical program involved UNG incubation at 37°C for 15 minutes, followed by polymerase activation at 95°C for 5 minutes, then 2-step cycles of denaturation at 95°C for 10 seconds, then primer annealing and extension at 62°C for 40 seconds, and finally a temperature increase of 0.1°C/1second from 62°C to 95°C with continuous fluorescence measurement to generate melt-curves of the products. Products were purified by a QIAquick PCR purification kit (Qiagen, Hilden, Germany) and both strands were sequenced by capillary electrophoresis at Eurofins/MWG Operon (Ebersberg, Germany) using an ABI 3730xl DNA Analyzer (Life Technologies). Samples that yielded sequence results for the Ychromosome are shown in table S1. An average of one non-template control (NTC) was run for every 6.6 samples (including mocks). No DNA was amplified in either NTCs or mocks, demonstrating that no detectable equid DNA was introduced during sample preparation or was present in reagents.

Samples from the two individuals from Umm-el-Marra with the best preserved DNA identified via qPCR (UMM4 and UMM9) were selected for shotgun sequencing, and 24 double-stranded libraries using dual-barcodes were prepared from DNA purified from three or four different areas of each bone using the protocol described in Massilani, et al., 2016 (45). Seven of these libraries were treated with UDG to reduce the presence of cytosine deamination damage in the resulting sequences. Double-stranded libraries of the two E. h. hemippus (hemippe) museum specimens were constructed using the NxSeq ampFREE Low DNA Library kit (Lucigen, Middleton, WI, USA) following the protocol and the modifications described in Bennett, et al., 2019 (46). Barcodes were added during an amplification reaction using dual-barcoded single-stranded library adapters [START_REF] Gansauge | Single-stranded DNA library preparation for the sequencing of ancient or damaged DNA[END_REF] as primers, rather than those in the kit, where 20 μl of eluted library was added to 25 μl OneTaq 2x Master Mix (Roche) and 0.6 μM of each adapter for 50 μl total volume, and amplified with the following protocol: 5 minutes 95°C, 30 cycles of 30 seconds 95°C, 30 seconds 60°C, 45 seconds 68°C, followed by a 5 minute cycle at 68°C. A library for a hair and skin sample belonging to specimen NMW5493/B 3625 was constructed using the NEBNext® Ultra™ II DNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA).

Libraries from the Göbekli Tepe petrous bone GT64 extracts were constructed using the NEBNext® Ultra™ II DNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) after a pretreatment with USER enzyme mix (NEB, Ipswich, MA, USA).

Dual-barcoded libraries were then purified and size-selected using NucleoMag beads (Macherey-Nagel) for two rounds of purification following the supplied protocol at a ratio of 1.3x beads per reaction volume and eluted in 30 μl EBT.

All libraries were quantified with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific), Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and by qPCR. Screening by shotgun sequencing of Umm-el-Marra samples and of the two hemippe samples was performed on an Illumina MiSeq system using a v3 reagent kit for 2x75 cycles. The libraries constructed from sample SP345, which came from a molar belonging to specimen NMW6048/ST345 (1864 hemippe), a hair and skin sample belonging to specimen NMW1308/B4690 (1892 hemippe), and the two libraries from the Umm-el-Marra 9 (UMM9) sample purified using the 5M40 and 2M70 buffers were selected for deep genomic sequencing.

Sequencing of the 1864 hemippe genome was performed on an Illumina NextSeq using a NextSeq 500/550 High Output Kit v2 (2x75 cycles). The custom sequencing primer CL72 (47) was substituted for the read 1 primer sequencing steps, which is compatible with the single-stranded adapters used for these samples. Sequencing of the 1892 hemippe and of two of the UMM9 libraries was performed first on an Illumina MiSeq using a v3 reagent kit for 2x75 cycles and then on an Illumina Novaseq damaged reads typical of ancient DNA libraries, in particular when mapping reads from a noncaballine equid to the horse reference genome because these genomes are expected to differ markedly in repeat location and sequence variability. This filtration reduced the variant list to 22 million. We then called variants from this list on the modern equid genomes using bcftools (v1.9; (55)) mpileup -B -q30 -Q30 and bcftools call -m. The vcf file was imported in plink (v1.9; (56)), and filtered to include only SNPs, removing invariant and multiallelic positions. The final curated list contains 15.5 million SNPs.

Shotgun reads for the Umm-el-Marra (UMM9) and the Göbekli Tepe (GT64) samples were merged with leeHom (57) using the ancientdna option, while the historical hemippe reads were trimmed with cutadapt (v1.18) [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Fragments smaller than 28 bp were discarded, and the remaining reads were aligned to the E. caballus reference genome (eqCab2.0) using the BWA (v0.7.17) (51) aln program with parameters "-n 0.01 -l 0" followed by samse (UMM9 and GT64) or sampe (hemippe). PCR duplicates were removed using picard MarkDuplicates (v2.20.0)(52), and reads aligning to the reference genome below a mapping quality score of 20 and a length of 28 bp were removed. To reduce the increase in spurious alignments from shorter reads described in [START_REF] De Filippo | Quantifying and reducing spurious alignments for the analysis of ultra-short ancient DNA sequences[END_REF], mapped reads less than 35 bp containing indels were also removed using an awk script. The ancient nature of the UMM9 and GT64 sequences was confirmed by analyzing the damage profile using mapDamage2 (59) of libraries generated from extracts not treated with USER-enzyme (Fig. S3). To remove the C->T mutations at the end of the molecules that escaped the USER treatment (Fig. S3), the base quality was rescaled at the last 2 bases using mapDamage2 [START_REF] Jónsson | mapDamage2.0: fast approximate Bayesian estimates of ancient DNA damage parameters[END_REF]. Since all Umm el-Marra samples had very low levels of endogenous equid DNA (zero to 0.18% of reads), they were additionally aligned to the human and bovine genome reference sequence (GRCh37 and ARS_UCD1.2 respectively). Only libraries that had at least 5-fold more reads mapping to the horse than to the cow or human genome when a seed length of 18 was used during bwa aln mapping were kept. Summaries of the sequencing results are given in table S5. Hemippe reads were additionally aligned to the kiang mitochondrial genome (NC_016061.1) (60) using bwa aln and bwa mem. The resulting E. h. hemippus correction was disabled. The tree obtained with one migration/admixture event is represented in Fig. 4D and the residuals are plotted in Fig. S5. From these 40.4K SNPs, we identified those where all donkeys are identical and all hemippes that have the corresponding positions covered are identical and distinct from the donkeys (4738 SNPs) and counted the SNPs where the UMM9 SNPs are identical to either the donkey-or the hemippe-specific SNPs. These counts are represented in Fig.

4C.

The genome phylogeny shown in Fig. S7 was obtained using the 738.5K SNPs shared between all three hemippes and equids, after calculation of the pairwise distance matrix between all equids using plink [START_REF] Chang | Second-generation PLINK: rising to the challenge of larger and richer datasets[END_REF] and construction of the phylogenetic tree using fastme with nni optimization [START_REF] Lefort | FastME 2.0: A Comprehensive, Accurate, and Fast Distance-Based Phylogeny Inference Program[END_REF].

  ©K. Maekawa). (B) Detail from the Standard of Ur shows an equid team pulling a four-wheeled wagon in battle (photo: ©The British Museum Images). (C) An image of a rein-ring with decorative equid from a royal grave at Ur, contemporary and similar to those visible in the Standard of Ur. (D) Nineveh panel: "hunting wild asses" (645-635 BCE) (British Museum, London). Fig. S8 shows additional panels attesting that the equids depicted are non-caballine. C and D, British Museum, London, photos: E. Andrew Bennett.

Figure 2 .

 2 Figure 2. Map of 3 rd millennium BCE Syro-Mesopotamia showing the major historical and

Figure 3 .

 3 Figure 3. Median joining network of equid sequences. (A) 324 bp of the hypervariable region of

Acknowledgments

We are very grateful to Kazuya Maekawa for discussion of the cuneiform terms and for providing a drawing of the cuneiform tablet BM 23836, conserved in the British Museum, London. We are also grateful to Frank Zachos and Alexander Bibl from the "Naturhistorisches Museum Wien" for providing the hemippe specimens. We thank Gerhard Heindl from the Geschichtsforschung & Dokumentation department of the Schönbrunner Tiergarten GmbH, Vienna, Austria, for help with the search for photos of the last hemippe. We thank Elodie Trunet Réserve Africaine de Sigean, France, for providing the sample of the Somalian ass. We thank Tanja Kovaleva for translations of Russian articles, Caitlin Martin for critical reading of the manuscript and Olivier Gorgé for assistance with some of the sequencing.

Funding

The paleogenomic facility of the Institut Jacques Monod obtained support from the University Paris Diderot within the program "Actions de recherches structurantes". The sequencing facility of the Institut Jacques Monod, Paris, is supported by grants from the University Paris Diderot, the Fondation pour la Recherche Médicale (DGE20111123014), and the Région Ile-de-France (11015901).

Moreover, we acknowledge support from the French national research center CNRS. We are grateful to the Directorate-General of Antiquities and Museums, Syria, for its support of the Umm el-Marra project. The excavations at Umm el-Marra were funded by the National Science Foundation (Grants BCS-0137513 and BCS-0545610), the National Geographic Society, the Metropolitan Museum of Art, the Arthur and Isadora Dellheim Foundation, and the Johns Hopkins University. Faunal research at Göbekli Tepe was funded by the Deutsche Forschungsgemeinschaft (DFG) under Grants PE 424/10-1-4 to JP.

Data availability

Sequence data generated for this study are available from EBI European Nucleotide Archive PRJEB47929. Syrian wild ass (hemippe) mitochondrial sequences are available on Genbank MN990427, OK393913 and OK393914.

6000 using a S2 flow cell for 2x50 cycles. Sequencing of the GT64 libraries was performed on a Novaseq 6000 using a S4 flow cell for 2x75 cycles.

Paleogenetic data analyses

Sequences from PCRs were manually curated, assembled and aligned using the Geneious software suite [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF]. Median-joining network analysis [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF] was performed on mitochondrial sequences covering 357 bp of the hyper variable region generated in this study by PCR combined with those previously reported [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF] (accession numbers given in table S2), and Y-chromosome sequences generated in this study combined with those publicly available (samples and sources shown in table S1). Maximum likelihood (ML) analyses of the complete mitochondria of the two hemippes combined with donkey and hemione complete mitochondria after deletion of the 11-bp-tandem repeat in the HVR were computed using RAxML [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF] with a GTR nucleotide substitution model, a gamma-distributed rate of variation among sites with four rate categories, and invariant sites (i.e., GTR-GAMMA-I) (Fig. S6). We used 100 bootstraps to estimate node robustness.

Genomic analyses

Fastq reads from six modern horse genomes [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF], six domestic donkey genomes [START_REF] Schmidt | Sie bauten die ersten Tempel[END_REF], three Mongolian khulan genomes [START_REF] Schmidt | Sie bauten die ersten Tempel[END_REF][START_REF] Orlando | Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse[END_REF] (an E. hemionus subspecies from the Gobi), two kiang genomes [START_REF] Schmidt | Sie bauten die ersten Tempel[END_REF][START_REF] Wang | Donkey genomes provide new insights into domestication and selection for coat color[END_REF]) (E. hemionus kiang or E. kiang), and one Persian onager genome [START_REF] Wang | Donkey genomes provide new insights into domestication and selection for coat color[END_REF] (an E. hemionus subspecies from Iran), (table S4), were trimmed with cutadapt (v1.18) (50) and aligned to the E. caballus reference genome (eqCab2.0) using the BWA (v0.7.17) (51) mem program. PCR duplicates were removed using Picard MarkDuplicates (v2.20.0) (52), and reads aligning to the reference genome with mapping quality score below 30 were removed using samtools 1.9 [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF].

SNPs selection:

We curated the 36 million biallelic variant list used to differentiate equids in the Zonkey workflow [START_REF] Schubert | Zonkey: A simple, accurate and sensitive pipeline to genetically identify equine F1-hybrids in archaeological assemblages[END_REF] to filter out variants found in repeated sequences using a EqCab2 genome repeat mask downloaded from the UCSC browser https://genome.ucsc.edu/cgi-bin/hgTables . The rationale for this filtration was that these variants would be less reliably called using the short, mitochondrial genomes had a mean coverage of 52x (1864 hemippe), 40x (1892 hemippe) and 51x (Göbekli Tepe GT64). Complete mitogenome sequences were generated by consensus calling of the bases using Geneious [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF]. To obtain a full-length mitogenome, gaps were filled using both targeted PCR data of the HVR [START_REF] Bennett | Taming the late Quaternary phylogeography of the Eurasiatic wild ass through ancient and modern DNA[END_REF] and by analyzing, at the boundary of the gaps, the soft clipped reads resulting from mapping with bwa mem rather than bwa aln.

Nuclear SNPs were called from the Umm-el-Marra and hemippe bam files using the samtools [START_REF] Li | A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data[END_REF] mpileup command with the following parameters: -B -A -Q20 and specifying only the 15.5 million SNP positions described above. Calling and selection of a single allele for all heterozygous sites were performed using pileupCaller (61). This resulted in 6.8 million shared positions between extant equids and the 1864 hemippe, 2.2 million shared with the 1892 hemippe, 10.9 million with the Göbekli Tepe GT64 sample, and 45,604 with the Umm-el-Marra 9 sample (table S5).

Principal component analysis (PCA) was performed using EIGENSOFT SmartPCA (v16000) [START_REF] Price | Principal components analysis corrects for stratification in genome-wide association studies[END_REF][START_REF] Patterson | Population Structure and Eigenanalysis[END_REF] by projecting the samples with partial coverage onto eigenvectors calculated from all shared positions of well covered equids (projectlsq: YES). For the PCA represented in Figs. 4A andS4 D-E, we used the 15.2K SNPs shared between UMM9 and both the GT64 and the 1864 hemippe, and only the 1892 hemippe was projected. For the PCA represented in Figs. S4 A-C, we used the 4.1M SNPs shared between the extant equids used and both the GT64 and the 1864 hemippe, with both the UMM9 and the 1892 hemippe being projected. Admixture (v1.3.0) (64) was used to estimate ancestry of the six horses, six donkeys, three Mongolian khulans, two kiangs, the onager, the 1864 and the 11,000 BP GT64 hemippe and UMM9 using the 15.2K SNPs shared between UMM9 and both the GT64 and the 1864 hemippe, and a four population model (K=4). Fig. 4B represents the admixture bar graph obtained in 70% of the 30 iterations (90% showed the UMM9 sample as a 1:1 admixture of onager/hemippe and donkey). The bifurcating tree with gene flow was performed using treemix [START_REF] Pickrell | Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data[END_REF] with the 40.4K SNPs shared between UMM9 and either of the three hemippes, and considering the following equid groups (number of individuals): Horse ( 6), Donkey [START_REF] Heimpel | Towards an understanding of the term SiKKum[END_REF], Hemippe (3), UMM9 (1), Onager (1), Mongolian khulan (3) and Kiang [START_REF] Ur | Cycles of Civilization in Northern Mesopotamia, 4400-2000 BC[END_REF]. Horses were used as the outgroup and sample size

Author contributions

EMG and JW initiated the project; EMG and TG conceptualized and supervised the project; JW, JP and GS provided material; EAB, WB, SC, EMG and TG performed the laboratory work; EAB, TG and EMG analyzed the data; EAB, TG and EMG wrote the paper with input from JP and GS.

Competing Interests. The authors declare that they have no competing financial interests that might have influenced the work described in this manuscript. Horses were used as the outgroup and sample size correction was disabled. The tree obtained with one gene flow event is represented. The residuals with no or one gene flow event are plotted in Fig. S5.

Legends to Figures

Supplementary Material