Shedding Light on the Dark Ages: Sketching Potential Trade Relationships in Early Medieval Romania through Mitochondrial DNA Analysis of Sheep Remains

: Southeast Europe has played an important role in shaping the genetic diversity of sheep due to its proximity to the Danubian route of transport from the Near East into Europe, as well as its possible role as a post-domestication migration network and long tradition of sheep breeding. The history of Romania and, in particular, the historical province of Dobruja, located on the shore of the Black Sea, has been influenced by its geographical position at the intersection between the great powers of the Near East and mainland Europe, with the Middle Ages being an especially animated time in terms of trade, migration, and conflict. In this study, we analyzed the mitochondrial control region of five sheep originating from the Capidava archaeological site (Dobruja, Southeast Romania), radiocarbon dated to the Early Middle Ages (5th–10th century AD), in order to better understand the genetic diversity of local sheep populations and human practices in relation to this particular livestock species. The analyses illustrate high haplotype diversity in local medieval sheep, as well as possible genetic continuity in the region. A A.A., C.Z.-T., E.G., D.P. O.G.; investigation, A.A., C.M., C.Z.-T., E.G., I.R., O.G., C.D., curation, A.A. writing—original preparation, A.A.; writing—review I.R., C.M., E.G., D.P., O.G., C.U., Z.K.P., A.P. and B.K.; visualization, A.A. and I.R.; A.P. and B.K.; project administration, I.R. B.K.; funding acquisition, A.A.,


