Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence.

Domestic cattle were brought to Spain by early settlers and agricultural societies. Due to missing Neolithic sites in the Spanish region of Galicia, very little is known about this process in this region. We sampled 18 cattle subfossils from different ages and different mountain caves in Galicia, of which 11 were subject to sequencing of the mitochondrial genome and phylogenetic analysis, to provide insight into the introduction of cattle to this region. We detected high similarity between samples from different time periods and were able to compare the time frame of the first domesticated cattle in Galicia to data from the connecting region of Cantabria to show a plausible connection between the Neolithization of these two regions. Our data shows a close relationship of the early domesticated cattle of Galicia and modern cow breeds and gives a general insight into cattle phylogeny. We conclude that settlers migrated to this region of Spain from Europe and introduced common European breeds to Galicia.

Introduction

Neolithization is the transition of human societies from hunting and gathering to settlement and agriculture [1]. During this process, which started around 10,500 BP in the Middle East, selected plants and animals were domesticated by humans for food, clothing, or other resources these species provided [2]. Genetic data of ancient and modern samples from domestic cattle Bos taurus suggest that the domestication of its wild ancestor, the aurochs Bos primigenius, also began in the Middle East [3]. Domestic cattle spread from the Middle East to Anatolia, Europe and Africa, and local domestication or maternal interbreeding with B. primigenius seems to have had little effect on today’s gene pool [3, 4], although there is clear evidence that it occasionally happened [5-12].

Nowadays, all domestic cattle belong to four major mitochondrial haplogroups, T, Q, R, and I. Haplogroup Q originated together with T in the Middle East, while haplogroup R has so far only been found in Italy; it has been argued to be either the product of maternal interbreeding of European aurochs or to go back to a second independent domestication event [7, 13, 14]. However, a recent paleogenome study found an R haplotype in an almost 9,000 year old B. primigenius sample from Morrocco [14], leaving its origin open to discussion. A separate domestication event in southern Asia gave rise to zebu or indicine cattle, B. indicus, and this group is assigned to haplogroup I [15, 16]. Based on mitochondrial DNA data, a founding population of 29 to 783 females has been estimated for the beginning of the taurine cattle domestication in the Middle East [4, 17]. Following this bottleneck at the onset of the Neolithic, domestic stocks spread across East Africa, Anatolia, and Europe in co-migration with humans. Most of these taurine cattle nowadays belong to six sub-haplogroups of the mitochondrial haplogroup T [3]. Sub-haplogroups T1, T2, and T3 are derived from ancestral haplogroup T, and today sub-haplogroup T3 is dominant across the European mainland and Britain, while T1 mainly occurs in Africa [3]. T2 is currently common in Europe and the Middle East, but less frequent than T3 in Europe. The remaining sub-haplogroups have a much more narrow distribution, with T4 restricted to Asian breeds [6], T5 to Italian breeds [7], and T6 has only survived in a few Balkan breeds [9].

A region where the European mainland and Africa are in close proximity to each other is the Iberian Peninsula, to which pioneering farmers may have arrived from both the African and European coasts of the Mediterranean sea [18, 19]. Indeed, ancient domestic cattle are thought to have reached Spain at two separate times and following two separate routes. From the Near East, taurine cattle spread north-westward into Anatolia and south-westward into North Africa [14, 20]. From Anatolia, the early domestic cattle expanded into Europe following two different ways: One, associated with the LBK (Linearbandkeramik) culture and running along the Danube river, the other along the Mediterranean coast, associated with the Cardial culture and likely involving a smaller number of individuals transported by boat. In agreement with this hypothesis, there is a cline of genetic diversity from Eastern to Central-Western Europe, with the Iberian and Southern French cattle being the ones showing the lowest genetic diversity within Europe. Later on, African taurine cattle were probably introduced into Iberia from the Maghreb, from where they expanded into Southern Europe with a subsequent increase of the African haplogroup T1 in these areas [21]. Although the Iberian Peninsula has been at the center of studies trying to elucidate the times and ways of expansion of the Neolithic into Europe [22-24], most of the efforts have focused on the Northeastern and Mediterranean regions, while less is known for the Northwest or Atlantic Spain [18]. In the Cantabrian region, which includes the northern coast, Neolithization is dated to 8,000-4,000 BP, and there it is the result of a complex mosaic process in which hunters and gatherers lived along with settlers [25]. The first domestic animals of this region date to around 7,000 BP and were mostly ovicaprines, which were later followed by pigs and cattle [1]. However, much less is known about the region adjacent to Cantabria, Galicia, because no Neolithic sites have been identified in this region yet [26]. The oldest genetic information for Galicia is from the Chalcolithic period. It dates to 5,950-5,050 years BP and therefore to a time period in which settlers with domestic livestock were much more common than other forms of human societies [26]. To address this gap of knowledge and gain insights into the process of Neolithization in Galicia, we analyzed the mitochondrial DNA from 18 samples of cattle mostly found in mountain caves in this region, of which eleven yielded enough DNA for phylogenetic analysis. Their ages range from before the Neolithic up to recent (last 60 years), and they include the domestic cattle B. taurus as well as its pre-Neolithic wild ancestor B. primigenius. We used these data to retrace the genetic development of cattle in this area and compare them to their European, African, middle Eastern, and Asian relatives.

