Did the historical range of the European bison (Bison bonasus L.) extend further south?—a new finding from the Yenikapı Metro and Marmaray excavation, Turkey

The origin of the European bison (Bison bonasus, Linnaeus, 1758) has been widely discussed and investigated in recent years. The species had a wide historic geographic distribution throughout the European continent during the middle and late Holocene, ranging from France in the west to the Caucasus in the east. However, archaeological evidence is needed to resolve the southern extent of the European bison distribution. We discovered one bison skull fragment during archaeological excavations in 2008 in the area of Yenikapı Metro and Marmaray (Turkey). Radiocarbon dating indicated the skull was deposited during the Byzantine period (seventh to eighth century AD). Mitochondrial genome analyses provided clear evidence that the skull was from a European bison. This is the first unambiguous evidence of the presence of this species in southeastern Europe during Byzantine times, which validates the historical written records of a potentially wider range of the European bison in historical times.


Introduction
During the early Holocene, European bison (Bison bonasus) were widespread throughout the European continent, but did not occupy Scandinavia. B. bonasus appeared in Denmark and Sweden only at the end of the last glaciation period (c. 10,000 BP) (Pucek 2004). During the middle and late Holocene (last 6000 years), the population of the European bison was widely distributed from France in the west to Russia's Volga River and the Caucasus in the east (Benecke 2005;Pucek 2004). It was proposed that bison were also present in other regions (modern Armenia, Azerbaijan, Georgia, Northern Iran and Asiatic parts of Russia) although evidence is needed to substantiate this hypothesis (Pucek 2004;Flint et al. 2002). About 100 years ago, B. bonasus reached an extinction threshold due to habitat degradation and fragmentation linked to agricultural and logging activities, as well as unlimited hunting and heavy poaching (Pucek 2004). Fortunately, extinction was prevented through several conservation initiatives and controlled breeding in recent years, which resulted in a population increase (Pucek 1986(Pucek , 2004. Information on the presence and distribution of the B. bonasus in historical times has been mainly based on ancient written sources (Benecke 2005). However, names used to describe European bison have varied wildly, which creates some confusion and weakens the reliability of historical written sources (Benecke 2005). For example, Aristotle used the word Bbonasus^to describe European bison in his reports on the history of animals (Aristotle and Balme 1991). In later years, between 1450 and 1850, Bauerochs^was used as the official name of wisent/bison/żubr. There were attempts to eliminate this confusion in subsequent accounts by assigning the names Baerox^for aurochs and Bbisont^for wisent (Ahrens 1921). Yet, some twentieth century and medieval descriptions use the name Baurochs^for European bison (Avebury 1913;Von Lengerken 1953).
Natural habitats of B. bonasus include large deciduous forests, bushes, shrubs, meadows and marsh areas (Benecke 2005;Pucek 2004). Consequently, it is expected that European bison could have occupied an area from the northern European plains to interior Asia based on the suitability of niches (Filean 2006). In fact, archaeological evidence indicates the widespread presence of European bison in Europe at least during the sixteenth century (Aleksandrowicz 1999). However, interpretation of the archaeological record is problematic: domestic and wild cattle remains could easily be incorrectly identified as European bison because of their similar size and morphological features (Benecke 2005). In particular, domestic male cattle and wild female cattle were possibly described as bison in some written sources based on their basic features (Filean 2006). For example, a bison bone was allegedly reported in faunal remains from Sos Höyük (on the border of Turkey and Iran), although the finding was classified as Bdoubtful^ (Howell-Meurs 2001). Similarly, two bones were described as bison remains in the faunal assemblage of Pınarbaşı, Central Anatolia, but the accuracy of the identification of these remains is also doubtful because of their possibility of being mixed with other bovid species (Carruthers 2005).
A bison skull fragment was found in Yenikapı Metro and Marmaray (Turkey) excavation in 2008. This was one of the most important findings from the Byzantine period found at this site. The 58-km 2 excavation area is located in the Yenikapı area of Istanbul, about 1.5 km from the Marmara Sea (Başaran et al. 2007) (Fig. 1).
Relatively large number of animal species have been found during subsequent excavations (Onar et al. 2013a, b). Amongst these, one identified European bison skull fragment was excavated during the 2008 season. Geographically, the  (Benecke 2005). This paper reports a unique finding of the bison remains from the Istanbul site is of utmost importance in relation to the addition of new geographical ancient distribution range.

