Introduction
Archaeological faunal remains have been widely used to address various questions in biology and social sciences. The scope of these range from the demographic history of wild populations, which can inform about ecological dynamics and conservation biology, to animal management and breeding practices, providing insights into the subsistence strategies and lifeways of prehistoric human societies that exploited animals1–6. A key step here is the accurate taxonomic identification of animal remains. However, distinguishing morphologically similar species in zooarchaeological material is a prevailing challenge, constrained by the high level of similarity between skeletal elements, the fragmented state of excavated specimens (possibly with missing fragments), and the absence of morphological markers in subadults7,8. The need for an effective approach to identify species’ remains accurately has thus led to the development of several alternative methods, including isotope analyses, protein fingerprinting, and ancient DNA (aDNA) analyses9–14.
The majority of non-human aDNA data today is produced using shotgun DNA sequencing on Illumina platforms15. Beyond species identification, such data from well-preserved zooarchaeological samples can yield a wealth of information to study demographic and evolutionary history. However, relatively old (e.g. >1000 years old) zooarchaeological samples from regions with humid, temperate or warmer environments are mostly poorly preserved16. Cooking and other forms of heat treatment before human consumption may additionally degrade organic material17. In such poorly preserved samples, the proportion of endogenous DNA among the total DNA read pool will be low, down to 1% or even lower18. Accordingly, most experiments can produce only low amounts of DNA sequence data, if any, from zooarchaeological samples from temperate regions within reasonable budgets; such genomic data frequently remain at genome-wide depths of coverage <0.1x per sample19,20.
Theoretically, even 0.1x coverage genome data could allow accurate taxonomic identification by comparative alignment, i.e. mapping reads to the reference genomes of alternative candidate species, such as sheep versus goat, or horse versus donkey, and comparing coverages or mismatch rates. However, this only applies to situations where both species have assembled nuclear reference genomes (e.g. no such reference is available for the donkey). Even in cases where nuclear genomes from both species are available (e.g. sheep and goat), the limited amount of shotgun sequencing data available from poorly preserved samples, quality differences between the reference genomes, the highly fragmented nature of aDNA hence short read lengths, and postmortem damage can introduce high levels of uncertainty in the alignment process. The problem is further exacerbated when the sequence similarity between candidate reference genomes is high or when there exists strong differences in the genome assemblage qualities.
These call for new approaches for species identification with aDNA data. For instance, the Zonkey21 pipeline was developed for distinguishing horse, donkey and their hybrids by using nuclear aDNA variants with a clustering approach, but is only applicable for equid taxa from which there exists large datasets of genetic variation. Here we present a broadly applicable method, MTaxi, designed for distinguishing pairs of any closely related species using low amounts of shotgun aDNA sequencing data, whenever mitochondrial DNA (mtDNA) reference sequences are available. Our method focuses on mtDNA owing to its short size, haploid nature, having a lower rate of decay than nuclear DNA16, and having multi-copies per cell, which increases its availability relative to autosomes22, facilitating analyses. For example, across n=310 shotgun sequenced ancient DNA libraries from human, sheep, goat, horse and donkey generated by our group, each of which contained ≥0.01 endogenous DNA, the average ratio of mitochondrial DNA to nuclear DNA coverage was 87:1 (data not shown). The greater number of informative sites due to the high mitochondrial mutation rate is an additional advantage for taxon identification of closely related species23. Finally, the availability of mitochondrial reference sequences for a larger number of taxa (compared to a limited number of high quality reference genomes) allows our approach to be applied to a wider number of species, including extinct lineages. For instance, as of December 24 2021, the genome resources database from NCBI includes only 175 nuclear genomes for mammals24, compared to 1453 mitochondrial genomes25, an 8-fold difference. The MTaxi26 approach uses roughly the same information as competitive mapping in species assignment, but has the advantage that it calculates statistical significance for the result.
To exemplify the use of MTaxi we chose the case of sheep (Ovis) versus goat (Capra), two closely related species belonging to the same subfamily Caprinae. The aforementioned constraints on morphological identification causes a large proportion of sheep and goat remains to be only identified at the subfamily level as Caprinae8, and ambiguity which can significantly constrain zooarchaeological analyses, especially in the study of animal husbandry. Here we first estimate MTaxi’s accuracy using 1200 ancient mitogenome simulations with six different coverages from both sheep and goat. We then test its performance with n=9 ancient sheep samples19,20 and n=3 goat samples27. We further test MTaxi on the horse and donkey, a pair which also suffer from difficulties in morphological differentiation, using n=3 ancient horse, n=2 modern domestic horse, and n=5 modern domestic donkey samples adopted from the literature28–32.
Methods
Overview of the method
MTaxi makes use of the mismatch positions, i.e. putative substitutions, between two alternative candidate taxa, such as sheep and goat. We call these “target sites” and are obtained from pairwise alignment of mitochondrial reference genomes. Each read harbouring the target sites is classified according to the genotype compositions, and we identify the taxon using a binomial test for the read proportion of the sample (Figure 1).
Figure 1. Overview of the MTaxi pipeline.
Flowchart and representations of the steps to determine the sample taxon. Here sheep and goat stand for the candidate species pair, but MTaxi can be applied to any pair of species where mitochondrial DNA reference sequences are available. Target sites represent mismatches (candidate substitutions) between the reference genomes, restricted to transversions. Reads are assigned to either taxon based on target sites. Reads may be assigned to the wrong taxon due to homoplastic mutations, technical error, or incomplete lineage sorting.
Target sites
The method involves compiling a list of mtDNA target sites, representing likely substitutions between the species. To generate this list for sheep and goat, we first generated a pairwise alignment between sheep (Oar_v3.1) and goat (ARS1) mtDNA reference genomes via the R package Biostrings v.2.65.033 using default parameters, which yielded n=1699 single nucleotide substitutions. We then restricted these to transversions to avoid (a) confounding effects due to postmortem damage-induced transitions in ancient DNA, and (b) homoplasies that could arise by high-frequency transitions34. This yielded a set of n=197 transversion substitutions, which we refer to as the target site list 1.
We also created a subset of this, that we call target sites type 2, by removing polymorphisms in either species, which we reasoned might increase power by avoiding ambiguities. For this, we obtained a list of polymorphic sites using the software snp-sites v.2.4.135 from a data set assembled by Shi and colleagues36, which contains pairwise alignments (each sequence aligned to the reference genome) for mtDNA sequences belonging to n=47 domestic sheep and n=35 domestic goats. In this dataset, we identified n=57 and n=40 polymorphic single nucleotide positions overlapping with the n=197 target sites in domestic sheep and goat, respectively. After eliminating these polymorphisms we were left with n=120 positions, which we refer to as the target site list 2.
We applied the same procedures to horse and donkey by using mtDNA references NC_001640.1 and NC_001788.1. This resulted in n=1264 substitution sites, and restricting these to transversions yielded n=117 positions (Figure 2). The positions are concentrated around the D-loop, but are also represented across the mitochondrial genome following similar patterns between the two species.
Figure 2. Distribution of target sites along reference mitochondrial genomes.
The figure shows the position of target sites along (a) sheep and (b) goat (c) horse (d) donkey reference mitochondrial genomes. The sites represent transversion-type substitutions (n=197 for sheep and goat, and n=117 for horse and donkey). The figure was generated through CGview42.
Alignment and genotyping
Ancient DNA data processing. AdapterRemoval v.2.3.137 was used for trimming residual adapter sequences and merging paired-end sequencing reads with an overlap of at least 11 bp between the pairs. Whole genome sheep and goat data (FASTQ files) were aligned to both sheep (Oar_v3.1) and goat (ARS1) reference genomes using BWA aln v.0.7.1538 with parameters: “-n 0.01 -o 2” and by disabling the seed with “-l 16500”. Mitochondrial goat data (Gilat10, Shiqmim9, Kov27, Uiv17) were aligned to both sheep (NC_001941.1) and goat (NC_005044.2) mitochondrial references with the same parameters as well. To prevent the influence of the PCR duplicates, reads with identical start and end positions were removed using “FilterUniqSAMCons_cc.py”39. After the removal of PCR duplicates, reads with mapping quality scores (MAPQ) lower than 30 were filtered out using samtools (v.1.9)40. The reads mapping to the reference genome with >10% mismatches and having a length <35 bp were filtered out. Damage patterns which are characteristics of aDNA were estimated using PMDtools41 “--deamination” parameter. The reads aligned to mtDNA were extracted from whole genome alignments using samtools (v.1.9)40. For each alignment, we called the genotypes at our sites of interest using pysamv v.0.16.0.1, which also runs samtools (v.1.9)40; genotyping was performed with parameters “-B” and “-A”. As default, MTaxi uses both the reads that aligned only to one of the species’ references and the ones that aligned to both species’ references in the analysis (which we refer to as “all reads” below). Additionally, we included an option (“shared reads”), by which the reads are restricted to those that are aligned to both species’ references; this is a conservative approach that could eliminate the possible effects of quality differences between the two reference genomes. Using pybedtools v.0.8.143,44, we obtained the reads overlapping with the target sites. We note that aligning reads to both nuclear and mitochondrial genomes is superior to alignment only to the mitochondrial genome, because the latter can cause misalignment of nuclear mitochondrial DNA sequences (NUMTs) to the mitochondrial genome.
The alignment and genotyping procedures for whole genome ancient horse data were applied in the same way as described for sheep and goat data. However, for the alignment, due to lack of a nuclear reference genome for the donkey, equid reads were mapped only to mtDNA references of the two species.
For the comparison of whole genome mapping frequencies, the total number of bases aligned to both reference genomes were calculated using samtools stats40, and we calculated the number of mismatches for each alignment.
Modern DNA data processing. After removing residual adapter sequences using AdapterRemoval v.2.3.137, we mapped the whole genome data of modern horse and donkey at pair-ended mode to both horse (NC_001640.1) and donkey (NC_001788) mitochondrial reference genomes using BWA mem (version 0.7.15)38 module with the parameter ‘-M’, and sorted the output using samtools (v.1.9) sort40. Duplicates were removed using Picard MarkDuplicates. Reads with mapping quality scores lower than 20 were filtered out using samtools (v.1.9)40. Libraries from the same individual were merged using samtools (v.1.9) merge40, and then the same filtering and genotyping procedures described in Ancient DNA data processing section process were applied on all modern equid data.
Taxon assignment
Each read can carry either reference or alternative genotypes at its target sites. MTaxi uses this genotype data to assign reads to either taxon, species 1 (SP1) or species 2 (SP2). This is done by first calculating the frequency of alternative alleles per read. If an SP1 read was aligned to the SP1 genome, we expect no alternative alleles at target sites, and if aligned to the SP2 genome, we expect all alternative alleles. MTaxi retains reads with an alternative allele frequency of either 1 or 0, thus excluding reads with inconsistent genotypes (i.e. alternative and reference alleles mixed) at target sites. Such inconsistent variants could represent PCR or sequencing errors, convergent mutations, or incomplete lineage sorting. Having thereby assigned reads as SP1 or SP2, MTaxi uses the proportion of these two classes of reads to determine the sample taxon using a two-tailed binomial test with the null hypothesis of p=0.5.
Ancient mitogenome simulations
We simulated 1200 ancient sheep and goat mitochondrial genomes (100 sheep and 100 goats for each coverage) at six different coverages (0.5x, 1x, 2x, 3x, 4x, 5x) using gargammel45, and tested the accuracy of the method. The sequencing error was set to ~1% using the parameters “qs -10 qs2 -10”. The simulations for horse and donkey (100 horse and 100 donkey for each coverage) were run with the same parameters above, again at six different coverages (0.5x,1x,2x,3x,4x,5x). The same alignment and genotyping procedures described in Ancient DNA data processing section were applied to the simulated data, except that they were mapped only to the mtDNA references of the species.
Ancient and modern samples
We used published FASTQ files to study the performance of MTaxi on real data. For sheep, we used FASTQ files of n=5 ancient sheep individuals (TEP03, TEP62, TEP83, ULU26, ULU31) from Yurtman et al.19 downloaded from the European Nucleotide Archive (ENA) database (Table 1), and n=4 (OBI013,OBI014, OBI017, OBI018) ancientsheep individuals from Taylor et al.20 downloaded from ENA (Table 1). All data had been produced with Illumina sequencing using either whole genome shotgun sequencing or using SNP capture followed by sequencing. For goat, we used n=5 (Acem1, AP45, Azer3, Direkli1, Direkli6) ancient whole genome FASTQ files, and n=4 (Gilat10, Shiqmim9, Kov27, Uiv17) ancient mitochondrial capture FASTQ files, produced by shotgun Illumina sequencing and mtDNA capture-sequencing, and published by Daly et al.27, downloaded from ENA (Table 1). For equids, we used n=3 ancient28 and n=2 modern domestic horses29,30 and n=5 modern domestic donkey31,32 FASTQ files, downloaded from ENA (Table 1). We randomly downsampled the equid files to mtDNA coverages ranging from ~0.3x to ~4x using samtools view with the option “-s”40.
Table 1. Genome data used in the study.
The table lists ancient and modern-day genomes downloaded from European Nucleotide Archive (ENA), with the study accession IDs and sample aliases.
| Sample Aliases | Species | Age | Study Accession | Publication |
|---|
| TEP03, TEP62, TEP83, ULU26, ULU31 | Sheep | Ancient | PRJEB36540 | (19) |
| OBI013,OBI014, OBI017, OBI018 | Sheep | Ancient | PRJEB41594 | (20) |
Acem1, AP45, Azer3, Direkli1, Direkli6,
Gilat10, Shiqmim9, Kov27, Uiv17 | Goat | Ancient | PRJEB26011 | (27) |
| Au6, Et1, Ke14, Sp5 | Donkey | Modern | PRJNA431818 | (32) |
| Willy | Donkey | Modern | PRJEB24845 | (31) |
| FM1798 | Horse | Modern | PRJEB10098 | (30) |
| Twilight | Horse | Modern | PRJNA205517 | (29) |
| VHR031, VHR102, CdY2 | Horse | Ancient | PRJEB31613 | (28) |
Results
Application to simulated ancient mitogenomes
We first studied the performance of MTaxi using ancient-like mtDNA read data simulations. We produced n=1200 mtDNA read datasets at varying coverage, n=600 for sheep and n=600 for goat (Materials and Methods). Using n=197 transversion substitutions (target sites type 1), MTaxi assigned BAM files to their respective taxa with 100% precision (i.e. no false positives) across all mtDNA coverages from 0.5x–5x using the default (“all reads”) approach (Figure 3a,b). All simulated data had a recall (i.e. true positive rate) of 100% and no false positives, even at mtDNA coverages ≥0.5x (Figure 3a).
Figure 3. Results of the method applied to simulated ancient genomes at different coverages.
Binomial test p-values for comparing the proportions of reads assigned to (a–d) sheep versus goat, and (e,f) and horse versus donkey. For sheep and goat, results are based on transversion substitutions without (a,b) and with (c,d) the exclusion of polymorphic sites. Results for both pairs are obtained through default (a,c,e) and “shared reads” approaches (b,d,f). n refers to the number of simulated genomes in each case (100 for each species in a pair). The height of the blue bar represents the number of simulated goat/donkey and sheep/horse genomes identified correctly with p<0.05, and the height of the red bar represents the number of unidentified cases. No cases were misidentified. The height of the grey bar represents the trials that did not contain any reads aligned to target sites, and thus could not be evaluated.
We also tested the performance of two more conservative approaches. First, we tried the “shared reads” option, which uses only a subset of reads aligned to both genomes; here the recall was >50% at 0.5x, but reached >95% at 1x coverage (Figure 3b). This low recall appears to be caused by the lack of power to reject the null hypothesis due to the majority reads being eliminated by the “shared reads” approach.
Second, we repeated the analysis after eliminating polymorphic sites from the target substitution set (target sites type 2 with n=120 positions), and obtained 98% precision at 0.5x coverage (Figure 3c). Through the “shared” reads approach and eliminating polymorphic sites (target sites type 2), we could assign taxa with 99% precision only at 3x coverage (Figure 3d). This was again apparently caused by reduced power due to using fewer sites and fewer reads.
We also performed the same analysis for n=1200 datasets of horse or donkey. Using n=117 transversion substitutions (target sites type 1), we again achieved 100% precision and recall in taxonomic assignment (Figure 3e). Using the “shared reads” option, we obtained 97% recall at 0.5x, and 100% recall at >0.5x coverage (Figure 3f).
Overall, the simulations suggest that MTaxi can achieve full accuracy even at mtDNA coverages ≥0.5x. Conversely, limiting the analysis to subsets of reads aligned to both genomes or to non-polymorphic substitution positions reduces power, and does not increase accuracy.
Application to samples of known species identity
Sheep and goat. We tested MTaxi on n=9 published ancient sheep (Ovis aries / Ovis orientalis) samples with mtDNA coverages >0.2x, and n=9 published ancient goat (Capra hircus / Capra aegagrus) samples with mtDNA coverages >0.8x (Table 2). The samples were produced in three different laboratories and varied in their mtDNA coverage. MTaxi yielded 100% accuracy for all 18 samples using the default approach (“all reads”) (Table 2). The probability of correct assignment by chance across all 18 MTaxi-classified specimens would be only 0.0003%, indicating the overall accuracy of our method.
Table 2. MTaxi results on sheep/goat genome data of known species identity (default approach).
The analysis was performed with n=197 transversion substitutions between sheep and goat, (target sites type 1). “Taxon” stands for known identity based on full genome data of the same sample (Table 1); “mtDNA coverage” shows coverage when mapping reads to the mtDNA reference of the original species (e.g. mtDNA coverage using the sheep reference for sheep data) after the duplicates have been removed; “Total assigned reads” refers to reads that could be mapped to both mtDNAs and the ones that could map only to one of the species’ references with high quality, overlapped target sites, and could be assigned to either species; “Sheep reads” and “Goat reads” show the number of reads that could be unambiguously assigned to either species; “p-value” shows the two-sided binomial test p-value for the proportion of sheep and goat reads being equal, and “Identified taxon” shows the final taxon assignment.
| Samples | Taxon | mtDNA coverage | Total reads | Sheep reads | Goat reads | p value | Identified taxon |
|---|
| TEP03 | Sheep | 3.73327 | 285 | 279 | 6 | <0.001 | Sheep |
| TEP62 | Sheep | 45.4534 | 3181 | 2555 | 626 | <0.001 | Sheep |
| TEP83 | Sheep | 21.7152 | 1474 | 1236 | 238 | <0.001 | Sheep |
| OBI014 (OB20-06) | Sheep | 9.97996 | 1043 | 1037 | 6 | <0.001 | Sheep |
| OBI018 (OB20-04) | Sheep | 2.39606 | 235 | 232 | 3 | <0.001 | Sheep |
| OBI013 (OB20-01) | Sheep | 1.25211 | 110 | 108 | 2 | <0.001 | Sheep |
| ULU26 | Sheep | 1.39023 | 178 | 171 | 7 | <0.001 | Sheep |
| ULU31 | Sheep | 1.01956 | 112 | 107 | 5 | <0.001 | Sheep |
| OBI017 (OB21-06) | Sheep | 0.21726 | 29 | 29 | 0 | <0.001 | Sheep |
| Acem1 | Goat | 77.4073 | 7969 | 751 | 7218 | <0.001 | Goat |
| AP45 | Goat | 2.00415 | 217 | 0 | 217 | <0.001 | Goat |
| Azer3 | Goat | 49.7898 | 5351 | 144 | 5207 | <0.001 | Goat |
| Direkli1 | Goat | 58.1783 | 5579 | 488 | 5091 | <0.001 | Goat |
| Direkli6 | Goat | 106.625 | 10536 | 928 | 9608 | <0.001 | Goat |
| Gilat10 | Goat | 1.71545 | 171 | 1 | 170 | <0.001 | Goat |
| Shiqmim9 | Goat | 0.799651 | 84 | 1 | 83 | <0.001 | Goat |
| Kov27 | Goat | 2.99934 | 325 | 5 | 320 | <0.001 | Goat |
| Uiv17 | Goat | 0.909304 | 98 | 1 | 97 | <0.001 | Goat |
As observed in the simulations, using the "shared reads'' approach did not improve accuracy, and we could correctly assign only 15 samples, while 3 samples with the lowest coverage had too few reads for assignment at p<0.05 (Table S1, Extended Data). One sheep sample (ULU31), with mtDNA coverage at 1x, had no reads overlapping the target sites, and thus could not be analysed at all.
Interestingly, we observed 1-26% of reads misassigned with the default (“all reads”) approach. These could represent homoplasy, shared polymorphism, or PCR/sequencing error. However, they do not influence the final outcome.
Since the positions that remain polymorphic within species can introduce noise in downstream analyses, for sheep and goat, we also studied the performance of the method using target sites type 2 (excluding polymorphisms; n=120 sites). Using the default approach (“all reads”), this yielded 100% accuracy for all samples except one sheep, which had a coverage lower than 0.3x (Table S2, Extended Data). With the “shared reads” approach, 100% accuracy was achieved for only n=4 sheep and n=8 goat samples (Table S3, Extended Data). Meanwhile, one goat and three sheep samples did not have any reads aligned to sheep and goat references that contained the target sites. We also noted that species-misassigned reads identified in the samples were not eliminated by this procedure (Table 2, Table 3, Table S4 & S5 in Extended Data). This result resonates with the above result from simulations, that removing polymorphic sites lowers statistical power but does not improve accuracy, at least in the case of sheep/goat assignment.
Horse and donkey. Applying MTaxi on n=5 horse and n=5 donkey samples, our method yielded 100% accuracy with both approaches (Table 3, Table S4 in Extended Data). The overall rate of correct assignment in this sample set appears significant (one-sided binomial test p=0.001).
Table 3. MTaxi results on horse/donkey genome data of known species identity (default approach).
The analysis was performed with n=117 transversion substitutions between horse and donkey (target sites type 1). “Taxon” stands for known identity based on full genome data of the same sample (Table 1); “mtDNA coverage” shows coverage when mapping reads to the mtDNA reference of the original species (e.g. mtDNA coverage using the horse reference for horse data) after the duplicates have been removed; “Total assigned reads” refers to reads that could be mapped to both mtDNAs and the ones that could map only to one of the species’ references with high quality, overlapped target sites, and could be assigned to either species; “Horse reads” and “Donkey reads” show the number of reads that could be unambiguously assigned to either species; “p-value” shows the two-sided binomial test p-value for the proportion of horse and donkey reads being equal, and “Identified taxon” shows the final taxon assignment.
| Sample ID | Taxon | mtDNA
coverage | Total assigned
reads | Horse
reads | Donkey
reads | p value | Identified
taxon |
|---|
| Au6 | Donkey | 0.843071 | 54 | 1 | 53 | <0.001 | Donkey |
| Et1 | Donkey | 1.96281 | 123 | 1 | 122 | <0.001 | Donkey |
| Ke14 | Donkey | 4.03893 | 208 | 5 | 203 | <0.001 | Donkey |
| Sp5 | Donkey | 0.981224 | 57 | 2 | 55 | <0.001 | Donkey |
| Willy | Donkey | 0.490042 | 27 | 1 | 26 | <0.001 | Donkey |
| FM1798 | Horse | 1.83667 | 97 | 97 | 0 | <0.001 | Horse |
| Twilight | Horse | 0.608403 | 30 | 29 | 1 | <0.001 | Horse |
| VHR031 | Horse | 0.319088 | 47 | 45 | 2 | <0.001 | Horse |
| VHR102 | Horse | 0.72599 | 86 | 83 | 3 | <0.001 | Horse |
| CdY2 | Horse | 1.29874 | 174 | 168 | 6 | <0.001 | Horse |
Whole genome comparative alignment
Comparative alignment can theoretically be a simple alternative to MTaxi when nuclear reference genomes are available. Here we explored the performance of comparative alignment using sheep/goat assignment as a model.
First, we observed that among ancient sheep BAM files used in this study, mapping results revealed inconsistencies in terms of the total number of bases aligning to each reference genome (Figure 4). Out of 9 sheep datasets with known species identity, only 4 showed a higher number of bases aligning to the sheep reference relative to the goat reference. However, we did not observe a similar inconsistency for the ancient goat samples, all of which had a higher number of bases mapped to the goat reference genome, most likely due to higher assembly quality of the goat reference. Unsurprisingly, the number of bases aligned to the nuclear genomes may not be an appropriate statistic for taxon identification between closely related taxon pairs (see Discussion).
Figure 4. Whole Genome Comparative Alignment Results.
Total number of bases aligned to sheep versus goat reference genomes (left panels) and the mismatch proportions (right panels) for whole genome ancient sheep (upper panels) and goat (lower panels) samples.
We then analysed mismatch proportions of reads aligned to either genome. In all n=10 sheep and n=5 goat samples, including the lowest coverage samples with 0.0004x nuclear coverage, we found lower proportions of mismatches to their own reference genome. Again, this result is unlikely to happen by chance (one-sided binomial test p=3e-05). We note, however, that 3 of the sheep samples (TEP83, ULU26, ULU31) differ only marginally (by ~0.7%) in their sheep vs. goat mismatch proportions. This suggests that comparative alignment can be an alternative to mitochondrial analysis in species identification, although our results imply that its success may not be guaranteed in all circumstances.
Discussion
We showcase a simple but elaborate method to distinguish between closely related taxa using low coverage ancient DNA data, utilising mtDNA substitutions. Focusing on mtDNA is advantageous both in terms of high copy number and of greater variability (see Introduction). Also, the mtDNA/nuclear DNA ratio has been stated to correlate positively with a decrease in endogenous DNA content22, suggesting that it should be more likely to obtain higher amounts of mtDNA than nuclear DNA in particularly poorly preserved samples.
Simulation results showed that MTaxi can distinguish sheep vs. goat with full accuracy at mtDNA coverages ≥0.5x. We also obtained 100% correct results with 18 ancient samples of known identity, (one sided binomial tests p = 3e-06). Likewise, simulations of ancient horse and donkey data yielded 100% accurate results at mtDNA coverages ≥0.5x, while downsampled modern and ancient domestic equid samples (n=10) of known species identity were also assigned fully correctly (one-sided binomial test p=0.001). Overall, MTaxi appears as a simple and efficient tool for correct taxon identification using ultra-low coverage shotgun sequencing data.
Meanwhile, our results suggested that conservative modifications of the pipeline that involve limiting the analysis to “shared reads” or excluding polymorphic sites did not improve performance, but on the contrary reduced statistical power and recall.
We note that MTaxi is successful even at mitochondrial coverages of 0.5x, which is a level frequently reached in low-coverage sequencing experiments when there exists 1% endogenous DNA. For example, in the aforementioned shotgun sequencing dataset (see Introduction), we had n=226 ancient mammalian samples with 1-10% endogenous DNA (median = 3%), and each library sequenced to a size of up to 50 million total reads per sample (median = 415,146 reads); within this set 50% of the libraries reached mitochondrial DNA coverages ≥0.5x, sufficient for effective identification by MTaxi (while 62% and 84% of the libraries reached ≥0.3x and ≥0.1x, respectively).
Our observations on comparative alignment were also notable. The comparison of total number of bases mapped to sheep and goat reference genomes showed that mapping frequencies can be deceiving, even when analysing whole genome data. A sheep FASTQ file can align more widely to the goat reference genome, and the degree to which this occurs seems to vary among samples. Meanwhile, all the goat samples had higher numbers of bases mapped to the goat reference genome. The reason for the observed differences between the performance of goat and sheep samples in alignment of their respective genomes could be related to variability in reference genome qualities and/or polymorphism between the species (indeed, the N50 of the goat genome ARS1 is 26,244,591 while that for the sheep genome Oar_v3.1 is 40,376). More generally, this result indicates that taxon identification using only the number of aligned bases in comparative alignment is not reliable, although competitive mapping (aligning to the both genomes simultaneously) or using mismatch frequencies could be a more effective alternative.
Indeed, the comparison of mismatch proportions in comparative alignment appears to be a more robust approach based on our empirical sample of 15 sheep and goat samples, even at a nuclear coverage of 0.0004x (also a value typically displayed in low coverage sequencing experiments). This could be a simple solution for taxon identification if reference genomes are available for both taxa. Still, our observation that mismatch proportions can vary only marginally in some sheep samples mapped to goat (e.g. TEP83 and ULU26 in Figure 4), calls for caution in using this strategy. Moreover, MTaxi has the advantage over simple competitive mapping or mismatch frequency comparison approaches, such that it involves calculating statistical significance.
MTaxi would be expected to perform on any species pairs with a degree of divergence comparable to that of sheep and goat, and would be particularly convenient when reference nuclear genomes of one of the species is lacking, which precludes comparative alignment. Candidate taxa that pose challenges for zooarchaeological identification include several mammal species in families Cervidae (deer), Leporidae (rabbit/hare) and Bovidae (cattle/bison), and birds46–48. Horse and donkey are another such pair, on which we checked the performance of our method. Compared to the Zonkey pipeline21, designed to classify ancient equid samples, MTaxi does not require a reference panel and is solely based on mitochondrial DNA data, hence an easier and faster method of classification.
In summary, the performance of MTaxi will depend on various factors, including evolutionary divergence and reference genome qualities of the species pairs, but we expect it to be an effective tool in various settings, as long as mitochondrial introgression can be excluded. We also note that its parameters and the data processing steps can be fine tuned to adjust for particularities of the species in question, such as the exclusion of polymorphic sites.
Software availability
Source code available from: https://github.com/goztag/MTaxi
Archived source code at time of publication: https://doi.org/10.5281/zenodo.680630726.
Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).
Data availability
Underlying data
European Nucleotide Archive: Archaeogenetic analysis of Neolithic sheep from Anatolia suggests a complex demographic history since domestication. Accession number: PRJEB3654019.
European Nucleotide Archive: Evidence for early dispersal of domestic sheep into Central Asia. Accession number: PRJEB4159420.
European Nucleotide Archive: Ancient Goat Genomics. Study Accession: PRJEB2601127.
European Nucleotide Archive: Dezhou Donkey de novo sequencing and assembly. Accession number: PRJNA43181832.
European Nucleotide Archive: Improved de novo Genomic Assembly for the Domestic Donkey. Accession number: PRJEB2484531.
European Nucleotide Archive: Evolutionary genomics and conservation of the endangered Przewalski’s horse. Accession number: PRJEB1009830.
European Nucleotide Archive: Equus Genome sequencing and assembly. Accession number: PRJNA20551729.
European Nucleotide Archive: Tracking five millennia of horse management with extensive ancient genome time-series. Accession number: PRJEB3161328.
Extended data
Zenodo: Extended Data for MTaxi, https://doi.org/10.5281/zenodo.675905349.
This project contains the following extended data within the file “Supplementary_tables.pdf”:
- Table S1 : MTaxi results on sheep/goat genome data of known species identity using target sites type 1 using the “shared reads” approach
- Table S2 : MTaxi results on sheep/goat genome data of known species identity using target sites type 2 with the default approach (“all reads”)
- Table S3 : MTaxi results on sheep/goat genome data of known species identity using target sites type 2 with the “shared reads” approach
- Table S4 : MTaxi results on horse/donkey genome data of known species identity using target sites type 1 with the “shared reads” approach
Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).