Materials and methods

Bone samples were collected from nine sites in Galicia, which are shown in Fig 1. Two of them were archaeological sites, but most of the bones were found in mountain caves, that were natural traps for the cattle. We tentatively assigned them to species based on their morphology and determined the age of the bones by direct 14C dating or by 14C dating of other bones found in the same cave (stratigraphic). All specimens are housed in the “Instituto Universitario de Xeoloxía, Universidade da Coruña”, the director of which is one of the co-authors (Aurora Maria Grandal D’Anglade). Therefore, no permits were required. For all samples, DNA was extracted from bone powder following a published protocol [27]. Genetic libraries were prepared from 10 μl (for the samples from Cobas, A Rocha Forte, Minas Rubiais, and A Tara) or 20 μl (for the remaining samples) extract following a single-stranded library preparation protocol [28]. Uracil-DNA glycosylase (UDG) was used to remove uracil. To enrich the resulting libraries for mitochondrial DNA, hybridization capture was performed according to a solid-state capture protocol, using the Agilent SureSelect 244k MicroArray. The baits on the array were designed based on the complete mitochondrial genome from an aurochs (B. primigenius; GenBank Acc Nr NC_013996.1) using 60 bp bait length and 3 bp tiling space. Hybridization capture was carried out in two consecutive capture rounds according to a published protocol [29]. The resulting product was then sequenced on the Illumina NextSeq 500 platform using 75 bp paired-end chemistry.

Fig 1

Map of Galician mountain caves.

A. The black rectangle shows the location of Galicia within the Iberian Peninsula. B. Map of Galicia, corresponding to the black rectangle region in A. C. More detailed map of the region corresponding to the black rectangle in B); it shows the location of caves in which bone samples described in this study were found; stars denote cave locations. (The map was originally created by author Aurora Grandal-d’Anglade specifically for this study.).

Sequence reads were adapter-trimmed and quality-filtered at a phred quality score cutoff of 28 using cutadapt-1.12 [30], forward and reverse reads were merged with flash-v-1.2.10 [31] were possible, and reads below a length of 30 bp were excluded. Mapping of merged reads against a B. taurus mitochondrial reference genome (GenBank Acc Nr.: AY526085) was done with the program bwa-v0.7.8 aln with default parameters [32]. We filtered out all alignments with a quality score below 30 using samtools-v0.1.19 (program: view) [33]. The same samtools version was used to remove duplicates (program: rmdup) and to collect statistics about the alignments (programs: index, idxstats, and depth). We used mapDamage 2.0.7 to assess damage patterns of the ancient DNA samples [34] and computed consensus sequences based on a minimum depth of 3 reads with angsd-v0.914 using the majority rule [35]. Consensus sequences were checked manually. If possible, samples were assigned to haplogroups based on diagnostic single nucleotide polymorphisms (SNPs) in the D-Loop region of the mitochondrial genome of B. taurus as they are described in previous studies [3, 6, 9, 36].

With the consensus sequences, we generated three multiple sequence alignments using MAFFT v7.310 and its FFT-NS-i algorithm with 1000 iterations [37]. All alignments were checked by eye for any obvious misalignments. First, we generated a multiple sequence alignment that contained data from most common cattle haplogroups and included all of our samples. As references, we used mitochondrial genome sequences of the species Bos taurus, Bos indicus, Bos grunniens and Bos primigenius from a published study [38]. The complete dataset then covered five cattle haplogroups (T, R, I, P, Q) and five sub-haplogroups of T (T1, T2, T3, T4, T5). The resulting alignment was used for a Bayesian phylogenetic analysis with the program BEAST v1.10.4 [39]. For this, we transferred the annotation of the reference mitochondrial genome EU177854 to our alignment, and starts and ends of the features were moved according to gaps in the alignment. Then, we predicted the best models of molecular evolution for each of the features in the annotation with the program partitionfinder-2.1.1 [40]. The analysis with BEAST was done using the best model predicted for each partition, a strict molecular clock, and a coalescent model, assuming a constant population size. It ran with a chain length of 100 million and a sampling every 1,000 generations. 10% of the generations were used as a burn-in, and the convergence of the analysis was checked using Tracer v1.7.1 [41].

Next, two species-specific alignments were generated. For B. taurus, we included our eight newly sequenced B. taurus samples and all publicly available sequences for which we could determine the haplotype based on the same eight diagnostic SNPs used to assign our samples. The B. primigenius-specific alignment contained our samples as well as all publicly available sequences of European wild aurochs. A full list of public data used in our study is given in S3 Table. The two species-specific multiple sequence alignments were used to generate statistical parsimony networks using the R-script TempNet [42]. Networks we generated were not temporally separated, instead they separated our samples from publicly available data.

Results