Material and methods
A single skull fragment identified as European bison, from level −4.00 of trench 1Hc2 from the 2008 excavation season of Istanbul Yenikapi Metro and Marmaray excavation area, was used in this study (Fig. 2). The skull fragment was sent   Shapiro et al. (2004) to the Oxford Radiocarbon Accelerator Unit, Oxford University, for radiocarbon dating.
Ancient DNA work Samples were sent to the Australian Centre for Ancient DNA (University of Adelaide, Australia) where sample preparation and DNA extraction were performed in an off-campus, purpose-built, physically isolated, ancient DNA facility equipped with positive air pressure. The laboratory follows strict clean-room conditions to minimise the impact of DNA contamination. All rooms and workbenches are routinely decontaminated with a 3% sodium hypochlorite solution (bleach) and UV irradiated overnight. Consumable packages are discarded where possible or extensively wiped with 3% bleach and UV irradiated before entering the facility. Workers use sterile and disposable hooded coveralls, face masks, dedicated shoes and two pairs of gloves. A second laboratory on campus accommodates all post-DNA-amplification work. A strict one-way movement of materials and personnel is implemented, from pre-to post-amplification laboratory. Methodological precautions include total removal of sample surface, repeated irradiation with UV light, mock extractions, DNA amplification reaction blanks and replication of DNA amplification.
Sample preparation and DNA extraction Potential DNA contamination on the surface of the bison skull fragment was reduced by UV irradiation (260 nm) of all surfaces for 15 min, thorough wiping with 3% bleach and removal of ∼1 mm of the sample surface by abrasion using a dremel tool and a disposable carborundum disk. The sample (200 mg) was powdered using a Mikrodismembrator (Sartorius) and immediately decalcified and lysed overnight in 4 mL 0.5 M EDTA pH 8.0, 200 μL SDS and 40 μL proteinase K at 55°C under constant rotation. After lysis, DNA was purified using an in-solution silica-based method previously described (Brotherton et al. 2013) and resuspended in 200 μL of 10 mM Tris-1 mM EDTA buffer.

Sequencing of the mitochondrial control region
Polymerase chain reaction (PCR) was prepared in the ancient DNA laboratory. To test if the sample contained bison DNA, the mitochondrial DNA D-loop was amplified in six overlapping fragments. PCR was performed using 0.4 mM of each primer (see Table 1 for primer sequences), 2.5 mM MgSO 4 , 0.25 mM dNTPs, 1 mg/mL rabbit serum albumin (Sigma, fraction V) and 1.25 U Platinum Taq Hi-Fi polymerase (Invitrogen). Thermal cycling parameters were 94°C for 2 min, 50 cycles of 94°C for 20 s, 58°C for 30 s and 68°C for 30 s, then 68°C for 10 min. Resulting PCR products were visualised on a pre-stained 3.5% agarose gel under UV light. For successful PCR, amplicons were purified using Agencourt Ampure (Beckman Coulter) according to the manufacturer's protocol. Sanger sequencing was performed at the Australian Genome Research Facility Ltd (AGRF). For successfully sequenced amplicons, forward and reverse sequences were aligned to the bison control region, cleaned from primer sequences, and a consensus was made using the software Geneious v7 (Kearse et al. 2012).
Sequencing of the whole mitochondrial genome Doublestranded Illumina libraries were built from 20 μL of DNA extract using partial uracil-DNA glycosylase (UDG) treatment (Rohland et al. 2015) and truncated Illumina adapters with dual 7-mer internal barcodes, following the protocol from (Soubrier et al. 2016). The complete mitochondrial genome of the Marmaray sample was captured using hybridisation with commercial RNA baits (MYcroarray) (Soubrier et al. 2016). The enriched library was sequenced in paired-end reactions on an Illumina MiSeq machine. Sequencing reads were processed using the pipeline Paleomix v1.0.1 (Schubert  (Lindgreen 2012) was used to trim adapter sequences, merge the paired reads and discard reads shorter than 25 bp . BWA v0.6.2 was then used to map the processed reads to the reference mitochondrial genome of the European bison (GenBank ref. NC_014044). Minimum mapping quality was set at 25, seeding was disabled, the maximum number or fraction of gap opens was set to 2 and the edit distance was relaxed. MapDamage v2 (Jónsson et al. 2013) was used to verify the presence of expected contextual mapping and damage patterns considering the partial UDG treatment of the library (i.e. only the first and last base of sequencing reads should show higher levels of C-to-T and G-to-A transitions, respectively). Finally, the Marmaray consensus sequence was called using samtools/bcftools ) with minimum base quality at 30 and minimum depth of coverage at 3, and the Paleomix script Bvcf_to_fasta^. To place the studied individual within the modern bovid assemblage, a maximum-likelihood phylogenetic tree was constructed using the program PhyML v3 , HKY+G6 as the model of nucleotide substitution selected by comparison of Bayesian information criterion (BIC) scores in ModelGenerator v0.85 (Keane et al. 2006), NNI and SPR rearrangements to search for the tree topology and approximate likelihood-ratio tests to establish statistical support for the internal branches. The resulting tree, comparing the Marmaray skull DNA with DNA of other Bovini tribe representatives, is shown in Fig. 4.