Introduction
Sheep (Ovis aries) were domesticated approximately 11,000 years ago in the Fertile Crescent region [1] and are classified into five mitochondrial DNA (mtDNA) haplogroups (A, B, C, D, E) that diverge phylogenetically and probably descend from multiple wild Moufloniformic populations [2][3][4]. From the Near East, where they were first domesticated, they did not spread across Europe in a single migration episode, but rather in multiple, separate migration events related to human movements [5]. Three main colonization routes can be traced in Europe following domestication: the Mediterranean route, the Danubian route, and the Northern European route [6][7][8][9].
As a consequence of livestock migration routes into Europe, as well as of the long history of sheep breeding in the region, Southeast Europe was found to be a genetic diversity hotspot [10][11][12][13][14][15] and may have served as a post-domestication migration node [15]. As such, detailed knowledge of the genetic background of local sheep can be useful in guiding conservation efforts and tracing genetic resources [16,17].
The genetic make-up of sheep populations is defined by a complex history due to a continuous domestication process, shaped by introgression from wild species, strong gene flow mediated by human migrations, as well as trade and artificial selection practices [15,[18][19][20][21]. Therefore, the resulting genetic background will be a reflection of the culture manipulating it [22].
Mitochondrial DNA, a useful tool for tracing maternal origins and revealing the migration patterns of humans and the domestic animals they brought with them [16], has proved to be a suitable genetic marker for understanding the genetic history of sheep. Studies focused on ancient sheep mtDNA hypervariable region (HVR), particularly from the Copper Age [23] and Bronze Age [24][25][26][27][28][29] have assisted in unfolding the origin of domestic lineages in China, Russia, Turkey, the Alps, and Estonia, while those that discuss sequences dated to the Iron Age, Medieval, and Post Medieval Times have illustrated different selection and trade practices. For example, modern selective breeding practices did not change the haplogroup distribution of Italian sheep already established in the Middle Ages [30], while extensive research in the Northern European region has shown genetic continuity in Estonia for the past three millennia [28], no sign of population replacement in Finland since the Iron Age [8], as well as signs of expanding trade in the Middle Ages and strong artificial selection in the Modern Era [22,28,29]. Unfortunately, studies that aim to illustrate sheep mitochondrial genetic variation across different historical periods and in comparison to modern sheep, have so far been restricted to these few geographical areas.
The present paper is focused on investigating the mitochondrial genetic diversity of a small group of sheep consisting of five individuals originating from Early Middle Ages Dobruja (Capidava, Southeast Romania, 5th-10th centuries AD), a region underrepresented in published literature but of great importance both in terms of historical ties [31] and sheep genetic diversity [10][11][12][13][14][15]. Given the period of intense raids and increasing trade in the region throughout the 5th-6th century, as well as the rising power of the First Bulgarian Empire, which occupied Dobruja from the late 7th century and its documented ambivalent relationship with the Byzantine Empire-alternating between times of peace, alliance and lucrative trade and times of war [31]-this paper also aims to reveal the genetic relationships between the Capidava cohort and other published ancient and modern sequences.

Sample Collection, Archaeological Context, and Radiocarbon Dating
The samples (Ovis_Cap) consisted of a total of 7 sheep mandibles with embedded teeth found in the archaeological fill excavated from the Capidava archaeological site in South-East Romania during the campaigns of 2010-2014, representing the entire ovine material found in that respective fill. Initial morphological inspection indicated that the remains belonged to domestic sheep, and prior to DNA and collagen extraction, the bones were safely stored in a dedicated storage room.
The archaeological campaigns focused on the retrieval of human skeletal remains led to the recovery of individuals dated through grave goods, funerary rites, and radiocarbon analysis as part of a 10th century AD, probably Christian, population [32,33]. However, the same clear stratigraphic context could not be inferred for the livestock remains uncovered, and, as such, radiocarbon dating was performed for all 7 Ovis_Cap samples.
Prior to collagen extraction, the bone was tested by elemental analysis on a Thermo Scientific FlashEA 1112 Series Elemental Analyzer to determine if the samples were suited for Accelerator Mass Spectrometry (AMS) dating. Generally, a ratio of C:N ranging from 2.9-3.5 was considered appropriate, according to Harvey et al. 2016 [34].
Bone collagen extraction was performed using an adapted acid-alkaline-acid technique, according to [35,36]. The bone was ground to a powder and homogenized using liquid nitrogen. About 0.7-1.2 g of bone powder was transferred to 15 mL tubes, 10 mL of 1N HCl were added, and the samples were let to rest for 30 min. The samples were centrifuged for 5 min at 3000 rpm, and the supernatant was removed. Next, the samples were brought to pH 7 by washing with distilled water (DW). The washed pellet was resuspended in 7 mL of 0.1% NaOH for 30 min, and then the samples were centrifuged for 7 min at 3000 rpm and the supernatant removed. Next, the samples were brought to pH 7 by washing with DW. Atmospheric CO2 absorbed during alkali treatment was removed by adding 7 mL of 1N HCl for 30 min and centrifuging 3 min at 3000 rpm to facilitate sedimentation. The pellet was washed until the supernatant reached pH 3. A few milliliters of solution were left over the pellet and gelatinization was accomplished by incubating the samples at pH 3, 56 °C, for 16 h. Following gelatinization, the samples were centrifuged at 2500 rpm for 3 min for better separation, filtered with 45-90 µm Ezee-filters (Elkay, Hampshire, UK) to remove insoluble residue, and lastly ultrafiltered by centrifuging in Vivaspin 15 30 KDa MWCO ultrafilters (Sartorious, Epsom, UK), until ~0.5-1 mL volume remained in the ultrafilter.
The radiocarbon measurements were carried out at the 1 MV AMS facility at the RoAMS laboratory from IFIN-HH [37].
Collagen was kept at −4°C before being freeze-dried at −45°C overnight. After this, step samples were weighed into tin capsules and burned in an Elemental Analyzer (EA, varioMicro Cube, Elementar, Langenselbold, Germany). The resulting CO2 from sample combustion was adsorbed on the zeolite trap of the Automated Graphitization Equipment (AGE, ETH Zurich, Switzerland). Finally, the pure CO2 was thermally released into an AGE reactor with an iron catalyst on the bottom, where it was reduced to graphite in the presence of excess hydrogen. The final product was represented by a homogeneous mixture of carbon and iron. The entire process was controlled by the PC software Lab VIEW [38]. In the end, all of the prepared targets, containing approximately 1 mg of C and 5 mg of Fe, were pressed in Al cathodes and inserted in the negative ion source carousel.
The samples were measured together with the standard and blank materials for 50 min, consisting of 10 individual measurements divided into 10 blocks of 30 s each. Oxalic acid II (OXII [39]) was used as the modern standard, obtaining 134.076 ± 0.607 pMC (determined as the average and standard deviation of 5 cathodes). Anthracite coal was used as a chemical blank from which resulted the process background of 0.267 ± 0.038 pMC, corresponding to a radiocarbon age of 47.659 ± 1222 yr BP. Data analysis and radiocarbon ages calculation were performed using the BATS software, developed at ETH Zurich [40].
Conversions from radiocarbon ages (yr BP) to calendar ages (yr calBP) were made using OxCal version 4.4 with IntCal20 [41,42] and were quoted as cal. AD.

DNA Extraction, PCR Amplification and mtDNA Sequencing
All steps were performed in sterile laboratory environments dedicated to ancient DNA research at the Institute for Interdisciplinary Research in Bio-Nano-Sciences, "Babeș-Bolyai" University, Romania, following standard guidelines and proper workflow protocols [43]. During the entire DNA extraction procedure, care was taken to minimize contamination from external DNA sources. Firstly, teeth were physically extracted from mandibles, their external layer was mechanically removed using a dental drill and were UVirradiated for 40 min to eliminate surface contaminants. Next, the dentin powder was obtained using sterile dental drills. Approximately 30 mg of dentin powder was obtained for each sample.
The decalcification and enzymatic digestion procedure involved adding over the obtained powder 1 mL of Digestion Solution (EDTA 0.5M pH 8, SDS 0.5%, TRIS-HCl 50 mM, pH 8) and 5 µl of Proteinase K (Jena Bioscience, Jena, Germany) at a final concentration of 0.1 mg/mL and incubating that overnight (approximately 16 h), with agitation, at 56 °C [44]. Following its preparation, the Digestion Solution was subjected to UV light treatment for 45 min, as were all other solutions made in-house [45].
DNA extraction was performed using an adapted protocol consisting of 2 main steps. It combined a silica-based column extraction step [46] with the washing, drying, and elution steps according to the Nucleospin Tissue DNA extraction Kit (Macherey-Nagel, Düren , Germany).
During the first step, the samples were centrifuged for 5 min at 2000× g, transferred to a clean 2 mL tube and then centrifuged again for 5 min at 2400× g. The supernatant was divided into ~200 µl aliquotes (in 2 mL tubes) and 5 volumes (~1 mL) of Binding Buffer (Guanidine hydrochloride 5 M, Isopropanol 40%, Tween 20 Df = 10x 0.05%, deionized water) were added. Approximately 750 µL of the resulted sample were added to a silicabased column and the samples were centrifuged for 1 min at 2400× g. The flowthrough was discarded and this process was repeated until the entire volume of the sample had been loaded.
Next, the samples were washed with Washing Buffers BW and B5, dried and eluted according to the instructions provided in the Nucleospin Tissue DNA extraction Kit (Macherey-Nagel, Düren, Germany).
Negative extraction and amplification controls were performed for every sample to screen for potential contamination.
A 598 base pairs (bp) sequence of the Ovis aries mtDNA hypervariable region was targeted ([GenBank: NC001941] positions 15,978-16,501) using a set of 5 overlapping primers as described by [8], but using adjusted annealing temperatures. In order to test for the presence of the Y-chromosome single nucleotide polymorphism (SNP) marker, a part of the 5′-promoter region of the SRY gene was also analyzed-a 130-bp sequence, spanning positions 58-178 of the Ovis aries Y-chromosome [GenBank: AY604734] [8].
Successfully amplified target sequences were cloned (6 clones/PCR product) using the pJET1.2/blunt cloning kit (Thermo Fisher Scientific, Waltham, USA). DNA sequences were obtained by Sanger sequencing (Macrogen Europe, Amsterdam, The Netherlands), and the resulting consensus sequences were deposited under GenBank Accession Numbers MW928520-MW928524.

Sequence and Population Genetics Analyses
In order to analyze the degree of mtDNA variation of the Capidava cohort, comprising of 5 samples, a summary statistic analysis was performed on a 529 bp region (positions 16027-16556), namely the number of segregating sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and average number of nucleotide differences (k). To better understand the result in a broader geographical and historical context, we performed the same analysis on 18 more sets, comprising 3 to 7 sequences of approximately 530 bp in length originating from the same country and historical period, from a compiled database of sequences deposited in international repositories. The results of the nucleotide diversity values were plotted against the haplotype diversity values. The detailed results of the summary statistics and sequences employed in the analysis are available in Table S1.
Using a comparative dataset of ancient [8,27,28,30], modern [3,8,[47][48][49][50], and reference sequences [4,49] previously deposited in international databases (Table S2), the following analyses were performed: 1. In order to investigate the genetic relationship between the Capidava cohort and other ancient sequences representative for different geographic regions and historical periods, as well as to estimate the degree of admixture and determine potential trade routes and migration patterns, a Median Joining Network [51] was constructed and pairwise genetic distances were calculated using the Ovis_Cap sequences and 41 sequences available in the international databases (GenBank Acc. No. available in Table  S2). All sequences were aligned and trimmed a 304 bp region between the nucleotide positions 16099 to 16402 (NC_001941) to allow for a homogenous set. 2. To estimate the degree of genetic continuity and similarity to modern-day breeds, a Median Joining Network [51] was constructed, and pairwise genetic distances were calculated on the Ovis_Cap sequences and 37 and 59 sequences, respectively, representing a variety of breeds originating from different countries (GenBank Acc. No. available in Table S2). All sequences were aligned and trimmed to a 474 bp region, between the nucleotide positions 16027 to 16500 (NC_001941). 3. As a means to analyze whether the Capidava cohort shared a higher percentage of haplotypes with the ancient or modern data set, a shared haplotype analysis (SHA) was performed by counting the absolute and relative number of shared haplotypes between the Ovis_Cap sequences and 88 other sequences, combining the previously defined ancient and modern data sets (GenBank Acc. No. available in Table S2). All sequences were trimmed to a 304 bp region, between the nucleotide positions 16092 to 16396 (NC_001941).
Median-joining networks [51] were constructed using PopArt v. 1.7. [52] with ε = 0. Pairwise genetic differences (Fst) were determined, and the shared haplotype analysis was performed in the Arlequin v. 3.5.2. software [53]. The model parameters used for pairwise genetic distance calculation and shared haplotype analysis were determined using the jModelTest 2 software [54,55]. For both the ancient and modern datasets used for comparative analysis, as well as for the combined ancient and modern set, the resulting model was  with a Gamma-value of 0.032. Heatmaps of pairwise genetic differences, shared haplotypes, as well as the scatterplot of haplotype and nucleotide diversity of the comparative dataset were visualized using the Rmcdr, lattice, grid, and RColorBrewer packages of the R v. 4.0.3 statistical environment [57]. Geographical visualization of pairwise genetic differences was achieved using the QGIS v. 3.16 software [58].

Radiocarbon Dating, PCR Amplification and Identified SNPs
All seven samples, previously determined to be suited for AMS analysis, were successfully radiocarbon-dated. Genetic analysis was successful for five out of the seven samples. The entire mtDNA target sequence (598 bp) was reconstructed for one sample, whereas for four other samples, a slightly shorter (529 bp) but still informative HVR fragment was obtained. Amplification of the 130 bp region of the 5′-promoter of the SRY gene was inconclusive, as was the amplification of mtDNA from the oldest sample (Ovis_Cap_7) and one of the samples originating from the 5th-6th century AD (Ovis_Cap_3).
Sequence analysis led to the identification of five different unique mitochondrial haplotypes and by comparing them to the reference sequence (NC_001941), 13 polymorphic sites were found, both indicators of a diverse cohort. Besides mutational motifs, the alignment of multiple clones from the same PCR product generally showed the presence of miscoding lesions, mostly in the form of C-T and G-A transitions, which occur due to post-mortem DNA damage, pattern representative of ancient sequences [59]. For all of the samples processed, both the negative extraction and amplification controls were clean, further supporting the authenticity of the results. Data concerning PCR amplification results, identified SNPs, as well as radiocarbon dating results for each sample, are presented in Table 1.

Summary Statistics Analysis
Due to the small sample sizes of the 18 compiled sequence sets (chosen as such in order to eliminate potential bias), the majority of sample sets have a haplotype diversity value of 1, however certain clusters can be distinguished based on the combined values of Hd and π (Figure 1). The nucleotide diversity value of the Capidava cohort was relatively low compared to the modern sequence sets and the Medieval and Post Medieval Finnish sequences, however, it was higher than that of the Estonian sequences originating from both the Late Bronze/ Iron Age and Late to Post Medieval times, as well as the Finnish Iron Age set. The Capidava sample set has a Hd value of 1 and a π value of 0.00644, thus forming a cluster with the modern set from Romania and Hungary, the modern set from the region of Karelia and the Estonian sequences of Late to Post Medieval times.
These results served as a first indication of similarity between the Capidava cohort and historical Nordic sequences, as well as potential genetic continuity in the current territory of Romania.

Haplotype Similarity Analysis
The level of haplotype similarity between the sequences was assessed by performing a Median Joining Network (Figure 2). The result illustrates the fact that the entire Capidava cohort belongs to Haplogroup B, in accordance with the reported geographical distribution of haplogroups [9,18,21,24], and has a higher degree of haplotype similarity with sequences from Northern Europe, namely Finland and Estonia. There was no specific grouping according to historical era in this case, with the Ovis_Cap sequences showing haplotype similarity with Finnish and Estonian sequences from the Late Bronze Age to Post Medieval periods. A lower degree of similarity was indicated between the Capidava cohort and the Medieval Italian dataset, with only 2 out of the 5 sequences sharing the same haplotype for the analyzed region. The group that displayed the lowest similarity was the one representative of the Altai region, with sequences forming part of the Haplogroup A cluster.
While Ovis_Cap_1 and Ovis_Cap_6 have the same haplotype as 22 other sequences present in the database, Ovis_Cap_4, Ovis_Cap_5, and Ovis_Cap_2 have unique haplotypes, residing one and two mutation steps away from the main haplotype, respectively.

Genetic Distance Estimates
Pairwise genetic distances were calculated in order to estimate the degree of admixture between the sequences, and the results of Slatkin's linearized Fst were displayed in heatmap form, with the lowest value indicating the shortest genetic distance (Figure 3). The exact Slatkin's Fst values and their corresponding p-values are available in Table S3. The results of the pairwise genetic distances calculations were in agreement with those of the MJ Network and allowed for differentiation according to the historical period. As such, the closest sequence set to the Capidava cohort was that representative of the Esto-

Haplotype Similarity Analysis
The degree of haplotype similarity between the populations was assessed through a Median Joining Network (Figure 4). For this analysis, only the haplogroup B sequences (37) of the data sets were chosen in order to better illustrate the position of the Capidava set compared to modern sequences belonging to the same haplogroup and facilitate result visualization. The MJ Network containing both the haplogroup B and haplogroup A sequences was, however, available as Figure S1.
The first important aspect illustrated by the network was the fact that all five Ovis_Cap sequences were at least one mutation step away from all of the modern sequences, sharing no identical haplotypes with current breeds. Thus, the resulting network position was representative of an ancient, diverse cohort.
Moreover, the Capidava sheep display similarity with breeds from Romania and Hungary, Turkey, the UK, Poland, the Russian province of Karelia and Finland, in accordance with the results of the ancient dataset analysis and further suggesting possible genetic continuity in the region of Romania. As such, Ovis_Cap_6 has two additional polymorphisms compared to the Turkish breed Karayaka, Ovis_Cap_4 and Ovis_Cap_5 were one mutation step away from the main haplotype shared by the Turkish Tuj, Finnish Landrace, British Oxford Down (Finnish subpopulation), Polish Olkuska and the Karelian Viena breed. Ovis_Cap_1 had one additional polymorphism compared to the Turkish Hemsin and Serbian Pramenka breeds, as well as to the main haplotype and the Gyimesi Racka/Turcana breed of Romanian and Hungarian origin. Ovis_Cap_2 lies two mutation steps away from Ovis_Cap_5.

Genetic Distance Estimates with Focus on Different Geographic Regions
While the MJ Network illustrates the level of haplotype similarity between the Capidava sequences and individual breeds, the results of the pairwise genetic distance analysis offer a broader picture of the degree of genetic distance, with values indicating closeness to a population of diverse breeds from a particular geographic region. The exact Slatkin's linearized Fst values and their corresponding p-values are available in Table S4. Below, in Figure 5, they are displayed in heatmap form, with the lowest value indicating the shortest genetic distance. Overall, they illustrate which countries have modern sheep populations that more closely resemble the Capidava population.
As such, the results show that the populations from Romania and Hungary, the UK, and the Netherlands were the closest to the Capidava cohort (Slatkin's Fst = 0), followed closely by the Karelian population (Slatkin's Fst = 0.14991). The Turkish, Finnish, Polish, and German populations displayed similar degrees of genetic distance, with values ranging from 0.4 to 0.5, while the Chinese and Central Asian populations showed the highest degree of genetic distance.
Interestingly, there was a relatively high level of genetic distance between the Capidava and South Caucasus populations, which may not have been expected given the regions' geographical proximity and access to the Black Sea. However, a comparable pattern was revealed by [11], who showed, based on nuclear markers, that the Caucasian Bozakh breed did not cluster with the Southern European Zackel breeds.

Analysis of Shared Haplotypes between the Capidava, Ancient, and Modern Datasets
When comparing the relative degree of haplotype sharing between the Capidava, ancient and modern datasets, a clear distinction can be made, with the ancient sets displaying higher levels of haplotype sharing with the Capidava set. The results of the relative number of shared haplotypes are presented below (Figure 6), in heatmap form, while the absolute numbers are given in Table S5.
Out of the ancient sequences, a historical, as well as a geographical distinction, can again be made, with the Finnish Iron Age, Estonian Iron Age, and Estonian Late and Post Medieval sequences showing the highest level of haplotype sharing, followed by the Estonian Medieval set (80%) and only then by the Italian Medieval set (60%). These results reinforced the similarity between the Ovis_Cap and Nordic, over Italian, sequences, but provide a different angle regarding which Nordic group shows the closest genetic affinity to the Capidava cohort than the results of the pairwise genetic difference analysis.
Out of the modern sequences, the results of the MJ Network ( Figure 4) were supported, with the Karelian, Turkish, and Serbian sets displaying the highest level of haplotype sharing. Interestingly, the Romanian/Hungarian set showed one of the lowest levels of haplotype sharing (14.3%). What was also noteworthy was the fact that the Karelian set shared a higher percentage of its haplotypes with the Capidava set (83.3%) than the Estonian Medieval set.
When referring strictly to the Capidava cohort, 3 out of the 5 Capidava haplotypes, meaning 60%, were specific to the cohort, further indicating a diverse set which, while showing similarity to other cohorts, maintained a degree of uniqueness.

Discussion
Taken together, these results illustrate the relationship between the Capidava group and other sequences representative for different geographic regions and historical periods, providing insight into the degree of admixture and serving as a starting point for further investigation into potential trade routes and migration patterns.
As such, when compared to other ancient datasets, the Capidava sequences show the highest affinity to Nordic sequences. Similar patterns were identified in a previous study of mtDNA variation in ancient Finnish sheep, which showed close genetic connections with Southern and Eastern European breeds [8]. However, no clear, universal distinction between historical periods can be made, which may be a consequence of the observed continuity and consistency of genetic diversity through time in maternal lineages in the region [8,28].
The fact that the Capidava cohort is more tightly associated to Nordic, rather than Italian and shows no affinity to Altai sequences, indicates a possible tendency for North to South sheep trade and migration routes, rather than East to West (Figure 7). Yet, given the limitations of using a small sample size and only a partial mtDNA marker, this can only be treated as a hypothesis and would need to be further confirmed through genomewide assessments performed on larger sample sizes. Next, despite the fact that the Altai sheep populations are rather diverse, possibly due to strong cultural interactions in the region [27], the high degree of genetic distance between the Capidava and Altai populations serves as a prospective indicator of a lack of interaction along the Asian route at that time. Again, however, a potential Asian interaction may not have been revealed in this context due to the low number of ancient Ovis samples.
The results of the comparison with modern sequences facilitate the depiction of a slightly clearer image. The affinity of the Capidava group towards Nordic sequences is visible in the modern Finnish sequences as well, while the observed high level of genetic similarity with the sequences from the province of Karelia, in the Russian Federation, is in accordance with previous findings regarding the historical genetic influence of Russian breeds in those originating from Northern Europe [8,9]. Our findings are different than those of [11] in regards to the genetic affinity between Nordic and Southeast European breeds, however, this might be a consequence of the different genetic markers used. Nuclear SNPs reflect biparental ancestry and have higher resolution, while mitochondrial DNA only reflects maternal ancestry. As such, one possible explanation is that the close genetic distance and haplotype similarity between the Ovis_Cap samples and Nordic sequences might be the result of gene flow focused on ewes.
Another notable similarity is the one between the Capidava cohort and modern Turkish breeds, previously shown to have weak genetic structure, possibly as a consequence of intense human-mediated gene flow [60], and high haplotype diversity [2,26]. The entire cohort displays genetic similarity to the modern Turkish breeds, with Ovis_Cap_5 (9th-10th century AD) behaving no differently from three of the other samples (8th-9th century AD). The four samples originating from the 8th-10th century AD have one or two additional polymorphisms than the Turkish breeds. Ovis_Cap_2, however, dated to 5th-6th century AD, lies further away. Since the Ovis_Cap samples span the period between the 5th and 10th century AD, these results need to be analyzed from the perspective of Dobruja's relationship with various invading nomadic tribes, as well as with the Byzantine Empire, which at the time occupied territory belonging to modern-day Turkey. The 5th-6th century saw Dobruja under loose Byzantine control, plagued by invasions from tribes such as the Slavs and Avars, but also the site of blooming trade, particularly in the 6th century [31]. The historical province became part of the First Bulgarian Empire at the end of the 7th century AD [31] and this occupation appears to be visible in the Capidava cohort as well. The four samples dated to the time of the Empire, which had strong ties to the Byzantines, show higher similarity to native Turkish breeds than Ovis_Cap_2, which is highly polymorphic but less similar to the Turkish sequences. While the historical record portrays the relationship between the Bulgarian and Byzantine Empires as complex and alternating between periods of alliance and war [31], our results indicate the possibility of a simpler, more constant dynamic for the everyday man, with no significant changes in local trade relationships for the last two centuries of the millennium. However, the skewed representation of time periods, with three Ovis_Cap samples originating from the 8th-9th century and only one originating from the 9th-10th century, could be hiding potential changes as a result of the tumultuous relationship between the two empires.
The results of the comparison with the modern Chinese sequences complete those of the comparison with the ancient Altai sequences and support the hypothesized lack of interaction along the Asian route, with the Chinese dataset showing the highest level of genetic distance out of the modern datasets. Moreover, 4 out of the 5 Chinese breeds analyzed, namely Tong, Small Tail Han, Hu, and Mongolian, originated from Mongolian sheep [61] (cited by [25]), [62], and, when combined with the geographical position of the Altai mountains at the border between Siberia, Mongolia and China, a distinct absence of Mongolian influence in Capidava sheep can be noted. This corresponds with the dating results, which place the Capidava samples centuries before the Mongol Invasion of Europe [31].
Another noteworthy aspect is the fact that the change in local occupation is visible in the Capidava cohort, with the Ovis_Cap_2 sample being more genetically distant from the modern breeds and ancient sequences than the rest of the Capidava samples. The Ovis_Cap samples originating from two different time frames during the Bulgarian occupation (8-9th and 9th-10th AD) behave as part of the same, albeit diverse, group, and no difference in genetic closeness to both modern and ancient sequences can be observed. Ovis_Cap_5, even though it originates from a later period, shows the same pattern of affinity to ancient Nordic and modern Finnish, Turkish, Russian, and Polish as the other three sequences, which could indicate no major changes in the genetic make-up of local sheep at that time.
In order to further assess the possibility of genetic continuity in the region, the Ovis_Cap sequences were compared with sequences originating from both Romania and Hungary and belonging to the Gyimesi Racka and Turcana populations, considered to be one transboundary breed according to genetic data [47]. The sequences used for comparison belong to highly diverse populations of one of the most abundant and genetically diverse breed groups in Europe, the South-East European Zackel, namely the Pramenka subgroup [10,14,47]. In this case, the combined results of the MJ Network, pairwise genetic distance analysis and shared haplotype analysis paint an interesting picture. While the MJ Network shows no identical haplotypes between the ancient and modern populations, the ancient haplotype Ovis_Cap_1 has only one and two additional polymorphisms from two of the modern sequences, respectively. The results of the SHA indicate a similar tendency, with only 14.3% of the modern haplotypes being shared with the ancient population. The results of the pairwise genetic distance analysis, however, show that the Capidava cohort is closest to the Romanian/Hungarian dataset out of the modern datasets. As such, the absence of Ovis_Cap haplotypes in modern populations could suggest their possible loss due to artificial selection [8], while the close genetic distance between the ancient and modern datasets could suggest the fact that the process of artificial selection was based on local sheep populations.

Conclusions
This study provides new insight into the genetic history of sheep in Southeast Romania and facilitates the understanding of trade practices in Early Medieval Dobruja. The increased haplotype diversity of the Capidava cohort, coupled with a variable affinity towards sequences originating from different geographic regions, serve as a potential indicator of expanding sheep trade in Southeast Romania at that time.
The results of the comparisons with previously deposited ancient and modern sequences indicate a higher tendency for North to South sheep exchange, rather than East to West and serve as a potential confirmation of the documented lucrative trade between the First Bulgarian Empire and the Byzantine Empire. The samples cluster according to historical context, with the Bulgarian occupation being visible in the cohort.
The four samples dated to the 8th-10th century AD behave like one diverse cohort with no change in the way they cluster according to the period of origin (8th-9th vs. 9th-10th century AD), which indicates no significant changes in the genetic make-up of local sheep at that time. The comparison with modern Romanian/Hungarian sequences suggests that, while certain haplotypes may have been lost in time due to artificial selection, this process may have been based on ancient, local, populations.
However, as a consequence of the small size and the use of only a partial mtDNA marker, our study represents only a small part of the picture and can only put forward hypotheses regarding livestock treatment in historical Romania. Further research concerning larger sample sizes and spanning a larger time frame is needed in order to complete the picture of Romanian medieval trade practices.

Supplementary Materials:
The following are available online at www.mdpi.com/1424-2818/13/5/208/s1, Table S1: Summary statistics, detailed results, and sequences sets used; Table S2: GenBank Accession Numbers of sequences used in comparative analyses (ancient/modern/reference); Table S3: Slatkin's linearized Fst values and corresponding p-values for the ancient dataset; Table S4: Slatkin's linearized Fst values and corresponding p-values for the modern dataset; Table  S5: Shared haplotype analysis results; Figure S1: Median Joining Network of the entire modern dataset. Data Availability Statement: The mtDNA sequences presented in this study were submitted to GenBank (https://www.ncbi.nlm.nih.gov/genbank/) under the accession numbers MW928520-MW928524. All other relevant data are within the paper and its Supplementary Materials.

Acknowledgments:
This study has been supported by the College for Advanced Performance Studies, Babeș-Bolyai University, Cluj-Napoca, Romania.

Conflicts of Interest:
The authors declare no conflict of interest.