We were able to generate consensus sequences from 11 (out of the 18) bone samples that were based on a coverage between 3.34x and 72.93x. Those 11 were all found in caves. The remaining seven samples were discarded due to low mapping coverage and sequencing quality. The 11 consensus sequences covered between 44.88% (7,333 bp) and 99.99% (16,337 bp) of the cattle reference mitochondrial genome (GenBank accession number: AY526085.1). An overview over dating, mapping and SNP analysis is shown in Table 1. More detailed dating and mapping results, as well as a figure containing SNP positions are in S1 and S2 Tables, and in S1 Fig, respectively. All raw sequencing data is deposited in the NCBI SRA under BioProject ID PRJNA705960. Addtionally, the partial mitochondrial genomes are made public in NCBI GenBank under accession numbers MW689247-MW689255, except mitochondrial genomes of samples Bt_MR_(3736BP) and Bt_CO_(60BP), which cannot be made available at NCBI GenBank, since they do not meet the GenBank submission criteria (less than 50% ambiguous characters). The analysis of ancient DNA patterns using mapDamage [34] shows terminal transitions indicating cytosine deamination ranged from 0.25 to 0.67% in mean per sample. These numbers are reduced due to UDG-treatment of the libraries. Mean fragment lengths ranged from 48.05 to 53.23 bp, consistent with an ancient origin of the sequences.

Table 1

Summarized dating and mapping results for 11 samples of Galician cattle subfossils.

Samples are referred to by their sample code in the text, which includes information about the name of the cave in which the subfossil was found, the dating, and the species assignment. Context refers to the historical context based on dating; results of subfossil dating are based on carbon dating and calibration. Mapping depth and and the percentage of the reference genome covered give an overview of the mapping results. Species assignments are based on morphological and genetic data. The haplotype was defined based on SNP analysis.

Sample	Cave	Context	Age cal BP (2-sigma)	Mapping depth	% of reference genome covered	Species	Haplotype
Bt_SA_(1200BP)_1	Saballeiros	High Middle Age	1125 ± 119	18.39	93.08	Bos taurus	T3
Bt_SA_(1290BP)_2	Saballeiros	High Middle Age	1233 ± 59	72.93	99.99	Bos taurus	(T3)
Bt_MR_(3736BP)	Mina Rubiais	Chalcolithic	4082 ± 63	6.49	44.88	Bos taurus	?
Bt_AR_(1200BP)	Arcoia	High Middle Age	1143 ± 88	30.95	98.75	Bos taurus	T3
Bt_AR_(2460BP)	Arcoia	Iron Age	2535 ± 175	33.36	99.68	Bos taurus	T3
Bp_CL_(9100BP)_1	Chan do Lindeiro	Mesolithic	9295 ± 170	10.07	66.89	Bos primigenius	P
Bp_CL_(9100BP)_2	Chan do Lindeiro	Mesolithic	9299 ± 167	41.34	99.51	Bos primigenius	P
Bp_CL_(9100BP)_3	Chan do Lindeiro	Mesolithic	9216 ± 185	8.81	68.19	Bos primigenius	P
Bt_CO_(60BP)	Chan do Lindeiro	Recent	100% modern C	7.24	45.99	Bos taurus	(T3)
Bt_AT_(3370BP)	A Tara	Bronze Age	3592 ± 101	17.78	94.83	Bos taurus	(T3)
Bt_AT_(1600BP)	A Tara	Suevic kingdom	1481 ± 87	8.98	63.25	Bos taurus	(T3)

Based on eight diagnostic SNP positions, the three B. taurus samples Bt_AR_(1200BP), Bt_AR_(2460BP), and Bt_SA_(1290BP)_2 could unambiguously be assigned to European haplogroup T3. The samples (Bt_SA_(1200BP)_1, Bt_AT_(3370BP), Bt_AT_(1600BP), and Bt_CO_(60BP)) displayed 3-6 SNPs characteristic for T3 but were missing data at the remaining diagnostic positions. The final sample Bt_MR_(3736BP) could not be assigned to a haplogroup with any confidence due to missing data at seven diagnostic positions. We did not detect evidence for any T1 (African descendant) haplogroups in our data. Haplogroup assignments based on the SNP analysis, results of dating, and further information about the samples are given in Table 1.

Tentative morphological identification of three of our samples as B. primigenius was confirmed by phylogenetic analysis: As shown in Fig 2, Bp_CL_(9100BP)_1, Bp_CL_(9100BP)_2, and Bp_CL_(9100BP)_3 grouped within the clade containing known B. primigenius mitochondrial genomes. The B. primigenius clade is a sister to a clade containing B. taurus haplogroups T and Q. This phylogeny also confirms all preliminary species assignments for our B. taurus bone samples that were based on morphology. In the phylogeny, our B. taurus samples are assigned to a large and poorly resolved haplogroup T-specific clade that contains all sub-haplogroups that were included in this analysis, T1-T5. The samples from Arcoia are found within the T3 subclade of the large T clade, while all other samples are found either basal to the T3/4 clade or even basal to a clade of T3/4 and T1. Overall, this T clade shows only high support values at the branches that seperate T5 and T2 from the remaining T sub-haplogroups. The clade containing all B. taurus haplogroup R sequences is sister to a large clade containing the haplogroups T, Q, and P.

