BACKGROUND
the field of ancient dna  has reached 30 years of age, in which time it has progressed from the amplification of fragments of single loci to the complete sequencing of individual complex genomes  <cit> .
in its first two decades, adna research was primarily focused on pcr-based amplification and subsequent sanger sequencing of selected loci and organellar genomes, with results applied to analyses of population differentiation and phylogeography, phylogenetics, and even metagenomics .
in the last decade the field of adna has moved from the genetic to the genomic level with the advent of massively parallel sequencing platforms.
this has been accompanied by a concurrent shift in focus to full genome sequencing and assembly, and genome-scale analyses of population trends  <cit> .
as this field enters its fourth decade, it is inevitable that application of genome-level sequencing will become more commonplace.
this will facilitate the broader development of ancient population genomics and phylogenomics approaches  <cit>  in a wider range of taxa.
nonetheless, in order to fully benefit from genome-scale approaches, it is critical to understand how the unique characteristics of adna might impact results from high-throughput sequencing and read-mapping technologies.
in general, there are several important characteristics that distinguish ancient from modern dna samples.
first, not only are adna molecules highly fragmented  <cit> , but the fragmentation process itself is biased toward breakpoints bordered by 5′ purine residues   <cit>  which could result in biased coverage patterns.
further, deamination of cytosine residues at or near fragment ends ultimately leads to c-t and g-a misincorporations at 5′ and 3′ sequence ends, respectively  <cit> , and could confound adna read mapping and variant calling.
last, but not least, ancient samples typically feature relatively low proportions of endogenous dna  <cit>  such that sequencing multiple samples to higher coverage depths  to enable confident genomic assembly and variant calling, is often either cost-prohibitive or simply not possible due to limited sample or dna availability.
as a consequence of this, most reports applying high throughput sequencing to full genomes of ancient samples have been limited to one or several samples sequenced and at low or ultra-low  to moderate coverage depths  <cit> .
this combination can severely impact the accuracy of variant calls  <cit> .
while it is relatively simple to document the biases associated with mapping adna sequence reads , it is more difficult to measure the impacts of these biases on variant calling when analysing empirical data.
considering this, we designed a study to quantify the impact of adna fragmentation and misincorporation biases on variant calling using simulated data.
by utilizing simulated datasets, we were able to control for a number of potentially confounding factors that commonly feature in and could have an effect upon adna studies.
these include level of divergence between sample and reference, average read length, average sequencing depth and damage level .
in addition, by synthesizing in silico reference and sample sequences prior to incorporation of adna-associated damage patterns, it is possible to know the exact “pre-damage” nucleotide sequences of samples.
this allows direct comparison of called variants to the true variants, in terms of both accuracy and overall completeness.
RESULTS
to explore the impacts of post-mortem dna damage on variant calling, we first simulated triplicate pairs of reference and diploid sample sequences at each of three gc contents  and two divergence levels , and at ca.
<dig> mbp  lengths, as described in the methods.
short read sequence pools were artificially generated from the sample sequences at average read lengths of  <dig>   <dig> and 80 bp, and were further subjected to either zero, low or high levels of synthetic post-mortem damage reflective of previous reports.
this damage included both fragmentation bias and elevated 5′ c-t/3′ g-a transitions near sequence ends.
these pools were subsequently mapped back to their corresponding reference sequences using commonly applied alignment software and at a range of ultra-low to moderate sequencing coverage depths .
this method allowed us to directly compare known damage patterns and variant positions to those measured from mapped read pools, and ultimately to estimate whether and to what extent a variety of variables, including coverage depth, read length and damage level, might impact upon variant calling.
synthetic damage patterns
implemented patterns of fragmentation bias were clearly seen when mapping reads with both low and high levels of synthetic damage .
in low damage read pools, 5′ guanosine and adenosine residues were elevated by ca.
20% and 10%, respectively, over genomic levels at positions immediately 5′ to mapped low damage reads.
conversely, guanosine and adenosine residues were diminished by ca.
20% and 10% compared to genomic levels at positions immediately 3′ to mapped low damage reads.
the same pattern was seen in mapped high-damage reads, but at ca.
60% and 30% differences from reference genome levels.
for both low and high damage read pools, complementary pyrimidine residues at genomic positions immediately 5′ and 3′ of mapped reads showed a pattern directly opposite of purine residues.
fragmentation bias was not evident in mapped reads from no damage read pools.
misincorporation rates of mapped reads in low and high damage read pools followed expected trends based on patterns of synthetically introduced damage, and were similar to reports from empirical data for 5′ c-t and 3′ g-a transition rates  .
nonetheless, we also found significant differences between expected and observed misincorporation rates.
for example, misincorporation rates of mapped reads were generally lower than expected based on the frequencies of incorporated damage in low and high damage treatments , in some sample replicates dropping to ca.
40% of the expected value at sequence read ends.
while misincorporation frequencies at both 5′ and 3′ ends of reads were significantly and positively correlated to damage level , they were also significantly and negatively correlated to read length at 5′ ends and the 3′-most position of mapped reads .
in addition, 5′ c-t misincorporations were slightly but significantly less common than 3′ g-a misincorporations, with 5′ c-t misincorporations averaging ca.
84-98% of the frequency of 3′ g-a misincorporations across divergence and damage levels .
read mappability
gc content, read length and damage all had significant impacts on the proportion of read pools able to be mapped back to their reference at both levels of divergence, while coverage depth was not significantly correlated for either low or high divergence read pool mapping .
read length clearly had the most dramatic impact on read mappability, however, and decreases in successful read mapping were strongly correlated with increasing read length for all damage levels and gc contents at both levels of divergence .
the rate of decrease in read mapping success was very similar for low and high divergence read pools .
nonetheless, the overall lower success rate of read mapping in high divergence read pools resulted in an approximately three-fold decrease in the percent of reads mapped when read lengths were increased from  <dig> to 80 bp, as opposed to an approximately two-fold decrease observed in low divergence pools.
increasing average read length of a read pool also resulted in a greater decrease in average length of mapped reads across all damage levels and at both low and high divergence , such that average lengths of mapped reads at 40/60/80 bp expected size were  <dig> / <dig> / <dig>  and  <dig> / <dig> / <dig>  bp for low and high divergence read pools, respectively.
the impact of damage levels on read mappability was also clear, albeit lesser than that of read length, as the average percent of reads mapped decreased with increased damage at all read lengths, for both divergence levels .table  <dig>
impact of independent variables on the percent of reads mapped from read pools
divergence
%gc
read length
damage
coverage depth
slope coefficients resulting from multiple variable regression analysis are shown applying best-fit linear models, considering as dependent variables the percent of reads mapped from a pool, while gc content, read length, damage level and coverage depth were treated as independent variables.
values of  <dig>   <dig> and  <dig> were assigned to variables no-, low- and high-damage, respectively.
values not significant at pr <  <dig>  unless otherwise indicated; * signifies  <dig>  < pr < = <dig> ; *** signifies pr < = <dig> .
average proportions of simulated read pools mapped to references across treatments.
a) low divergence level and b) high divergence level.
white, grey and black circles represent no damage, low damage and high damage treatments, respectively.
error bars represent standard deviations of measurements.
for a) and b), best-fit linear relationship is shown for no damage , low damage  and high damage  samples; slopes for linear relationships in a) and b) range from - <dig>  to - <dig> ; r <dig> values range from a)  <dig>  to  <dig>  and b)  <dig>  to  <dig>  .
accuracy of variant calls
both snps  and indels  were generally called with high accuracy across coverage depths, gc contents, read lengths, and damage levels.
indels averaged  <dig> %  and  <dig> %  correct, and snps averaged  <dig> %  and  <dig> %  correct, for low and high divergence treatments, respectively.
nonetheless, snp and indel call accuracy were differentially impacted by damage levels.
the accuracy of indel calls was not substantially affected by damage level at either level of divergence .
in contrast, snp call accuracy was diminished by higher rates of damage, particularly at low divergence .
a significant and positive correlation between read length and variant call accuracy was also demonstrated for both snps and indels in high divergence treatments, and for snps in low divergence treatments, although the impact of this relationship was relatively small .
gc content had significant correlations to variant call accuracy only for snp calls  at high divergence ; the magnitude of this relationship was also relatively small.table  <dig>
impact of independent variables on variant call accuracy and variant call completeness
divergence
%gc
read length
damage
coverage depth
<dig> log***/- <dig> ***
- <dig> ***
<dig> ***
<dig> ***
<dig> log***/- <dig> ***
- <dig> ***
<dig> ***
<dig> ***
slope coefficients resulting from multiple variable regression analysis, considering as dependent variables the percent of indels/snps correctly called  and the percent of total indels/snps called , while gc content, read length, damage level and coverage depth were treated as independent variables.
values of  <dig>   <dig> and  <dig> were assigned to variables no-, low- and high-damage, respectively.
values not significant at pr <  <dig>  unless otherwise indicated; * signifies  <dig>  < pr < = <dig> ; ** signifies  <dig>  < pr < = <dig> ; *** signifies pr < = <dig> .
for each row, the independent variable with the strongest significant effect is highlighted in bold.
results solely for best-fit linear models are shown except for “%indels called correctly” at low and high divergence and “% snps called correctly” at low divergence.
for these cases, a slight improvement in correlation was seen using a logarithmic model, and slope coefficients for both models are shown as: log coefficient / linear coefficient.
effects of damage on accuracy of snp calls across all triplicate-averaged read pools.
box and whiskers plots for a) low and b) high divergence levels.
completeness of variant recovery
in contrast to variant call accuracy, variant call completeness was substantially impacted by coverage depth at both levels of divergence for both snps and indels .
for example, snp calls in low divergence samples varied from ca.
<dig> % complete at  <dig> × coverage depth to nearly  <dig> % complete at 16× coverage depth.
this trend was not quite as pronounced at high divergence, as on average snp recovery ranged from ca.
0%  to ca.
47%  for high divergence read pools.
indel call completeness was lower on average than snp call completeness for low divergence read pools, peaking at ca.
59% completeness for 16× coverage depth, but was nearly indistinguishable from snp call completeness for high divergence read pools.
both snp and indel call completeness at low and high divergence were also significantly and negatively correlated with read length, with the exception of snp call completeness at low divergence  .
however, this impact was more than an order of magnitude smaller than that of coverage depth in all cases.
damage level had a significant and negative influence on variant call completeness for both snps and indels, although this relationship did not have as strong an impact as coverage depth .figure  <dig>
effects of coverage depth at low and high divergence levels on variant call completeness.
a) percent of indels called correctly out of all indels for a reference-sample pair; b) percent of snps called correctly out of all snps for a reference-sample pair.
squares and circles represent data points for low and high divergence samples, respectively.
best-fit linear relationships shown for low  and high  divergence samples for a) and b).
equations and r <dig> values for best-fit linear relationships are a) y =  <dig> x +  <dig> , r2 =  <dig>  ; y =  <dig> x +  <dig> , r2 =  <dig>  ; b) y =  <dig> x +  <dig> , r2 =  <dig>  ; y =  <dig> x +  <dig> , r2 =  <dig>  .
slopes for low versus high divergence samples are not significantly different in a) or b) at p <  <dig> ; however, percentage of variants found regressed on divergence level is significant at p <  <dig>  and p <  <dig>  for a) and b), respectively.
recovery of homozygous versus heterozygous variants
we observed for both snp and indel recovery at both levels of divergence that ratios of homozygous to heterozygous variant calls approached their expected values with increasing coverage depth , but closely approached expected values only for snps at the highest coverage depth tested .
neither gc content, read length or damage level significantly affected ratios of homozygous to heterozygous variant calls, with the exception of a slight negative correlation between read length and snp homozygous to heterozygous ratios at high divergence.
the impact of read length in this case was more than an order of magnitude less than that of coverage depth.figure  <dig>
effects of coverage depth on ratio of homozygous to heterozygous variant calls.
ratios of observed homozygous:heterozygous to true homozygous:heterozygous variant call ratios.
squares and circles represent data points for low and high divergence samples, respectively, for a) indel variant calls and b) snp variant calls.
equations and r <dig> values for best-fit linear relationship are a) y = - <dig>  +  <dig> , r2 =  <dig>  for low divergence samples ; y = - <dig> x +  <dig> , r2 =  <dig>  for high divergence samples ; b) y = - <dig> x +  <dig> , r2 =  <dig>  for low divergence samples ; y = - <dig> x +  <dig> , r2 =  <dig>  for high divergence samples .
dashed grey line added at unity for observed hom:het / true hom:het for visual clarity.
discussion
as adna research broadens its utilization of genome-level sequencing, it is essential to understand how the unique characteristics of adna molecules might impact variant calling from short-read sequence data sets generated by massively parallel sequencing platforms.
in the current study, we quantified the impact of several variables, including the level of divergence between reference and sample sequences, levels of damage in “ancient” sample sequences, and coverage depth on variant call accuracy and completeness.
our use of simulated datasets enabled precise and relatively unbiased measurements of the effects of adna-associated damage on read mapping and variant calling through comparison of known variants to variants called from mapping of short reads with calculated damage patterns.
the variable settings used in this study, including the divergence level between reference and sample sequences, sequencing coverage depth and damage patterns, were chosen based on previously reported adna studies.
for example, in some cases highly similar, same-species reference genomes are available for read mapping for adna samples  <cit> .
our low divergence reference-sample pairs, at approximately  <dig> % divergence, reflect this scenario.
alternatively, in many cases it is necessary to apply more divergent reference genomes from extant genus- or family-level congeners when sequencing extinct taxa  <cit> .
our high divergence reference-sample pairs, at approximately 3% divergence, were used to estimate these effects.
in addition, results from these divergence levels may be used as proxy for predicting variant call success across regions of varying evolutionary rates within the same genome.
importantly, the divergence levels  employed in this study also capture the challenges associated with mapping reads harbouring multiple variant positions, as sequence reads from high divergence read pools averaged more than one variant position per read length.
gc content, while not expected to have an effect on variant calling a priori, was also varied in order to reflect the range of conditions that could be encountered across a nuclear genome, or when sequencing organellar, prokaryotic or viral genomes.
levels of fragmentation, fragmentation biases, and misincorporation rates used in the present study were also based on previously reported adna damage levels.
fragmentation is a universal process associated with post-mortem dna degradation  <cit> , with the result that average read lengths reported in adna studies are typically under 100 bp, and can range to under 50 bp  <cit> .
in turn, nucleotide frequencies at positions immediately adjacent to fragment ends can shift by approximately 10% to over 70% of expected values due to fragmentation bias .
although it may seem intuitive that levels of fragmentation and fragmentation bias would increase with sample age and level of damage, as misincorporation rates do, these relationships do not hold consistently, and may in some cases be reversed  <cit> .
nonetheless, our goal was to both encapsulate typical variation in adna damage levels, as well as model “worse-case” scenarios for comparative purposes.
similarly, the misincorporation frequencies used in this study were within the range of reported values for a range of adna samples.
misincorporation frequencies for late pleistocene neandertal, cave bear and mammoth sequences have been reported in the range of ca.
<dig> - <dig>  at sequence read ends  <cit> , while over  <dig>  for a holocene human sample  <cit> , and as low as ca.
<dig> - <dig>  for a quagga museum specimen and a late pleistocene hippidion fossil  <cit> .
importantly, although our damage levels fell within those previously reported in a range of studies, our data suggest that fragmentation and deamination patterns of adna molecules can result in misrepresentation of the characteristics of adna read pools as reported solely from mapped reads.
for example, read length, deamination rates at fragment ends, and increased divergence to reference all substantially impacted the efficiency of read mapping.
as a result, the proportion of sequence reads mapped successfully from a given read pool ranged from over 90% to under 20% , which in turn would result in up to ca.
five-fold under-estimation of the endogenous content of a read pool.
similarly, misincorporation rates recorded from mapped read pools tended to be lower than values set in read pools prior to mapping, and across treatments were significantly and negatively correlated to read length .
taken together, these results further suggest that the loss of information due to unsuccessful read mapping of adna fragments can in some cases be substantial.
in our study, this loss was exacerbated by relatively longer read lengths, increased divergence between reference and sample sequences and by increased misincorporation damage near the ends of sequence reads.
it is possible that the utilization of alternative mapping and variant calling pipelines could also impact these results.
in particular, other short read alignment algorithms designed for higher sensitivity, such as soap  <cit> , bwa-mem  or bwa-pssm  <cit> , might improve read mapping or variant call performance, particularly at higher divergence levels  <cit> .
alternatively, adjusting bwa mapping parameters alone could improve read mapping success rates in some cases.
however, exploration of bwa parameters on illumina adna sequence reads has been shown to result in very small gains  in the number of reads mapped to a reference  <cit> , and alignment algorithms beyond bwa have not yet been widely applied to adna studies.
finally, the magnitude of decrease in mapped reads with longer read lengths is likely dependent to some degree on the read-simulation algorithm used in this study  <cit> , and so might also vary with different reference-sample combinations and the sequencing platform employed in adna studies.
our results also suggest that damage patterns associated with adna molecules may often have relatively little impact on variant call accuracy at low to moderate average coverage depths.
simulated indels, and to a slightly lesser extent, snps, were generally called accurately across read lengths and damage levels, and the magnitude of the impact of damage on variant call accuracy was small for indels.
the absolute values for variant call accuracy reported here are also likely affected by the use of simulated data  and the applied variant calling pipeline, and may be higher in biological datasets or with alternative bioinformatic strategies .
however, it is also clear that snp call accuracy is impacted at high damage levels.
this trend is reinforced by the significant and positive correlation between damage level and misincorporation frequency in mapped reads .
it is somewhat surprising that the impact of damage on snp call accuracy was more pronounced at a lower level of divergence .
it is likely in these cases that uneven coverage  coupled with relatively higher levels of misincorporations, resulted in sufficient numbers of reads carrying identical ‘post-mortem’ damage to give erroneous snp calls.
in support of this, we found significant negative correlations  between damage levels and the proportion of mapped reads in a pool with unique 5′ mapping positions at both low and high divergence levels, suggesting that increased fragmentation bias  does lead to uneven distribution of mapped reads across a reference genome.
this trend may be lessened at higher divergence, as the introduction of high levels of ‘post-mortem’ damage into reads already containing multiple variants further decreases the probability of mapping reads with misincorporations.
at the same time, the significant and positive correlations between read length and variant call accuracy for both snps and indels  further underscore the potential for loss of information through unsuccessful mapping of longer and/or damaged reads noted above.
in this regard, it is also worth noting that a high proportion of the lowest triplicate-averaged variant call accuracy scores came from read pools with the shortest average read lengths.
for example, among snp and indel call accuracy scores less than  <dig> percent, approximately two thirds were from read pools with 40 bp average read lengths .
we also found that variant call accuracy and coverage depth were significantly correlated, particularly at high divergence .
unexpectedly, this relationship was negative, suggesting that under the variant calling algorithm applied here, variant call accuracy can suffer slightly as coverage depth increases across low coverage levels.
this impact was relatively minor, however, and was approximately an order of magnitude smaller than concurrent percentage gains in overall variant recovery .
in addition, we found the negative impact of increasing coverage depth on variant call accuracy was diminished when considering coverage depths greater than 1×, suggesting that optimal gains in variant calls may be encountered by increasing sample coverage from low to moderate levels as opposed to, for example, increasing coverage depths from ultra-low to low levels.
similar to variant call accuracy, overall variant recovery was not strongly impacted by damage levels, although we did find significant and negative relationships between damage and variant call completeness for both snps and indels .
in contrast, coverage depth clearly had a substantial impact on the completeness of variant recovery, with recovery increasing by ca.
<dig> - <dig> × for snps and indels for each 1× increase in coverage depth .
slight negative correlations between variant call completeness and read length  again likely reflect the difficulty in mapping longer reads, which results in lower coverage depth and decreased variant recovery.
in addition, divergence clearly impacted variant call completeness for both snps and indels , in particular at high coverage levels, suggesting that the availability of a closely related reference may be of substantial benefit in this regard.
as a general trend, the lower overall recovery rate for indels compared to snps evident in our low divergence results is not unexpected.
the combination of short read lengths and the requirement for efficiency in alignment algorithms leads to a general bias in calling small polymorphisms, such that confidently called snps are overrepresented relative to small insertions and deletions  <cit> .
at the same time, small indels have been implicated in phenotype development in a range of organisms  <cit> , and may play a stronger role than snps in early divergence of closely related taxa  <cit> .
considering this, and the fact that phenotype interrogation will likely become more prevalent with increased genomic sequencing of ancient samples  <cit> , researchers may consider increasing sequencing efforts beyond that required for efficient snp recovery in order to more fully capture small indels present in ancient genomes.
it is worth noting, however, that a logarithmic correlation between variant call completeness and sequencing coverage depth was nearly as well supported as the best-fitting linear relationships.
this suggests that there may be a point of diminishing returns for variant recovery in many cases.
finally, relative recovery rates of homozygous and heterozygous variants were significantly impacted by coverage depth for both snps and indels, with observed rates more closely reflecting expected rates as coverage depths increased .
nonetheless, ratios of homozygous to heterozygous variants closely approached expected values only for snps at 16× coverage depth.
the trends in homozygous versus heterozygous variant call success in our analyses suggest that the application of population demographic algorithms based on patterns of heterozygosity  may be challenging for adna samples at low average coverage depths.
CONCLUSIONS
overall, the trends in our data are promising for adna studies, as increasing sequencing capacity of high-throughput platforms and decreased per bp sequencing costs  <cit>  will continue to enable both broader and deeper sequencing of ancient samples.
in addition, targeted enrichment strategies, either focused on sub-genomic targets  <cit>  or entire genomes  <cit> , may further increase sequencing efficiency and enable greater per-sample coverage depths.
nonetheless, it is clear that under some circumstances variant calling is significantly impacted by post-mortem damage patterns and low dna availability typical of adna samples, particularly so for accuracy in snp calls at high levels of damage, and for overall indel and snp recovery at low coverage depths.
further, evolutionary analyses requiring accurate identification of homozygous versus heterozygous variants are likely to suffer when incorporating low coverage adna variant calls.
considering that adna-associated damage patterns are essentially unavoidable yet quantifiable features of ancient samples, it is recommended that researchers carefully consider read mapping and variant calling strategies individually for different projects or even for different samples within a project.
further, it is likely that relatively recent efforts in recalibration software based on adna damage patterns  <cit>  may also aid in increasing variant call accuracy, and ideally will continue to be developed and find broad application in the future.

 the field of ancient dna has reached 30 years of age, in which time it has progressed from the amplification of fragments of single loci to the complete sequencing of individual complex genomes. in its first two decades, adna research was primarily focused on pcr-based amplification and subsequent sanger sequencing of selected loci and organellar genomes, with results applied to analyses of population differentiation and phylogeography, phylogenetics, and even metagenomics. in the last decade the field of adna has moved from the genetic to the genomic level with the advent of massively parallel sequencing platforms. this has been accompanied by a concurrent shift in focus to full genome sequencing and assembly, and genome-scale analyses of population trends. as this field enters its fourth decade, it is inevitable that application of genome-level sequencing will become more commonplace. this will facilitate the broader development of ancient population genomics and phylogenomics approaches in a wider range of taxa. nonetheless, in order to fully benefit from genome-scale approaches, it is critical to understand how the unique characteristics of adna might impact results from high-throughput sequencing and read-mapping technologies. in general, there are several important characteristics that distinguish ancient from modern dna samples. first, not only are adna molecules highly fragmented, but the fragmentation process itself is biased toward breakpoints bordered by 5′ purine residues which could result in biased coverage patterns. further, deamination of cytosine residues at or near fragment ends ultimately leads to c-t and g-a misincorporations at 5′ and 3′ sequence ends, respectively, and could confound adna read mapping and variant calling. last, but not least, ancient samples typically feature relatively low proportions of endogenous dna such that sequencing multiple samples to higher coverage depths to enable confident genomic assembly and variant calling, is often either cost-prohibitive or simply not possible due to limited sample or dna availability. as a consequence of this, most reports applying high throughput sequencing to full genomes of ancient samples have been limited to one or several samples sequenced and at low or ultra-low to moderate coverage depths. this combination can severely impact the accuracy of variant calls. while it is relatively simple to document the biases associated with mapping adna sequence reads, it is more difficult to measure the impacts of these biases on variant calling when analysing empirical data. considering this, we designed a study to quantify the impact of adna fragmentation and misincorporation biases on variant calling using simulated data. by utilizing simulated datasets, we were able to control for a number of potentially confounding factors that commonly feature in and could have an effect upon adna studies. these include level of divergence between sample and reference, average read length, average sequencing depth and damage level. in addition, by synthesizing in silico reference and sample sequences prior to incorporation of adna-associated damage patterns, it is possible to know the exact “pre-damage” nucleotide sequences of samples. this allows direct comparison of called variants to the true variants, in terms of both accuracy and overall completeness.  to explore the impacts of post-mortem dna damage on variant calling, we first simulated triplicate pairs of reference and diploid sample sequences at each of three gc contents and two divergence levels, and at ca. mbp lengths, as described in the methods. short read sequence pools were artificially generated from the sample sequences at average read lengths of and 80 bp, and were further subjected to either zero, low or high levels of synthetic post-mortem damage reflective of previous reports. this damage included both fragmentation bias and elevated 5′ c-t/3′ g-a transitions near sequence ends. these pools were subsequently mapped back to their corresponding reference sequences using commonly applied alignment software and at a range of ultra-low to moderate sequencing coverage depths. this method allowed us to directly compare known damage patterns and variant positions to those measured from mapped read pools, and ultimately to estimate whether and to what extent a variety of variables, including coverage depth, read length and damage level, might impact upon variant calling. synthetic damage patterns implemented patterns of fragmentation bias were clearly seen when mapping reads with both low and high levels of synthetic damage. in low damage read pools, 5′ guanosine and adenosine residues were elevated by ca. 20% and 10%, respectively, over genomic levels at positions immediately 5′ to mapped low damage reads. conversely, guanosine and adenosine residues were diminished by ca. 20% and 10% compared to genomic levels at positions immediately 3′ to mapped low damage reads. the same pattern was seen in mapped high-damage reads, but at ca. 60% and 30% differences from reference genome levels. for both low and high damage read pools, complementary pyrimidine residues at genomic positions immediately 5′ and 3′ of mapped reads showed a pattern directly opposite of purine residues. fragmentation bias was not evident in mapped reads from no damage read pools. misincorporation rates of mapped reads in low and high damage read pools followed expected trends based on patterns of synthetically introduced damage, and were similar to reports from empirical data for 5′ c-t and 3′ g-a transition rates. nonetheless, we also found significant differences between expected and observed misincorporation rates. for example, misincorporation rates of mapped reads were generally lower than expected based on the frequencies of incorporated damage in low and high damage treatments, in some sample replicates dropping to ca. 40% of the expected value at sequence read ends. while misincorporation frequencies at both 5′ and 3′ ends of reads were significantly and positively correlated to damage level, they were also significantly and negatively correlated to read length at 5′ ends and the 3′-most position of mapped reads. in addition, 5′ c-t misincorporations were slightly but significantly less common than 3′ g-a misincorporations, with 5′ c-t misincorporations averaging ca. 84-98% of the frequency of 3′ g-a misincorporations across divergence and damage levels. read mappability gc content, read length and damage all had significant impacts on the proportion of read pools able to be mapped back to their reference at both levels of divergence, while coverage depth was not significantly correlated for either low or high divergence read pool mapping. read length clearly had the most dramatic impact on read mappability, however, and decreases in successful read mapping were strongly correlated with increasing read length for all damage levels and gc contents at both levels of divergence. the rate of decrease in read mapping success was very similar for low and high divergence read pools. nonetheless, the overall lower success rate of read mapping in high divergence read pools resulted in an approximately three-fold decrease in the percent of reads mapped when read lengths were increased from to 80 bp, as opposed to an approximately two-fold decrease observed in low divergence pools. increasing average read length of a read pool also resulted in a greater decrease in average length of mapped reads across all damage levels and at both low and high divergence, such that average lengths of mapped reads at 40/60/80 bp expected size were / / and / /  bp for low and high divergence read pools, respectively. the impact of damage levels on read mappability was also clear, albeit lesser than that of read length, as the average percent of reads mapped decreased with increased damage at all read lengths, for both divergence levels.table impact of independent variables on the percent of reads mapped from read pools divergence %gc read length damage coverage depth slope coefficients resulting from multiple variable regression analysis are shown applying best-fit linear models, considering as dependent variables the percent of reads mapped from a pool, while gc content, read length, damage level and coverage depth were treated as independent variables. values of and were assigned to variables no-, low- and high-damage, respectively. values not significant at pr <  unless otherwise indicated; * signifies  < pr < =; *** signifies pr < =. average proportions of simulated read pools mapped to references across treatments. a) low divergence level and b) high divergence level. white, grey and black circles represent no damage, low damage and high damage treatments, respectively. error bars represent standard deviations of measurements. for a) and b), best-fit linear relationship is shown for no damage, low damage and high damage samples; slopes for linear relationships in a) and b) range from - to -; r values range from a) to and b) to. accuracy of variant calls both snps and indels were generally called with high accuracy across coverage depths, gc contents, read lengths, and damage levels. indels averaged % and % correct, and snps averaged % and % correct, for low and high divergence treatments, respectively. nonetheless, snp and indel call accuracy were differentially impacted by damage levels. the accuracy of indel calls was not substantially affected by damage level at either level of divergence. in contrast, snp call accuracy was diminished by higher rates of damage, particularly at low divergence. a significant and positive correlation between read length and variant call accuracy was also demonstrated for both snps and indels in high divergence treatments, and for snps in low divergence treatments, although the impact of this relationship was relatively small. gc content had significant correlations to variant call accuracy only for snp calls at high divergence; the magnitude of this relationship was also relatively small.table impact of independent variables on variant call accuracy and variant call completeness divergence %gc read length damage coverage depth log***/- *** - *** *** *** log***/- *** - *** *** *** slope coefficients resulting from multiple variable regression analysis, considering as dependent variables the percent of indels/snps correctly called and the percent of total indels/snps called, while gc content, read length, damage level and coverage depth were treated as independent variables. values of and were assigned to variables no-, low- and high-damage, respectively. values not significant at pr <  unless otherwise indicated; * signifies  < pr < =; ** signifies  < pr < =; *** signifies pr < =. for each row, the independent variable with the strongest significant effect is highlighted in bold. results solely for best-fit linear models are shown except for “%indels called correctly” at low and high divergence and “% snps called correctly” at low divergence. for these cases, a slight improvement in correlation was seen using a logarithmic model, and slope coefficients for both models are shown as: log coefficient / linear coefficient. effects of damage on accuracy of snp calls across all triplicate-averaged read pools. box and whiskers plots for a) low and b) high divergence levels. completeness of variant recovery in contrast to variant call accuracy, variant call completeness was substantially impacted by coverage depth at both levels of divergence for both snps and indels. for example, snp calls in low divergence samples varied from ca. % complete at × coverage depth to nearly % complete at 16× coverage depth. this trend was not quite as pronounced at high divergence, as on average snp recovery ranged from ca. 0% to ca. 47% for high divergence read pools. indel call completeness was lower on average than snp call completeness for low divergence read pools, peaking at ca. 59% completeness for 16× coverage depth, but was nearly indistinguishable from snp call completeness for high divergence read pools. both snp and indel call completeness at low and high divergence were also significantly and negatively correlated with read length, with the exception of snp call completeness at low divergence. however, this impact was more than an order of magnitude smaller than that of coverage depth in all cases. damage level had a significant and negative influence on variant call completeness for both snps and indels, although this relationship did not have as strong an impact as coverage depth.figure effects of coverage depth at low and high divergence levels on variant call completeness. a) percent of indels called correctly out of all indels for a reference-sample pair; b) percent of snps called correctly out of all snps for a reference-sample pair. squares and circles represent data points for low and high divergence samples, respectively. best-fit linear relationships shown for low and high divergence samples for a) and b). equations and r values for best-fit linear relationships are a) y =  x + , r2 = ; y =  x + , r2 = ; b) y =  x + , r2 = ; y =  x + , r2 = . slopes for low versus high divergence samples are not significantly different in a) or b) at p < ; however, percentage of variants found regressed on divergence level is significant at p <  and p <  for a) and b), respectively. recovery of homozygous versus heterozygous variants we observed for both snp and indel recovery at both levels of divergence that ratios of homozygous to heterozygous variant calls approached their expected values with increasing coverage depth, but closely approached expected values only for snps at the highest coverage depth tested. neither gc content, read length or damage level significantly affected ratios of homozygous to heterozygous variant calls, with the exception of a slight negative correlation between read length and snp homozygous to heterozygous ratios at high divergence. the impact of read length in this case was more than an order of magnitude less than that of coverage depth.figure effects of coverage depth on ratio of homozygous to heterozygous variant calls. ratios of observed homozygous:heterozygous to true homozygous:heterozygous variant call ratios. squares and circles represent data points for low and high divergence samples, respectively, for a) indel variant calls and b) snp variant calls. equations and r values for best-fit linear relationship are a) y = -  + , r2 =  for low divergence samples; y = - x + , r2 =  for high divergence samples; b) y = - x + , r2 =  for low divergence samples; y = - x + , r2 =  for high divergence samples. dashed grey line added at unity for observed hom:het / true hom:het for visual clarity. discussion as adna research broadens its utilization of genome-level sequencing, it is essential to understand how the unique characteristics of adna molecules might impact variant calling from short-read sequence data sets generated by massively parallel sequencing platforms. in the current study, we quantified the impact of several variables, including the level of divergence between reference and sample sequences, levels of damage in “ancient” sample sequences, and coverage depth on variant call accuracy and completeness. our use of simulated datasets enabled precise and relatively unbiased measurements of the effects of adna-associated damage on read mapping and variant calling through comparison of known variants to variants called from mapping of short reads with calculated damage patterns. the variable settings used in this study, including the divergence level between reference and sample sequences, sequencing coverage depth and damage patterns, were chosen based on previously reported adna studies. for example, in some cases highly similar, same-species reference genomes are available for read mapping for adna samples. our low divergence reference-sample pairs, at approximately % divergence, reflect this scenario. alternatively, in many cases it is necessary to apply more divergent reference genomes from extant genus- or family-level congeners when sequencing extinct taxa. our high divergence reference-sample pairs, at approximately 3% divergence, were used to estimate these effects. in addition, results from these divergence levels may be used as proxy for predicting variant call success across regions of varying evolutionary rates within the same genome. importantly, the divergence levels employed in this study also capture the challenges associated with mapping reads harbouring multiple variant positions, as sequence reads from high divergence read pools averaged more than one variant position per read length. gc content, while not expected to have an effect on variant calling a priori, was also varied in order to reflect the range of conditions that could be encountered across a nuclear genome, or when sequencing organellar, prokaryotic or viral genomes. levels of fragmentation, fragmentation biases, and misincorporation rates used in the present study were also based on previously reported adna damage levels. fragmentation is a universal process associated with post-mortem dna degradation, with the result that average read lengths reported in adna studies are typically under 100 bp, and can range to under 50 bp. in turn, nucleotide frequencies at positions immediately adjacent to fragment ends can shift by approximately 10% to over 70% of expected values due to fragmentation bias. although it may seem intuitive that levels of fragmentation and fragmentation bias would increase with sample age and level of damage, as misincorporation rates do, these relationships do not hold consistently, and may in some cases be reversed. nonetheless, our goal was to both encapsulate typical variation in adna damage levels, as well as model “worse-case” scenarios for comparative purposes. similarly, the misincorporation frequencies used in this study were within the range of reported values for a range of adna samples. misincorporation frequencies for late pleistocene neandertal, cave bear and mammoth sequences have been reported in the range of ca. - at sequence read ends, while over for a holocene human sample, and as low as ca. - for a quagga museum specimen and a late pleistocene hippidion fossil. importantly, although our damage levels fell within those previously reported in a range of studies, our data suggest that fragmentation and deamination patterns of adna molecules can result in misrepresentation of the characteristics of adna read pools as reported solely from mapped reads. for example, read length, deamination rates at fragment ends, and increased divergence to reference all substantially impacted the efficiency of read mapping. as a result, the proportion of sequence reads mapped successfully from a given read pool ranged from over 90% to under 20%, which in turn would result in up to ca. five-fold under-estimation of the endogenous content of a read pool. similarly, misincorporation rates recorded from mapped read pools tended to be lower than values set in read pools prior to mapping, and across treatments were significantly and negatively correlated to read length. taken together, these results further suggest that the loss of information due to unsuccessful read mapping of adna fragments can in some cases be substantial. in our study, this loss was exacerbated by relatively longer read lengths, increased divergence between reference and sample sequences and by increased misincorporation damage near the ends of sequence reads. it is possible that the utilization of alternative mapping and variant calling pipelines could also impact these results. in particular, other short read alignment algorithms designed for higher sensitivity, such as soap, bwa-mem or bwa-pssm, might improve read mapping or variant call performance, particularly at higher divergence levels. alternatively, adjusting bwa mapping parameters alone could improve read mapping success rates in some cases. however, exploration of bwa parameters on illumina adna sequence reads has been shown to result in very small gains in the number of reads mapped to a reference, and alignment algorithms beyond bwa have not yet been widely applied to adna studies. finally, the magnitude of decrease in mapped reads with longer read lengths is likely dependent to some degree on the read-simulation algorithm used in this study, and so might also vary with different reference-sample combinations and the sequencing platform employed in adna studies. our results also suggest that damage patterns associated with adna molecules may often have relatively little impact on variant call accuracy at low to moderate average coverage depths. simulated indels, and to a slightly lesser extent, snps, were generally called accurately across read lengths and damage levels, and the magnitude of the impact of damage on variant call accuracy was small for indels. the absolute values for variant call accuracy reported here are also likely affected by the use of simulated data and the applied variant calling pipeline, and may be higher in biological datasets or with alternative bioinformatic strategies. however, it is also clear that snp call accuracy is impacted at high damage levels. this trend is reinforced by the significant and positive correlation between damage level and misincorporation frequency in mapped reads. it is somewhat surprising that the impact of damage on snp call accuracy was more pronounced at a lower level of divergence. it is likely in these cases that uneven coverage coupled with relatively higher levels of misincorporations, resulted in sufficient numbers of reads carrying identical ‘post-mortem’ damage to give erroneous snp calls. in support of this, we found significant negative correlations between damage levels and the proportion of mapped reads in a pool with unique 5′ mapping positions at both low and high divergence levels, suggesting that increased fragmentation bias does lead to uneven distribution of mapped reads across a reference genome. this trend may be lessened at higher divergence, as the introduction of high levels of ‘post-mortem’ damage into reads already containing multiple variants further decreases the probability of mapping reads with misincorporations. at the same time, the significant and positive correlations between read length and variant call accuracy for both snps and indels further underscore the potential for loss of information through unsuccessful mapping of longer and/or damaged reads noted above. in this regard, it is also worth noting that a high proportion of the lowest triplicate-averaged variant call accuracy scores came from read pools with the shortest average read lengths. for example, among snp and indel call accuracy scores less than percent, approximately two thirds were from read pools with 40 bp average read lengths. we also found that variant call accuracy and coverage depth were significantly correlated, particularly at high divergence. unexpectedly, this relationship was negative, suggesting that under the variant calling algorithm applied here, variant call accuracy can suffer slightly as coverage depth increases across low coverage levels. this impact was relatively minor, however, and was approximately an order of magnitude smaller than concurrent percentage gains in overall variant recovery. in addition, we found the negative impact of increasing coverage depth on variant call accuracy was diminished when considering coverage depths greater than 1×, suggesting that optimal gains in variant calls may be encountered by increasing sample coverage from low to moderate levels as opposed to, for example, increasing coverage depths from ultra-low to low levels. similar to variant call accuracy, overall variant recovery was not strongly impacted by damage levels, although we did find significant and negative relationships between damage and variant call completeness for both snps and indels. in contrast, coverage depth clearly had a substantial impact on the completeness of variant recovery, with recovery increasing by ca. - × for snps and indels for each 1× increase in coverage depth. slight negative correlations between variant call completeness and read length again likely reflect the difficulty in mapping longer reads, which results in lower coverage depth and decreased variant recovery. in addition, divergence clearly impacted variant call completeness for both snps and indels, in particular at high coverage levels, suggesting that the availability of a closely related reference may be of substantial benefit in this regard. as a general trend, the lower overall recovery rate for indels compared to snps evident in our low divergence results is not unexpected. the combination of short read lengths and the requirement for efficiency in alignment algorithms leads to a general bias in calling small polymorphisms, such that confidently called snps are overrepresented relative to small insertions and deletions. at the same time, small indels have been implicated in phenotype development in a range of organisms, and may play a stronger role than snps in early divergence of closely related taxa. considering this, and the fact that phenotype interrogation will likely become more prevalent with increased genomic sequencing of ancient samples, researchers may consider increasing sequencing efforts beyond that required for efficient snp recovery in order to more fully capture small indels present in ancient genomes. it is worth noting, however, that a logarithmic correlation between variant call completeness and sequencing coverage depth was nearly as well supported as the best-fitting linear relationships. this suggests that there may be a point of diminishing returns for variant recovery in many cases. finally, relative recovery rates of homozygous and heterozygous variants were significantly impacted by coverage depth for both snps and indels, with observed rates more closely reflecting expected rates as coverage depths increased. nonetheless, ratios of homozygous to heterozygous variants closely approached expected values only for snps at 16× coverage depth. the trends in homozygous versus heterozygous variant call success in our analyses suggest that the application of population demographic algorithms based on patterns of heterozygosity may be challenging for adna samples at low average coverage depths.  overall, the trends in our data are promising for adna studies, as increasing sequencing capacity of high-throughput platforms and decreased per bp sequencing costs will continue to enable both broader and deeper sequencing of ancient samples. in addition, targeted enrichment strategies, either focused on sub-genomic targets or entire genomes, may further increase sequencing efficiency and enable greater per-sample coverage depths. nonetheless, it is clear that under some circumstances variant calling is significantly impacted by post-mortem damage patterns and low dna availability typical of adna samples, particularly so for accuracy in snp calls at high levels of damage, and for overall indel and snp recovery at low coverage depths. further, evolutionary analyses requiring accurate identification of homozygous versus heterozygous variants are likely to suffer when incorporating low coverage adna variant calls. considering that adna-associated damage patterns are essentially unavoidable yet quantifiable features of ancient samples, it is recommended that researchers carefully consider read mapping and variant calling strategies individually for different projects or even for different samples within a project. further, it is likely that relatively recent efforts in recalibration software based on adna damage patterns may also aid in increasing variant call accuracy, and ideally will continue to be developed and find broad application in the future. 
