Contamination with exogenous DNA is a constant hazard to ancient DNA studies, since their validity greatly depend on the ancient origin of the retrieved sequences. Since contamination occurs sporadically, it is fundamental to show positive evidence for the authenticity of ancient DNA sequences even when preventive measures to avoid contamination are implemented. Recently the presence of wheat in the United Kingdom 8000 years before the present has been reported based on an analysis of sedimentary ancient DNA (Smith et al. 2015). Smith et al. did not present any positive evidence for the authenticity of their results due to the small number of sequencing reads that were confidently assigned to wheat. We developed a computational method that compares postmortem damage patterns of a test dataset with bona fide ancient and modern DNA. We applied this test to the putative wheat DNA and find that these reads are most likely not of ancient origin.https://doi.org/10.7554/eLife.10005.001
Ancient DNA, that is to say DNA extracted from fossils and ancient remains, provides a window into the past lives of humans, animals and plants. But working with ancient DNA is challenging; DNA decomposes with time, and so ancient DNA is often fragmented, damaged and present in tiny quantities. Furthermore, ancient DNA is also easily contaminated by modern DNA from those handling it and its surroundings. Researchers have therefore developed special protocols for working with ancient DNA and tests for its contamination.
One approach used to check that DNA is of ancient origin identifies a pattern of damage that is specific to ancient DNA. This damage changes the building blocks that make up DNA, causing one (called cytosine or C) to be misread as another (thymine or T). This substitution occurs most frequently at the ends of ancient DNA molecules, and occurs less often along its length, forming a detectable and characteristic pattern of damage.
A common way to analyse ancient DNA is to sequence it and then compare the resulting sequences to the genomes of modern organisms to identify its origins. In a study published earlier in 2015, investigators sequenced the DNA present in sediments obtained from a submerged archaeological site off the coast of the Isle of Wight in the United Kingdom. This previous study identified some DNA fragments that matched sequences in the wheat genome. This led the investigators to conclude that wheat was present in the British Isles around 8000 years ago, some 2000 years earlier than previously thought.
However, possibly owing to the small number of fragments that were found, the previous study did not check if the damage pattern matched that expected for ancient DNA. Now, Weiß et al. have developed a new computational method that tests whether DNA shows a typically ancient, or typically modern, pattern of C-to-T substitutions. When this test was used to assess the wheat sequences that were previously claimed to have ancient origins, it revealed that their pattern of DNA damage did not fit statistically with those of ancient DNA.
Weiß et al.'s findings contest those of the earlier study, and suggest that the new statistical method could be used to authenticate ancient DNA even when the number of available sequences is low.https://doi.org/10.7554/eLife.10005.002
The evolutionary reconstruction of the past has been greatly enriched by direct interrogation of ancient DNA (aDNA) from plants and animal remains (Shapiro and Hofreiter, 2014). Although a vast proportion of flora and fauna do not fossilize, traces of their DNA may be preserved in sediments allowing the characterization of past biodiversity (Pedersen et al., 2015). A challenge to exploiting such resources is the ubiquitous threat of contamination with exogenous DNA. Therefore, special sample preparation procedures have been developed to reduce DNA contamination (Cooper and Poinar, 2000). Nevertheless, it remains difficult to estimate how well preventive measures work. If contamination is a possible explanation for the result, it is crucial to exclude this possibility by giving positive evidence for the authenticity of aDNA (Prüfer and Meyer, 2015). Fortunately, a large number of full-length DNA sequences can be generated using next generation sequencing, which allows for the authentication of aDNA. In aDNA an excess of C-to-T (cytosine to thymine) substitutions occur at the 5′ and 3′ ends of molecules (or its mirror image G-to-A (guanine to adenine) at the 3′ end, depending on the library protocol employed). When considering the 5′ end of sequences, the excess of C-to-T substitutions is highest at the first base and decreases exponentially towards the center (Figure 1A). This pattern is the result of cytosine deamination to uracil in single stranded overhangs (Briggs et al., 2007). Since it is present in aDNA-derived sequences but absent in much younger samples, it has been used as an authentication criterion in aDNA experiments (Krause et al., 2010; Prüfer and Meyer, 2015).
Smith et al. analyzed sediments from Bouldnor Cliff, a submerged archeological site in the United Kingdom, and suggested the presence of domesticated wheat 8000 years ago based on sedimentary ancient DNA (sedaDNA). This is 2000 years earlier than expected based on archeological remains in the British Isles and 400 years earlier than in nearby European sites (reviewed in Smith et al.). Since Smith et al. did not find wheat pollen or archeological remains associated with wheat cultivation, they conclude that the wheat presence in Bouldnor Cliff was the result of trading.
In total they produced ∼72 million Illumina reads, of which they robustly assigned 152 to wheat (Triticum), with dozens more (160 reads) to higher taxonomic ranks that include wheat. Smith et al. took state-of-the art preventive measures to avoid contamination and exercised great effort to ensure the accuracy and robustness of their phylogenetic assignments. The authors attempted to authenticate the aDNA molecules based on the expected excess of C-to-T substitutions, but because of the very small number of reads assigned to wheat, they failed to do so using standard approaches. As a result of that, the authors did not present any positive evidence for the ancient origin of their reads. Here we present an approach that compares the pattern of C-to-T substitutions in a set of test reads with the distributions of C-to-T substitutions in reads from known ancient- and modern-DNA and apply this approach to sedaDNA from Smith et al.
Although the excess of C-to-T substitutions at the 5′ end occurs at different magnitudes in samples of different ages, the exponential increase of substitutions towards the end is a ubiquitous pattern in aDNA studies (Sawyer et al., 2012). In order to score the presence of this pattern in various datasets, we fitted an exponential function and evaluated the goodness of fit by using a one-sided t-test to test for significant exponential decay. As expected, true aDNA libraries show significant goodness-of-fit p-values (Figure 1A), whereas non-significant goodness-of-fit p-values, neither decay nor growth, are observed in libraries derived from modern DNA (Figure 1B). A given C-to-T damage pattern plot can thus be summarized by its goodness-of-fit p-value that when it is significant indicates C-to-T exponential decay at the 5′ end (Figure 1A).
We resampled (with replacement) 10,000 sets of 150 sequences from a library of historic Solanum tuberosum collected in 1846 (Yoshida et al., 2013a). The number was selected to be comparable to the 152 reads that Smith et al. assigned to wheat. An empirical distribution of goodness-of-fit p-values was generated by performing the goodness-of-fit test for each subsample (Figure 2A). When we evaluate the sedaDNA goodness-of-fit p-value, we find that it falls within the upper 3% of subsamples with the least good fit. We can therefore reject the null hypothesis that the sequences assigned to wheat are as ancient as the historic S. tuberosum library. We repeated the whole procedure using this time a modern wheat library to generate the distribution of goodness-of-fit p-values (Figure 2A) and find a better match (p = 0.83). Thus, we cannot reject the hypothesis that the sequences assigned to wheat are of modern origin.
We sought to investigate how the test behaves when the empirical distribution of goodness-of-fit p-values is generated from different aDNA libraries. For this purpose we used a set of samples from animal (Sawyer et al., 2012) and plant remains (Yoshida et al., 2013) with an age of 85–170 years before present, and scored the sedaDNA wheat sequences against distributions generated from these libraries (subsamples of 150 sequences again). We observed that the goodness-of-fit p-value for the libraries is positively correlated with the empirical p-value for the sedaDNA wheat sequences tested against them (Figure 2B). Using a significance level of 0.05, we rejected the hypothesis that the wheat sequences are of ancient origin with 7 out of 13 libraries used in our test (Figure 2B). Thus, the purportedly 8000-year old wheat sequences show a less pronounced deamination pattern than many plant and animal samples with an age of less of 200 years. Finally, we took a less conservative approach and scored the sedaDNA against a distribution of goodness-of-fit p-values (subsamples of 150 read) generated from a 7000-years-old human Mesolithic sample from la Braña site in Northern Iberia (Olalde et al., 2014). La Braña is a site with cold environment and stable thermal conditions that has yielded exceptionally well conserved human fossils with ∼50% of human endogenous DNA that reach a ∼15% C-to-T substitution rate at the 5′ end (Olalde et al., 2014a) (Figure 2—figure supplement 1). We could reject the null hypothesis that the sedaDNA reads are as ancient as the sample from la Braña (p = 0.0014), a sample that is closer in time with the allegedly 8000-year-old wheat reads (Figure 2—figure supplement 2). It is worth pointing out that almost all 10,000 subsamples from la Braña had a very low (close to 0) goodness-of-fit p-value, even though we subsample only 150 reads (Figure 2—figure supplement 2).
We assessed the statistical power of the test by testing both an aDNA (Figure 3A) and a modern DNA library (Figure 3B) against a distribution built from a bona fide aDNA library, while varying the number of sampled sequences. Whereas the hypothesis that a true aDNA library is ancient was never rejected (Figure 3A), the hypothesis that a modern library has ancient origin could be rejected only when sufficient number of sequences were used for the subsample test (in tests with more than 300 reads the median empirical p-value was always below 0.05) (Figure 3B).
Finally, we skipped the phylogenetic curation step applied by Smith et al. to reduce the number of false positive wheat alignments, and mapped all reads sequenced by Smith et al. to the wheat genome. After stringent filtering of sedaDNA mappings we repeated our test varying the size of the subsample sets from 100 to 1000 reads. The empirical p-value was dependent on the number of reads tested, and declined with an increasing number of tested reads for all layers of sediments sequenced in Smith et al (Figure 3C). This pattern resembled the one obtained from a modern DNA library (Figure 3B). As for the phylogenetic curated 152 sequences, we were able to reject the hypothesis that the mapped reads are of ancient origin (mean p-value < 0.05 for all tests with more than 400 reads for layers 1–2 and 4, and 800 reads for layer 3). Our analysis also shows that the 152 sequences after phylogenetic curation are not a biased subsample from the distribution of all wheat-matching sequences.
We were able to reject the hypothesis that the sequences assigned to wheat by Smith et al. are of ancient origin. This is true even when we compared the putative 8000 year old sequences with only century old samples that show much lower deamination signatures. This means that a scenario in which wheat was transported to the Bouldnor Cliff site 8000 years ago is unwarranted. Our approach for authentication of aDNA can be used even with a very small number of sequences, and we hope that it will proof useful to test for positive evidence of authenticity for ancient DNA studies whose conclusions rely heavily on the ancient origin of the analyzed sequences.
Reads from most of the samples were downloaded from the European Nucleotide Archive (Table 1 and Supplementary file 1), with the exception of the Gorilla gorilla reads that were provided directly by the authors (Sawyer et al., 2012). Adapters were trimmed for both paired- and single-end runs using the program Skewer (version 0.1.120) using default parameters (Jiang et al., 2014). For paired-end runs (Supplementary file 1) forward and reverse reads were merged requiring a minimum overlap of 10 base pairs (bp) using the program Flash (version 1.2.11) (Magoc and Salzberg, 2011). Merged or single-end reads were mapped as single-end reads against their respective nuclear or organellar genomes: S. tuberosum nuclear genome (Potato Genome Sequencing Consortium et al., 2011), Solanum lycopersicum nuclear genome (The Tomato Genome Consortium, 2012), Triticum aestivum nuclear genome (International Wheat Genome Sequencing C, 2014), G. gorilla mitochondrial genome (Xu and Arnason, 1996), Homo sapiens nuclear genome (Genome Reference Consortium Human Build 37). The mapping was carried out using BWA-MEM (version 0.7.10) with default parameters, which include a minimum read length of 30 bp (Li, 2013). PCR duplicates were removed after mapping using bam-rmdup (available at https://github.com/udo-stenzel/biohazard), which computes a consensus sequence for each cluster of duplicated sequences. Alignments were stored in the bam format (Li et al., 2009).
We used two different approaches to process the reads from sedimentary DNA (Smith et al., 2015).
Phylogenetic curated reads: we used a set of 152 reads assigned to tribe Triticeae and to genus Triticum by Smith et al. after phylogenetic curation. However, we consider the complete sequence and do not exclude the initial 10 nucleotides as was done in the original processing (Smith et al., 2015). Reads were then aligned to the wheat genome as described above.
All sedimentary DNA reads: we aligned independently reads from all four layers sequenced by Smith et al. to the T. aestivum nuclear genome (International Wheat Genome Sequencing C, 2014). Duplicates were removed and only alignments with mapping quality greater or equal than 30 were used for further analysis. Additionally, we include a sequence complexity filter based on entropy, which removed low complexity reads with entropy less or equal to 50. The entropy filtering was carried out with prinseq-lite (version 0.20.4) (Schmieder and Edwards, 2011).
For each set of aligned reads (complete libraries or subsamples) the C-to-T substitutions patterns along the 5′ end of the read were assessed using the program PMDtools (Skoglund et al., 2014). We fitted an exponential function to the frequency of C-to-T substitutions for the first 20 nucleotides at the 5′ end. The fitting was performed in R (http://www.r-project.org) using the nls function, which determines the nonlinear least squares estimates of the parameters in a nonlinear model. The fitting was carried out with the model formula: y ∼ N∗exp(−rate∗x). From the nls fitting we obtained the t-value and degrees of freedom for the rate parameter and then calculated a goodness-of-fit p-value by using a one-sided t-test.
Subsets of different alignment numbers were randomly sampled (with replacement) 10,000 times from alignments stored in the bam format (Li et al., 2009). The random sampling was performed using samtools view (Li et al., 2009). For every subset of alignments we assessed the fraction of C-to-T substitutions, fitted an exponential function and calculated a goodness-of-fit p-value as explained above.
Phylogenetic curated reads: we compare the goodness-of-fit p-value of our test set of 152 sedimentary DNA reads with distributions of goodness-of-fit p-values generated from bona fide modern and ancient DNA. For the distribution of goodness-of-fit p-values from aDNA, we count how many of them are equal or greater than the sedimentary DNA goodness-of-fit p-value. To calculate the empirical p-value of the test we subsequently divided this number by the total number of values in the empirical distribution. With this approach we test the null hypothesis that the test set of reads contains a signal of ancient DNA damage that is comparable or even more pronounced than the signal in the aDNA library used to generate the empirical distribution of goodness-of-fit p-values.
For the distribution of goodness-of-fit p-values from modern DNA, we count how many of them were smaller or equal than the sedimentary DNA goodness-of-fit p-value. We calculate the empirical p-value of the test by dividing this number by the total number of p-values in the empirical distributions. With this approach we test the null hypothesis that the test set of reads matches the absence of ancient DNA damage patterns seen in reads of modern origin.
All sedimentary DNA reads: We tested independently alignments from each of the layers sequenced by Smith et al. using a bona fide aDNA sample for the generation of the distribution of goodness-of-fit p-values. For each layer we tested 10 sets of different numbers of reads (from 100 to 1000 reads, with increments of 100 reads). For each layer and for each number of reads in the test set we repeated the test and calculated the empirical p-value 1000 times as described above.
Other ancient DNA and modern DNA libraries were tested using the same procedure.
Patterns of damage in genomic DNA sequences from a NeandertalProceedings of the National Academy of Sciences of USA 104:14616–14621.https://doi.org/10.1073/pnas.0704665104
Triticum aestivum strain: SynOpDH Genome sequencingTriticum aestivum strain: SynOpDH Genome sequencing, EBI European Nucleotide Archive, PRJNA250383. http://www.ebi.ac.uk/ena/data/view/PRJNA250383.
Aligning sequence reads, clone sequences and assembly contigs with BWA-MEMarXiv preprint arXiv:13033997.
Illumina HiSeq 2000 sequencing: Whole genome sequencing of the La Brana 1 sample, AmpliTaqGold libraryIllumina HiSeq 2000 sequencing: Whole genome sequencing of the La Brana 1 sample, AmpliTaqGold library, EBI European Nucleotide Archive, SRR1045127. http://www.ebi.ac.uk/ena/data/view/SRR1045127.
Quality control and preprocessing of metagenomic datasetsBioinformatics 27:863–864.https://doi.org/10.1093/bioinformatics/btr026
Separating endogenous ancient DNA from modern day contamination in a Siberian NeandertalProceedings of the National Academy of Sciences of USA 111:2229–2234.https://doi.org/10.1073/pnas.1318934111
Dates of Shotgun metagenomic study of sedimentary ancient DNA (sedaDNA), from four strata of sediment core taken from an Bouldner Cliff, a submarine archaeological site in the SolentDates of Shotgun metagenomic study of sedimentary ancient DNA (sedaDNA), from four strata of sediment core taken from an Bouldner Cliff, a submarine archaeological site in the Solent, EBI European Nucleotide Archive, PRJEB6766. http://www.ebi.ac.uk/ena/data/view/PRJEB6766.
A complete sequence of the mitochondrial genome of the western lowland gorillaMolecular Biology and Evolution 13:691–698.https://doi.org/10.1093/oxfordjournals.molbev.a025630
Resequencing Solanaceae (Potato and Tomato) 19th century samplesResequencing Solanaceae (Potato and Tomato) 19th century samples, EBI European Nucleotide Archive, PRJEB1877. http://www.ebi.ac.uk/ena/data/view/PRJEB1877.
Joseph K PickrellReviewing Editor; New York Genome Center & Columbia University,, United States
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for submitting your work entitled “Assessing ancient DNA authenticity with low-coverage data: a case study of wheat in the British Isles 8,000 years ago” for peer review at eLife. Your submission has been favorably evaluated by Mark McCarthy (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The authors use an approach based on testing whether C-to-T damage patterns seen in degraded DNA fits a model of exponential decay from sequence termini in order to test a recent claim of 8,000 year old wheat DNA at a submerged site in Britain (Smith et al. 2015, Science).
The reviewers identified two additional analyses that they considered essential:
1) The authors should show visually the damage patterns in the reads from Smith et al.; even if this plot is noisy, it will provide some visual sense of the data.
2) The authors should provide some analysis of whether their results are robust to the choice of aligner. BWA-MEM uses extensive soft-clipping, which distorts patterns of damage. The authors should consider BWA-ALN/samse as an alternative.
Other reviewer comments below are provided as suggestions but need not be addressed in a revision.
In this paper, the authors critically assess a claim from Smith et al. that wheat was present in the UK 8,000 years before present. Specifically, they argue that the sequencing reads claimed to come from ancient wheat samples in this study are instead most likely to be modern contaminants based on patterns of DNA damage.
1) The Smith et al. claim is based on 152 sequencing reads, and so the authors here have little to work with. However, they claim that the expected exponential increase in C->T substitutions in these reads is not present. It would have been nice to see the patterns of DNA degradation of these reads visually, rather than it being summarized by the p-value from their test. Why do the authors not show a figure like Figure 1A, except using the 152 reads from Smith et al? This would be extremely useful for visually determining if the expected damage patterns are present (even if the plot is extremely noisy due to small numbers of reads). The authors could show the pattern from the empirical data superimposed on damage patterns in their subsampled ancient and modern libraries.
Weiß et al. develop a pipeline to assess the authenticity of results of analyses of ancient sedimentary DNA, and use this pipeline to test recently published results. Such a pipeline is a useful tool for ancient sedimentary analyses, and is likely to be widely adopted in the field. I found the paper to be well written and straightforward. I have only one major query.
The authors choose to use BWA-MEM for alignment, which may not be the most appropriate aligner for this test and is likely to affect the results. This aligner (BWA-MEM) performs a lot of soft-clipping (unless the settings were modified from default, in which case the authors should state this), which causes two major issues. First, it makes misincorporation profiles unreliable. Second, it can cause a lot of random mapping. I strongly recommend the authors re-run these analyses using a more appropriate aligner, e.g. BWA-ALN/samse.https://doi.org/10.7554/eLife.10005.010
The reviewers identified two additional analyses that they considered essential:
1) The authors should show visually the damage patterns in the reads from Smith et al.; even if this plot is noisy, it will provide some visual sense of the data.
We thank the reviewers for the suggestion and agree that visualizing the damage patterns of the 152 reads assigned to Triticum by Smith et al. will facilitate the understanding of our analysis. We have modified Figure 1 and split it in two figures. In the new Figure 1 we have kept panels A and B, whereas in new Figure 2 we have original panels C and D. The new panel 2A now includes damage patterns from Smith et al. and also from different parts of the distributions of goodness-of-fit p-values of known ancient and modern DNA. We believe that the new figures provide an intuitive visual counterpart to different goodness-of-fit p-values.
2) The authors should provide some analysis of whether their results are robust to the choice of aligner. BWA-MEM uses extensive soft-clipping, which distorts patterns of damage. The authors should consider bwa aln/samse as an alternative.
We appreciate this very helpful comment and have performed a series of analyses to address it.
First we evaluate the amount of soft-clipping included in our BWA-MEM mapped reads. The percentage of mapped reads that are soft-clipped was very small for the reads assigned to Triticum in Smith et al. (5.4%), for the potato library reads used as known ancient DNA (aDNA) library (5.9%), and for the wheat library reads used as known modern DNA library (4%). Since our known aDNA library is affected by soft-clipping to a similar degree as the Smith et al. dataset, we do not expect a bias against an ancient signal. Nevertheless, we have repeated all analyses in an alternative approach using the software “PMDtools”, which handles soft-clipping by disregarding soft-clipped bases when estimating C-to-T damage patterns. Results were consistent using the two methods; in the manuscript we now report results based on “PMDtools”.
Second, we repeated the analysis by aligning all datasets with BWA-ALN using sensitive parameters that have been previously used for aDNA (no seed (–l 16500), allow more mismatches (–n 0.01) and allow up to two gap open events (–o 2)). Whereas we could map all 152 reads assigned to Triticum by Smith et al. to the Triticum aestivum reference genome using BWA-MEM, we could only map 78% of these reads (108 reads) using BWA-ALN (in Smith et al. reads were aligned using BLAST). We also aligned only a fraction of reads with BWA-ALN from the ancient potato (96%) and modern wheat libraries (73%). Despite this bias we found that the distribution of goodness-of-fit p-values from the known potato aDNA library with BWA-ALN resembled the distribution created with BWA-MEM (Author response image 1A). However, the wheat test dataset is not significantly different from the ancient potato reads (P=0.08) (Author response image 1B). This is most likely due to the reduced number of aligned reads; if we subsample 108 reads from the 152 reads in the original wheat dataset (10,000 mapping with BWA-MEM) and repeated the test the result is not significant in ∼48% of the cases.
Third, the result of the test also depends on the known aDNA library used to generate the empirical distribution of goodness-of-fit p-values. We took a rather conservative approach by using a S. tuberosum library that is only 169 year old, while the putative wheat reads from Smith et al. are claimed to be 8000 year old. Consequently, we have added a new analysis, where we used as known aDNA library from a 7.000-year-old Mesolithic human from La Braña site in Northern Iberia. We could again reject the null hypothesis that the test set of wheat reads are as ancient as the human fossil from La Braña (P=0.0014) (Figure 2–figure supplement 2). The distribution of goodness-of-fit p-values from La Braña shows that C-to-T signal is present in almost all subsamples of 150 alignments. The empirical p-values were significant independent on the choice of the aligner (Author response image 2).https://doi.org/10.7554/eLife.10005.011
- Michael Dannemann
- Kay Prüfer
- Michael Dannemann
- Clemens L Weiß
- Hernán A Burbano
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
We thank Moisés Expósito-Alonso, Janet Kelso, Rebecca Schwab and Detlef Weigel for discussions on data analysis and comments on the manuscript.
- Joseph K Pickrell, New York Genome Center & Columbia University,, United States
- Received: July 10, 2015
- Accepted: September 29, 2015
- Version of Record published: November 3, 2015 (version 1)
© 2015, Weiß et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.