Fig 2

Phylogenetic tree of mitochondrial genomes of different cattle species and haplogroups.

Result of a Bayesian analysis of mitochondrial genomes from five common cattle mitochondrial haplogroups (T,Q,P,R,I), which include the species Bos taurus (T,Q,R), Bos indicus (I), and Bos primigenius (P), and five subhaplogroups of T (T1,T2,T3,T4,T5). Reference mitochondrial genomes are coloured according to their haplogorup. Genomes generated in this study are denoted by their sample name. Squares on the nodes of the phylogeny mark a posterior probability larger than 0.7.

The statistical parsimony network that was computed from the B. taurus-specific alignment (Fig 3) shows a central hub representing mitochondrial genomes of different haplogroups. Most of these belong to T3, but all other haplogroups are present as well. Surrounding this central hub, nodes separate into haplogroup-specific satellite groups for T1, T3, and T. Where haplgroup T here refers to the sub-haplogroup T described in Troy et al., 2001 [3] and not to the general haplgroup. The samples Bt_AR_(2460BP), Bt_SA_(1290BP)_2, and Bt_AR_(1200BP) are deep within the T3 group. All other samples are closer to the central hub. Interestingly, our oldest B. taurus sample Bt_MR_(3736BP) was inferred as the central node of the central hub. However, the large amount of missing data in the sequence of this samples mitochondrial genome might influence its positioning in the network. The samples Bt_AT_(3370BP) and Bt_CO_(60BP) protrude into the sub-haplogroup T satellite group, and the remaining samples Bt_AT_(1600BP) and Bt_SA_(1200BP)_1 are in the outer ring of the central hub.

Fig 3

Statistical parsimony network of Bos taurus mitochondrial genomes.

Statistical parsimony network of 243 Bos taurus haplogroup T mitochondrial genomes, including the four sub-haplogroups T, T1, T2, and T3. Circle sizes are proportional to the number of identical mitochondrial genomes in the node, and numbers inside nodes are exact numbers of mitochondrial genomes. Marks on the edges represent the number of mutations by which mitochondrial genomes differ from each other. Nodes that represent samples of this study are described according to their names shown in Table 1 and are colored in the upper panel (A) of the figure. The lower panel (B) shows colors for publicly available data. Colors correspond to haplogroups as determined by SNP analysis. Nodes corresponding to sequence data of the present study are colored as T3, although haplogroup assignment based on SNP is not certain for all of them (see text and Table 1 for details).

Compared to the B. taurus data, the statistical parsimony network of the B. primigenius samples shown in Fig 4 appears rather different. Sequences are separated in the network by overall more mutations than the B. taurus sequences. Our samples are older than published data and differ from other B. primigenius by at least ten mutations. This level of sequence divergence is comparable to that found in published data. Our sample Bp_CL_(9100BP)_1 appears as the the central node of the network, while all other nodes are surrounding this sample. The two other samples of this study are directly connected to the central sample, but separated from it by 15 and 22 mutations.

Fig 4

Statistical parsimony network of Bos primigenius mitochondrial genomes.

Statistical parsimony network of seven B. primigenius mitochondrial genomes. Marks on the edges represent number of mutations by which mitochondrial genomes differ from each other. This figure includes three samples sequenced for this study, they are described by their names according to Table 1, and four published sequences, indicated by their NCBI accession number. Consistent with Fig 3, colors are used for the new samples of the current study in panel A and for publicly available data in panel B.

Discussion

Ancient DNA analyses of domestic animals provide not only information about the history of these lineages and species, but also about the people who kept these animals. Along with the domestication of animals, there was a change in human societies to settlement and agriculture [2]. Applying phylogenetic analysis to mitochondrial genomes recovered from eleven cattle subfossil bones, we were able to confirm our tentative species assignments based on morphology. Samples identified as B. taurus group together with other, mainly modern, B. taurus samples, while our three B. primigenius samples form a clade with available mitochondrial genomes of other B. primigenius. Furthermore, the 14C dating of our samples shows greater ages of the B. primigenius remains than for those of B. taurus, which is consistent with observations for other parts of Europe, where B. primigenius populations declined with the growth of human societies until their extinction in 1627 in Poland [43].

Results of our analyses allow insights into the time frame of Neolithization in Galicia. The samples that were identified as individuals of the wild aurochs B. primigenius were dated to around 9,000 BP (Table 1). Insight into the historic context of the age of these samples is given by comparison with the neighboring region Cantabria. The age of the B. primigenius samples dates to a period when the Mesolithic, the period before the Neolithic, just became common in Cantabria [25], and therefore to a time before the large transition to agricultural societies had occurred. Unfortunately, our data set contains a large temporal gap, as the earliest B. taurus samples of this study date to 4082 ± 63, which, in the Cantabrian regions, is at the end of the Megalithic, the latest time period of the Neolithic [25]. Thus, based on our data we can conclude that by about 4000 BP, domestic cattle had arrived in the north of Spain, while 9,000 years ago, aurochs were still roaming the region. Determining when exactly in the intervening 5,000 years cattle were introduced and aurochs went extinct will have to await further archaeological findings.