we simulated pairs of sequences to act as reference and sample genomes at varied gc contents and divergence levels. massively parallel sequencing platforms, featuring high throughput and relatively short read lengths, are well suited to ancient dna studies. mapping of sample read pools to reference sequences revealed several trends, including decreased alignment success with increased read length and decreased variant recovery with increased divergence. modest increases in sequencing depth resulted in rapid gains in total variant recovery, and limited improvements to recovery of heterozygous variants. this in silico study suggests adna-associated damage patterns minimally impact variant call accuracy and recovery from short-read alignment, while modest increases in sequencing depth can greatly improve variant recovery. short-read sequence pools were generated from sample sequences, and subjected to varying levels of “post-mortem” damage by adjusting levels of fragmentation and fragmentation biases, transition rates at sequence ends, and sequencing depths. variants were generally called with high accuracy, however identification of snps was less accurate for high damage/low divergence samples. variant identification from short-read alignment could be hindered, however, by low dna concentrations common to historic samples, which constrain sequencing depths, and post-mortem dna damage patterns.
 massively parallel sequencing platforms, featuring high throughput and relatively short read lengths, are well suited to ancient dna studies. variant identification from short-read alignment could be hindered, however, by low dna concentrations common to historic samples, which constrain sequencing depths, and post-mortem dna damage patterns.  we simulated pairs of sequences to act as reference and sample genomes at varied gc contents and divergence levels. short-read sequence pools were generated from sample sequences, and subjected to varying levels of “post-mortem” damage by adjusting levels of fragmentation and fragmentation biases, transition rates at sequence ends, and sequencing depths. mapping of sample read pools to reference sequences revealed several trends, including decreased alignment success with increased read length and decreased variant recovery with increased divergence. variants were generally called with high accuracy, however identification of snps was less accurate for high damage/low divergence samples. modest increases in sequencing depth resulted in rapid gains in total variant recovery, and limited improvements to recovery of heterozygous variants.  this in silico study suggests adna-associated damage patterns minimally impact variant call accuracy and recovery from short-read alignment, while modest increases in sequencing depth can greatly improve variant recovery. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords ancient dnavariant callingnext-generation sequencingcoverage depthissue-copyright-statement© the author 2015 