Results
Radiocarbon dating The skull fragment belonged to the seventh-eighth century AD as determined by radiocarbon dating (OxA-32358- Table 2). Calibration of the radiocarbon date was performed using OxCal v4.2.4 and the IntCal13 curve (Reimer et al. 2013), and the calibrated age range curve is given in Fig. 3.
Ancient DNA analysis No amplicons were detected for the PCR amplification negative controls. Fragments of expected sizes were detected for amplicons 1, 4, 5 and 6, a larger fragment than expected was detected for amplicon 2, and amplicon 3 could not be generated. Sanger sequencing results showed that fragments 1, 4, 5 and 6 corresponded to the bison control region, whereas fragment 2 was identified as cattle, therefore considered as a potential contaminant. Based on these results, hybridisation capture and sequencing of the whole mitochondrial genome were attempted. A total of 25,486 sequencing reads mapped uniquely to the European bison mitochondrial reference sequence (NC_014044), covering 100% of the reference at an average depth of 123×. Damage pattern analysis showed 13.8% of C-to-T at 5′ ends and 7.1% of G-to-A at 3′ ends, with less than 0.5% substitution level at all other positions, as expected from an ancient sample and the partial UDG treatment of the DNA library. The phylogenetic analysis showed that the Marmaray specimen is nested within the B. bonasus clade, with very high confidence (Fig. 4).

Discussion and conclusions
The B. bonasus skull fragment found in Istanbul, at the meeting point of Asia and Europe, is the only remain of this species ever found in this region. Previous reports of European bison bones in Anatolian sites such as Sos Höyük (Howell-Meurs 2001) and Pınarbaşı (Carruthers 2005) are debatable because of the highly likely possibility of confounding some morphological features with other bovid species. We propose that DNA analysis of these other Anatolian remains is required to identify them accurately, such as in the present study. Thus, the nearest unequivocal presence of B. bonasus archaeological remains is in Bulgaria and Romania, in southeastern Europe (Benecke 2005). In Byzantine times, the Theodosius Harbour, being the second largest port of the Mediterranean Sea, was used as the most important centre for grain shipments from Egypt and remained a large harbour until the twelfth century (Kocabaş 2008(Kocabaş , 2015. The discovery of seventh-eighth century European bison remains in this harbour area suggests that commercial activities related to this species may have spread through Bulgaria and Romania to interior areas of Thrace. Bison may have been traded from the Balkans/Carpates as food supply for people living in the Constantinople area, although there was no sign of cut marks on the Marmaray bison skull fragment. This is the 57th animal species found in the Yenikapı Harbour area, which exhibits an extraordinarily diverse faunal assemblage from Byzantine times (Onar et al. 2013a, b). More importantly, the presence of a B. bonasus skull fragment in the Yenikapı Metro and Marmaray area strongly suggests that the historic distribution of European bison may have extended in the south as far as Asia Minor. This finding supports historical reports of the presence of European bison in the Constantinople area during the Byzantine period (Ahrens 1921).