Our results also provide insights into cattle populations from historical times. The study area was quite populated during the Iron Age, through a network of small fortified settlements (castros) that were still active during the Roman conquest and settlement in the 1st century AD [44]. The Roman settlement in the Sierra was likely important for the empire as well as substantial, as evidenced by intensive exploitation of numerous gold mines in the area [45, 46]. Therefore, we would expect to find contributions from foreign livestock into the existing population. In many regions across Europe, the Roman occupation was followed by an increase in cattle size, a phenomenon that has been attributed to both the introduction of more robust animals and specific livestock management [47-49]. In Roman Galicia, an increase in the size of cattle has also been described [50]. If this increase in size was due to the introduction of foreign livestock, either the introduced animals were genetically indistinguishable from the local livestock, or only males were introduced, as we did not find any evidence of a genetic turnover of the mitochondrial lineages in the area. Further introductions of foreign livestock could have occurred during the period of migration (4th to 6th centuries), when successive Germanic peoples settled in ancient Roman Hispania: first the Suevi (409-585 AD), followed by the Visigoths, until the arrival of the Muslims at the beginning of the 8th century [51]. These peoples had been moving around the declining Roman Empire for more than a century. However, it is unlikely that they carried herds of cattle with them. Even less contribution is be expected from the Muslims conquerors of the Iberian Peninsula, as they never settled in Galicia [52]. Thus, it is not surprising that the cows of the Middle Ages studied here also do not reveal any genetic turnover compared to previous periods. However, those that are complete enough to be investigated with respect to their morphology (from A Tara cave) show a small size and dentition with strong mesowear, probably as a result of their life in the highlands [53]. Thus, it seems that it is the environmental conditions and the type of feeding that determined the size of historical Galician cattle, rather than any specific genetic lineage.

Maternal haplogroups of B. taurus provide information about cattle origin and migration, and they were historically defined based only a few diagnostic positions within the D-loop region of the mitochondrial genome [3, 6, 9, 36]. The network and Bayesian analysis of our study show that also whole mitochondrial genomes of this species group roughly according to those defined haplogroups (Figs 2 and 3). Our network analysis shows a very detailed picture of well resolved haplogroup clusters but also of an additional group of mitochondrial genomes in the center of the network that contains all different haplogroups and not only our ancient samples but also modern ones. Troy et al. describe the structure of haplogroups in this species as star-like, where the haplogroups T1-3 surround the haplogroup T, suggesting all other haplogroups coalesced from T [3]. Although the network analysis of this study shows a similar picture of haplogroups T1-3 surrounding a central group, that central group here contains a mixed of all haplogroups, with T forming a derived group from that. Our analysis, therefore, hints towards a more complicated origin and history of B. taurus.

For the samples sequenced in this study, the assigned haplogroups and the knowledge about their origin can be used to infer a possible descendance and migration of the cattle remains found in Galicia. In the SNP analysis, we detected only SNPs typical for sub-haplogroup T3 in our Bos taurus samples, which is common in Europe and Britain [3], although four of them have missing data in some of the diagnostic positions. A look at the analyses of their full mitochondrial genomes supports the hypothesis of their Europoean descent for the two samples from Arcoia and the sample Bt_SA_(1290BP)_2. In the phylogeny resulting from the Bayesian analysis they group clearly within in the T3 clade, which is also supported by the network analysis. In the our phylogeny, the sample Bt_SA_(1200BP)_1 is placed in the T3 clade as well, but the network analysis placed it in the mixed group from which the haplogroups descend in a star-like formation. A European origin of this sample is, therefore, still likely but not as clear as for the samples mentioned previously. A more complicated picture becomes evident for the remaining samples. In the Bayesian analysis they are found to be basal to the T3/4 clade, with one of them even basal to a T3/4 and T1 clade. However, the support of these branches is weak. In the network, the samples Bt_CO_(60BP) and Bt_AT_(3370BP) seem to be more closely related to genomes of the haplogroup T, which is the primarily domesticated haplogroup, and from which all others were derived, but which was was propagated into modern cattle [3]. This might be in accordance with the Bayesian analysis, in which this haplogroup is missing but must be assumed to be basal to the overall T clade. This is plausible for the older sample from A Tara, but not necessarily for the very young sample from Chan do Lindeiro. Its placement in the network, however, might be due to missing data in determining regions of the mitochondrial genome due to its low mapping coverage. Our sample Bt_MR_(3736BP) is displayed in the very center of the network. Carbon dating identified it as the oldest of our samples, and it therefore might be more closely related to the founder population than to modern samples. However, the lack of all diagnostic SNPs and the overall lower coverage of this sample make an inference of the haplotype impossible and question its placement in the network. Regions in its mitochondrial genome which we could not characterize in this study can have a large impact on the results of all our analyses for this sample.

Our cattle haplogroup assignments are in agreement with ancient DNA data of human remains from Spain, which showed that the transition from Mesolithic to Neolithic societies was mediated by demographic processes, involving admixture with pioneering farmers of Near Eastern origin [23, 24, 54, 55]. This genetic signal of Neolithic immigration in humans is consistent with the first appearance of cattle carrying the T3 haplotype, likely introduced by the incoming farmers. Noticeable, although recent studies have shown prehistoric gene flow between human population from African and southern/central Iberia around the Middle Neolithic/Bronze Age [24, 56], we found no evidence of African diversity in the domestic cattle analyzed in this study, as none of them was related to the African haplogroup T1.

This study presents the oldest available sequenced mitochondrial genomes of B. primigenius to date. In the network analysis, one of our samples that is around 9100 years old is placed in the center of the network, from which not only the other two samples of this study connect, but also published samples (Fig 4). In light of the age difference between our samples and the oldest of the publicly available samples (6700 BP), this seems plausible. The statistical parsimony network of B. primigenius shows that individual mitochondrial genomes separate into distinct nodes that are separated by many mutational steps of samples of very different ages but also samples of similar age. Even the number of mutations separating the Galician samples is not very different from numbers of mutations between other samples. This suggests an overall higher genetic variability of B. primigenius mitochondrial genomes, which is likely because the domestication of B. taurus can be seen as a bottleneck through which this species went, caused by the initial domestication process, breeding, and selection by humans. Such a process did not happen in the naturally occurring population of wild aurochs, which spread across Europe and other continents. However, in comparison to B. taurus, there are very few sequenced mitochondrial genomes available for B. primigenius. Therefore, further sequencing and research of this lineage is necessary in order to resolve the historical population structure of B. primigenius.

Conclusion

This study identified and characterized remains of 11 different Galician cattle roaming this region throughout different ages. Altogether, this adds to the picture of the history of different cattle species in Galicia, but, due to the close relationship of domesticated cattle with human, also to the history of human colonization and migration in that area. In accordance with human studies, we presented evidence that people with livestock migrated to that region from the European mainland. Over centuries and many other human migrations, domesticated cattle in this region seem to have remained of European descent. The analyses of full mitochondrial genomes of cattle also gave insight in to the evolution of these animals in general. We showed that haplogroups defined based on only a small part of the mitochondrial genome were confirmed in analysis results of the full sequences. Additionally, these analyses confirmed the species assignment of the so far oldest sequenced mitochondrial genome of the extinct wild aurochs.

Dating information.

Table including detailed information about the dating of the samples as well as about the sites the samples were found.

(XLSX)

Click here for additional data file.

Mapping information.

Table includes detailed information about the results of the sequencing and mapping experiments done with the samples of this study.

(XLSX)

Click here for additional data file.

Publicly available data.

Table contains NCBI accession numbers of all publicly available data sets used for the analyses displayed in Figs 2–4.

(XLSX)

Click here for additional data file.

Diagnostic SNPs.

Image shows the the D-loop region of the mitochondrial genomes for all samples aligned. Only the bases in the diagnostic positions are displayed, other bases are displayed as “-”.

(PNG)

Click here for additional data file.

29 Jan 2021

PONE-D-20-34180

Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence

PLOS ONE

Dear Dr. Hofreiter,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

All comments need to be addressed before re-submission.

Please submit your revised manuscript by Mar 15 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Peter F. Biehl, PhD

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2.In your manuscript, please provide additional information regarding the specimens used in your study. Ensure that you have reported specimen numbers and complete repository information, including museum name and geographic location.

If permits were required, please ensure that you have provided details for all permits that were obtained, including the full name of the issuing authority, and add the following statement:

'All necessary permits were obtained for the described study, which complied with all relevant regulations.'

If no permits were required, please include the following statement:

'No permits were required for the described study, which complied with all relevant regulations.'

For more information on PLOS ONE's requirements for paleontology and archaeology research, see https://journals.plos.org/plosone/s/submission-guidelines#loc-paleontology-and-archaeology-research.

3. We note that you are reporting an analysis of a microarray, next-generation sequencing, or deep sequencing data set. PLOS requires that authors comply with field-specific standards for preparation, recording, and deposition of data in repositories appropriate to their field. Please upload these data to a stable, public repository (such as ArrayExpress, Gene Expression Omnibus (GEO), DNA Data Bank of Japan (DDBJ), NCBI GenBank, NCBI Sequence Read Archive, or EMBL Nucleotide Sequence Database (ENA)). In your revised cover letter, please provide the relevant accession numbers that may be used to access these data. For a full list of recommended repositories, see http://journals.plos.org/plosone/s/data-availability#loc-omics or http://journals.plos.org/plosone/s/data-availability#loc-sequencing.

4. Please amend the manuscript submission data (via Edit Submission) to include author Stefanie Hartmann.

5.We note that Figure(s) 1 in your submission contain map images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

1) You may seek permission from the original copyright holder of Figure(s) 1 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

b)  If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

6. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.

Additional Editor Comments :

Please address the comments by reviewer 2 before re-submission.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The paper "Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence" provide new and useful data about cattle population dynamics in a non-sampled region of the Iberian peninsula. This data is important for the knowledge of the world cattle phylogeny and evolution as well as to understand ancient human history. The paper is very well written and developed and the analyses are appropriate and correct. I don't have any changes to highlight.

Reviewer #2: This work describes the retrieval of mtDNA genome sequence from 11 ancient Iberian Bos specimens. The data, as described, is of good quality; many of the genomes are of high coverage, the methods are sound, the laboratory track record is excellent. I think these data are of interest - three very old wild ox specimens are included for example, and the phylogenetic placing of some of the domestic sequences at basal phylogenetic positions is notable. Therefore my overall recommendation is for acceptance.

However I have a few comments, line numbers are given

20. the founding population estimate should be quoted as a range rather than the central estimate of 81 females.

table 1. the % coverage of three samples: Bt_MR, Bt_C, Bt_AT_1600BP and one Mesolithic are much lower than others. This needs to be borne in mind when one discusses the phylogenetic positions of these. For example is Bt_MR basal in position rather than simply incompletely described for several mutations that would remove it from this position?

174 relating to this point above, Bt_MR is inferred as the central node - but is it?

268, same point. Are we sure these are basal.

Fig 2, the BEAST tree. This analysis has the facility to consider nodes with different time horizons. A strong aspect of the data here is that these samples have direct C14 dates. I suggest that a better or additional analysis might only use the best covered mtDNA genomes and time-stamp the ancient samples. It would be interesting if the topology of the tree changed in this circumstance.

The discussion is somewhat lengthy and could be abridged. For example, I think that the relevance/comparison to the aDNA analysis of human remains is over-discussed. The only point (and it is worth making) is the absence of T1.

323 An archaeologist might query the claim that this presents “ the first evidence that people with livestock migrated into that region from the European mainland.” Is there no other non aDNA evidence for this?

Data availibility: It is welcome that all raw fastq data will be deposited. This is because processed data (eg mtDNA sequences) are subject to local parameter choices and subsequent users may wish to process these slightly differently.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Dan Bradley

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

19 Mar 2021

Dear Editor,

thank you and the reviewers for the comments to our manuscript entitled “Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence.” As requested, we have revised our manuscript and include a version with tracked changes. The point-by-point responses to the reviewers’ comments are below.

Editor Comments:

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

Response: We have used the most current version of the LaTeX template, available at https://journals.plos.org/plosone/s/latex. This template appears to be more recent than the listed PDFs.

2. In your manuscript, please provide additional information regarding the specimens used in your study. Ensure that you have reported specimen numbers and complete repository information, including museum name and geographic location. If permits were required, please ensure that you have provided details for all permits that were obtained, including the full name of the issuing authority, and add the following statement: 'All necessary permits were obtained for the described study, which complied with all relevant regulations.' If no permits were required, please include the following statement: 'No permits were required for the described study, which complied with all relevant regulations.'

Response: All specimens studied for our manuscript are housed in the "Instituto Universitario de Xeoloxía, Universidade da Coruña", the director of which is one of the co-authors (Aurora Maria Grandal D'Anglade). Therefore, no permits were required, and we have added the corresponding statement to the manuscript.

3. We note that you are reporting an analysis of a microarray, next-generation sequencing, or deep sequencing data set. PLOS requires that authors comply with field-specific standards for preparation, recording, and deposition of data in repositories appropriate to their field. Please upload these data to a stable, public repository (such as ArrayExpress, Gene Expression Omnibus (GEO), DNA Data Bank of Japan (DDBJ), NCBI GenBank, NCBI Sequence Read Archive, or EMBL Nucleotide Sequence Database (ENA)). In your revised cover letter, please provide the relevant accession numbers that may be used to access these data. For a full list of recommended repositories, see http://journals.plos.org/plosone/s/data-availability#loc-omics or http://journals.plos.org/plosone/s/data-availability#loc-sequencing.

Response: We deposited all raw sequence reads of this study in the NCBI Sequence Read Archive under BioProject accession number PRJNA705960. All generated mitochondrial genomes were deposited in NCBI GenBank. All accession numbers have now been included in the manuscript.

The samples Bt_MR_(3736BP) and Bt_CO_(60BP) cannot be made available at NCBI SRA (really – you mean GenBank, or? SRA should be fine), since the do not meet the GenBank submission criteria (less than 50 \\% ambiguous characters). These are available on request.

4. Please amend the manuscript submission data (via Edit Submission) to include author Stefanie Hartmann.

Response: The author Stefanie Hartmann is now included in the manuscript submission data.

5.We note that Figure(s) 1 in your submission contain map images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: \\url{http://journals.plos.org/plosone/s/licenses-and-copyright}.

Response: The map was created by coauthor Aurora Grandal-d'Anglade specifically for this study; it has not been used elsewhere. We have added this information to the legend of Figure 1 in our manuscript.

6. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: \\url{http://journals.plos.org/plosone/s/supporting-information}.

Response: We reorganized the data of our supporting information into several different files that follow the PLOS ONE supporting information guidelines. Captions and in-text citations in our manuscript were updated as requested.

Reviewer 2 comments:

20. the founding population estimate should be quoted as a range rather than the central estimate of 81 females.

Response: We addressed this comment in line 20 of our manuscript. The sentence now reads: "Based on mitochondrial DNA data, a founding population of 29 to 783 females has been estimated for the beginning of the taurine cattle domestication in the Middle East [4, 17]"

table 1. the % coverage of three samples: Bt_MR, Bt_C, Bt_AT_1600BP and one Mesolithic are much lower than others. This needs to be borne in mind when one discusses the phylogenetic positions of these. For example is Bt_MR basal in position rather than simply incompletely described for several mutations that would remove it from this position?

Response: We agree with this comment, especially for the samples Bt_MR and Bt_CO. For both of them, we only were able to retrieve sequence information that covered less then 50% of the mitochondrial genome. Therefore, we added additional explanations for their placements in the phylogeny and the network in the discussion. Specifically, the following short sections were added.

line 271: "This is plausible for the older sample from A Tara, but not necessarily for the very young sample from Chan do Lindeiro. Its placement in the network, however, might be due to missing data in diagnostic regions of the mitochondrial genome due to its low mapping coverage."

line 277: "However, the lack of all diagnostic SNPs and the overall lower coverage of this sample make an inference of the haplotype impossible and question its placement in the network. Regions in its mitochondrial genome which we could not characterize in this study can have a large impact on the results of all our analyses for this sample."

174 relating to this point above, Bt_MR is inferred as the central node - but is it?

Response: This point is related to previous one, but is referring to a line in the results. Because in this line we only wanted to describe the samples placement in the network, we did not make any changes there. However, we do agree with the overall point that the central position of this sample might might be due to missing data at potentially important sections of the mitochondrial genome.

268, same point. Are we sure these are basal.

Response: We addressed this comment specifically with the changes made in the discussion at line 271, shown two comments above.

Fig 2, the BEAST tree. This analysis has the facility to consider nodes with different time horizons. A strong aspect of the data here is that these samples have direct C14 dates. I suggest that a better or additional analysis might only use the best covered mtDNA genomes and time-stamp the ancient samples. It would be interesting if the topology of the tree changed in this circumstance.

Response: We thank the reviewer for this suggestion and tried to implement it. Unfortunately, the BEAST analysis did no coalesce even after a very long run time (> 1 billion steps). Therefore, we assume that our data set does not contain enough temporal information for this type of analysis. Consequently, we have omitted it from the revisions.

The discussion is somewhat lengthy and could be abridged. For example, I think that the relevance/comparison to the aDNA analysis of human remains is over-discussed. The only point (and it is worth making) is the absence of T1.

Response: In order to abridge the discussion, we shortened the paragraph about the comparison with human remains. The new paragraph starts in line 282 of the manuscript and now reads as follows:

"Our cattle haplogroup assignments are in agreement with ancient DNA data of human remains from Spain, which showed that the transition from Mesolithic to Neolithic societies was mediated by demographic processes, involving admixture with pioneering farmers of Near Eastern origin [23, 24, 44, 45]. This genetic signal of Neolithic immigration in humans is consistent with the first appearance of cattle carrying the T3 haplotype, likely introduced by the incoming farmers. Noticeable, although recent studies have shown prehistoric gene flow between human population from African and southern/central Iberia around the Middle Neolithic/Bronze Age [24, 46], we found no evidence of African diversity in the domestic cattle analyzed in this study, as none of them was related to the African haplogroup T1."

323 An archaeologist might query the claim that this presents “ the first evidence that people with livestock migrated into that region from the European mainland.” Is there no other non aDNA evidence for this?

Response: We agree with this comment, as well, and changed the concerned line 315 accordingly. It now reads: "In accordance with human studies, we presented evidence that people with livestock migrated to that region from the European mainland."

Data availibility: It is welcome that all raw fastq data will be deposited. This is because processed data (eg mtDNA sequences) are subject to local parameter choices and subsequent users may wish to process these slightly differently.

Response: The raw fastq data is deposited in the NCBI Sequence Read Archive under BioProject accession number PRJNA705960. Accession numbers have now been included in the manuscript.

Other changes:

While the manuscript was in review, we received new C14 dating information for three samples of this study. The concerned samples are Bt_AT_strat, Bt_SA_strat_1, and Bt_SA_strat_2. Their sample names were changed accordingly to Bt_AT_3370BP, Bt_SA_1200BP_1, and BT_SA_1290BP_2. Due to the new samples names, we made changes in Figure 2 and Figure 3 and also changed the names in the manuscript and Table 1.

Yours sincerely,

Marie Gurke, on behalf of all co-authors

22 Mar 2021

Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence

PONE-D-20-34180R1

Dear Dr. Hofreiter,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Peter F. Biehl, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

5 Apr 2021

PONE-D-20-34180R1

Insight into the introduction of domestic cattle and the process of Neolithization to the Spanish region Galicia by genetic evidence

Dear Dr. Hofreiter:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Peter F. Biehl

Academic Editor

PLOS ONE