Relative performance of two DNA extraction and library preparation methods on archaeological human teeth samples

ABSTRACT DNA extraction and library preparation are crucial steps in any ancient DNA study. Although palaeogenomic researchers are facing a growing choice of DNA extraction and sequencing library preparation methods, how their performance varies with DNA preservation remains unclear. To help elucidate this question, we compared the performance of two common DNA extraction and Illumina library preparation methods on a set of archaeological human samples, considered to contain ancient DNA of intermediate to good preservation (5–50% endogenous DNA). Results indicate that while the levels of contamination and endogenous DNA recovered are comparable for both silica-in-solution and silica-column based extractions, the ability of the former to accommodate larger starting quantities of sample material confers notable benefits with regards to library complexity, and furthermore seems to aid with the recovery of shorter endogenous DNA molecules. While our observations gained from comparing the single-stranded with double-stranded DNA library construction methods largely replicate earlier observations, the combination of our data with previously published datasets demonstrate that the benefits gained using single-stranded methods are inversely proportional to the endogenous DNA content in the ancient sample.


Introduction
The research potential offered by ancient DNA (aDNA) has expanded enormously over the last decade thanks to the introduction of High-Throughput Sequencing (HTS) technologies and associated developments in laboratory methodologies dedicated to the generation of palaeogenomic data from archaeological and historic samples. While these innovations have revolutionized the field by enabling genomicscale analyses (Shapiro and Hofreiter 2014;Der Sarkissian et al. 2015;Orlando et al. 2015), optimization of all steps in the data generation and analytical process can further reduce costs and increase data output.
DNA extraction and HTS library preparation methods stand at the beginning of any palaeogenomic study. In both cases, ideal methods should maximize the potential for sequencing endogenous DNA, at the exclusion of contaminants and exogenous sources of DNA. Although several promising protocols have been developed to address these goals (Yang et al. 1998;Rohland and Hofreiter 2007;Meyer et al. 2012;Dabney et al. 2013), validation of their efficacy has been principally limited to a small number of samples with DNA of either very low (<1%) or very high (>30%) endogenous content (Meyer et al. 2012;Bennett et al. 2014;Gamba et al. 2016; Barlow et al. 2016). Therefore, as argued by Bennett et al. (2014) and Barlow et al. (2016), there is a need to expand the comparative performance of different methods across a wider range of endogenous DNA preservation and investigate the different effects that methods have on data generation.
To help address this, we compared the performance of two common aDNA extraction and Illumina library preparation methods ( Figure 1) using a sample set of archaeological human teeth that, based on prior PCR-based analyses, were suspected to contain relatively high quality ancient DNA (>5% endogenous content) (Gilbert et al. 2006). The efficacy and effects of the extraction and library preparation methods was assessed with regards to endogenous DNA content (before and after mitochondrial DNA enrichment), sequence read clonality, average read length, contamination, microbial composition of unmapped reads, and time-dependent damage patterns.

Samples
Eleven human tooth samples were used to investigate the performance of two DNA extraction and Illumina library preparation methodologies. All samples corresponded to premolar teeth from individuals excavated in Trondheim, Norway between 1984-85, and were context dated to 500-1000 years BP (Long, C. D., 1975;Turner-Walker and Syversen 2001).
DNA extraction, library preparation and the first PCR amplification set-up (indexing PCR) were performed in dedicated clean labs separated from post-PCR laboratories at the Centre for GeoGenetics in Copenhagen, Denmark. Protocols followed established procedures to prevent contamination, including the cleaning of the sample surface with 10% bleach solution, UV-irradiation for 5 min on each side of the tooth, and the use of indexed primers during library amplification.

Experimental Setup -Extraction Test
The extraction comparison was performed on six teeth and had two main objectives: 1) to explore the relative success of a silica-column (based on Yang et al. (1998), full details below) versus silica-in-solution method (based on Rohland and Hofreiter (2007), full details below), and 2) to assess the microbial and contaminant fraction in two distinct tooth sections, the inner versus the outer root.
Both extraction methods make use of the properties of silica particles to bind to DNA under certain pH and salt conditions. The volume of the binding reaction and the chaotropic salt used are the main differences between the methods. While the silica-column method (following Yang et al. 1998) uses 5X volume of binding buffer based on guanidine hydrochloride, the silica-insolution method uses 4X volume of binding buffer based on guanidine thiocyanate. The differential effect of the two different chaotropic salts for DNA isolation was not investigated in this study.
Two subsamples of material were taken from each of the six teeth. The first subsample was obtained through drilling into the pulp with a Dremel® drill at 100 rpm, while the second subsample was obtained through complete pulverization of the remainder of the root using a Micro Dismembrator S (Braun Biotech, Melsungen, Germany) at 2000rpm. The first subsample was subject to the column-based extraction, the second to the silica-in-solution extraction. Post extraction, all these samples were converted into DNA libraries following a double-stranded DNA protocol (Meyer and Kircher 2010) as detailed below, then both directly shotgun sequenced as well as subject to target enrichment for the complete mitochondrial genome using an in-solution hybridization capture method as in Maricic et al. (2010).

Silica-column extraction
Tooth powder was digested in 1.3 ml of an EDTAbased digestion buffer containing 0.25 mg/mL Proteinase K and 10 M Urea and concentrated in ∼100 µl using an Amicon Ultra 4 (30 kDa) filter. DNA was then isolated using MinElute spin columns (Qiagen, Hilden, Germany) following the manufacturer's instructions.

Silica-in-solution extraction
Tooth powder was digested in 5 ml of an EDTA-based digestion buffer containing 0.25 mg/mL Proteinase K and 10% N-laurylsarcosyl, then purified using the silica-in-solution (SIS) based approach modified from Rohland and Hofreiter (2007) as follows. After digestion, samples were spun down for 5 min at 2000 rpm and the supernatant was transferred into 20 mL of binding buffer containing 5 M GuSCN, 25 mM NaCl, 50 mM Tris, 20 mM EDTA, and 1.3% TritonX100 (pH 4.0-5.0) and 50 µl of silica suspension. DNA was purified by binding to the silica for 3 h at room temperature with agitation, followed by two washing steps with 80% ethanol. Samples were eluted in 60 μl EB buffer. The amount of tooth powder used for each extraction is shown in table S.1.

Experimental Setup -Library Test
The library comparison was performed on 5 teeth, whose entire roots were ground to a powder, after which DNA was extracted using the above-described silica-in-solution method, then converted into both double-stranded (Meyer and Kircher 2010) and single-stranded libraries (Gansauge and Meyer 2013). As with the extraction test, the DNA content of these libraries was assessed through both direct shotgun, as well as mitochondrial target enrichment sequencing.

Library preparation
Double-stranded DNA (dsDNA) libraries were prepared using a modified protocol based around the NEB blunt-end library preparation kit (Cat. No. E6070,New England Biolabs,Ipswich,Massachusetts) and Illumina adapters modified following (Meyer and Kircher, 2010). Specifically, the initial nebulization step was skipped because of the fragmented nature of aDNA. End-repair was performed in 50 μl reactions with 42.5 μl of DNA extract. The end-repair cocktail was incubated for 20 min at 12°C and 15 min at 37°C and purified using MinElute silica spin columns (Qiagen) following the manufacturer's instructions and eluted in 30 μl of buffer EB. After end-repair the adaptors were ligated to the end-repaired DNA in 50 μl reactions. The reactions were incubated for 15 min at 20°C and purified using QiaQuick columns (Qiagen) before being eluted in 42 μl EB. The adapter fill-in reaction was performed in a final volume of 50 μl and incubated for 20 min at 37°C, followed by 20 min at 80°C to inactivate the Bst enzyme.

Amplification of Ancient DNA Libraries Index amplification of double-stranded libraries
Double-stranded libraries were PCR amplified and indexed in 100 μl reactions, using 49 μl of library template, AmpliTaq Gold Polymerase 0.1 U/µl, H 2 O, 1X AmpliTaq Gold Buffer, 2.5 mM MgCl 2 , 0.8 mg/ml BSA, 0.25 mM dNTPs, and 0.2 μM forward and reverse indexing primer (Meyer and Kircher 2010). Thermocycling conditions were 12 min at 95°C, followed by 10 cycles of 20 sec at 95°C, 30 sec at 60°C, and 40 sec at 72°C, and a final 5 min elongation step at 72°C. The amplified libraries were then purified following the manufacturer's instructions using QiaQuick PCR Purification Kit (Qiagen) and eluted in 30 µl Buffer EB.

Index amplification of single-stranded libraries
Single-stranded libraries were PCR amplified and indexed in a 50 µl PCR reaction, using 30 µl library template, AmpliTaq Gold Polymerase 0.1 U/µl, H 2 O, 1X Buffer, 2.5 mM MgCl 2 , 0.8 mg/ml BSA, 0.25 mM dNTPs, and 0.2 μM forward and reverse indexing primer (Meyer and Kircher 2010). Thermocycling conditions were 12 min at 95°C, followed by cycles of 20 sec at 95°C, 30 sec at 60°C, and 40 sec at 72°C, and a final 5 min elongation step at 72°C. The amplified libraries were then purified following the manufacturer's instructions using QiaQuick PCR Purification Kit (Qiagen) and eluted in 30 µl Buffer EB. Due to high concentrations of primer-dimers, the amplified libraries were additionally purified using Agencourt AMPure XP beads and eluted in 30 μl EB.

Re-amplification of libraries
In order to obtain sufficient concentration of DNA for sequencing, libraries were re-amplified in 100 μl PCR reactions using 10 μl of library template from above and 1X Phusion® High Fidelity DNA Polymerase Master Mix (New England Biolabs, Ipswich, United Kingdom), 0.8 mg/ml BSA, and 10 μM re-amplification primers IS5 and IS6 (Meyer and Kircher 2010). Thermocycling conditions were 30 sec at 98°C, followed by cycles of 20 sec at 98°C, 30 sec at 60°C, and 30 sec at 72°C, and a final 5 min elongation step at 72°C. The number of cycles given to each sample is shown in Table S.1.

Mitochondrial DNA enrichment
Enrichment for the mitochondrial genome was conducted using DNA bait produced from two longrange PCR products of modern human DNA and the hybridization of individual barcoded libraries, as described in detail in Maricic et al. (2010). Before enrichment, DNA concentration of libraries was quantified using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Between 500 ng and 1 µg of DNA for each library was used to enrich for the mitochondrial genome. After enrichment, libraries were amplified in a 100 μl PCR reactions using 15 µl of enriched library template from above and 1X Phusion® High Fidelity DNA Polymerase Master Mix, 0.8 mg/ml BSA, and 10 μM re-amplification primers IS5 and IS6. Thermocycling conditions were 30 sec at 98°C, followed by 14 cycles of 20 sec at 98°C, 30 sec at 60°C, and 30 sec at 72°C, and a final 5 min elongation step at 72°C. The amplified libraries were then purified following the manufacturer's instructions using QiaQuick PCR Purification Kit (Qiagen) and eluted in 30 µl Buffer EB.

Library Sequencing and Data Analysis
Samples were quantified before sequencing using a Bioanalyzer 2100 and indexed DNA libraries were pooled equimolar and sequenced together with other indexed samples on an Illumina HiSeq 2000. The extraction and mtDNA comparison experiments were sequenced on the Illumina HiSeq 2000 for 100 cycles on single-read mode. The library experiment was sequenced on the Illumina HiSeq 2000 for 95 cycles single-read mode (Figure 2A, 3A). The difference in sequencing cycles between experiments has no impact on the results, observations and conclusions.
Basecalling was performed using the Illumina software CASAVA 1.8.2 and the output fastq files were manipulated using AdapterRemoval (Lindgreen 2012) to remove sequencing adapters as well as N's and low-quality bases in the 3 ′ end of reads. Reads shorter than 30 bp were discarded. After filtering, reads were mapped to the hg18 human reference genome using BWA v0.7.5a-r405 (Li et al. 2009), quality and size filters were then applied to exclude sequences with a mapping quality <25. Duplicate reads were removed with SAMtools rmdup (Li et al. 2009). Fragmentation and damage patterns were calculated and plotted using mapDamage2 (Ginolhac et al. 2011;Jonsson et al. 2013), and library complexity was estimated for each sample using preseq (Version 1.0.1, Daley and Smith 2013). To avoid biases, all the sequencing analyses were done on both the complete set of reads and a random subset of 3 million reads for each experiment, drawn after filtering and before mapping to the human genome.
Mitochondrial contamination estimates were generated following the approach described in Fu et al. (2013). Metagenomic analysis to estimate the species level abundance in each of the samples for the two extraction methods compared was performed using Kraken (Wood and Salzberg 2014). We ran Kraken v0.10.5-beta with default parameters and using the standard Kraken database modified to include the human reference genome, for an estimation of endogenous content.

Extraction comparison
Although the two extraction methods used different subsamples of the teeth, we observed no statistically significant difference in the percentage of unique endogenous DNA reads recovered from the shotgun DNA sequence data (paired t-test = −0.3492, df = 5, p-value = 0.6294) ( Table 1). We hypothesised that the different fractions of the tooth could contain different proportions of microbial contamination (due to for example higher exposure to the environment or DNA from handling on the outside). However, no consistent differences were observed between either fraction, whether for mtDNA contamination estimates or microbial metagenomic signatures (Table 1, Table S.6).
Comparisons of mapped fragments showed a statistically significant difference in endogenous read length distribution, with shorter average endogenous read lengths observed for the SIS method (paired t-test = 3.3333, df = 5, p-value = 0.010) (Figure 2, Table S.2). Furthermore, library complexity (as measured both using the preseq software package and through sequence clonality on data normalised for sequencing depth) was considerably higher for SIS than for silica-column DNA extracts (Figure 2, Table  S.1). This directly correlated with the number of amplification cycles required to amplify the libraries prior to sequencing, and is of course expected given the larger starting tissue volumes that the SIS method enables ( Figure S.1 Figure 2. Frequency of read length distributions for sequenced reads (uniquely mapped) for all samples used in extraction comparison (A), and estimated library complexity curves for the samples in the extraction comparison (B and C). Curves were generated using the preseq software (Daley and Smith 2013). The y-axis shows the estimated amount of reads generated after one lane (180 million reads) and 20 lanes of sequencing (x-axis). Shadowed areas represent confidence intervals.
Overall therefore, the principal advantage of SIS lies in its benefit of enriching for shorter endogenous DNA fractions, and maximising the recovery of ancient template molecules.

Library preparation comparison
We subsequently explored the performance of two common aDNA Illumina library preparation methods across five samples (Table 1). Previous studies have provided evidence suggesting that the ssDNA method improves endogenous DNA yields and library complexity (Meyer et al. 2012;Prüfer et al. 2014;Bennett et al. 2014). This is likely due to its ability to (i) help recover DNA molecules with single-stranded breaks and those with end modifications that are entirely lost during dsDNA construction, and (ii) because it avoids DNA loss during purification steps that are an integral part of the dsDNA construction protocol (Bennett et al. 2014).
Although differences were observed in the percentage of endogenous DNA sequences retrieved for the different library preparation methods (Table 1), at the relatively low sequencing depth of our experiment these did not consistently favour one approach over the other, and differences were not found to be significant (paired t-test = -1.4773, df = 4, p-value = 0.8932). Recently Wales et al. (2015) have suggested that a threshold exists beyond which ssDNA library preparation tends to enrich endogenous content as foldincreases are only observed when dsDNA libraries contained <3% endogenous DNA. Given that endogenous content of the samples tested here ranges between 4.77 and 50.2%, our observations are in accordance to this. We speculate that the observed disparities between samples might be explained by sample-specific differences in DNA damage that relate to the library construction methods, as the ssDNA library protocol is advantageous for samples with strand breaks on one or both strands of the DNA molecule (Gansauge and Meyer 2013). Miscoding lesions observed were consistent with previous observations (Gansauge and Meyer 2013), demonstrating an elevation of C to T substitutions at both ends of sequences derived from ssDNA, while dsDNA showed the reverse complement, with an excess of G to A substitutions in the 3 ′ ends, as reported by Briggs et al. (2007) (Figure S.3).
Consistent with previous observations (Gansauge and Meyer 2013;Bennett et al. 2014;Ávila-Arcos et al. 2015), we observed a major difference with regards to sequence complexity, with the ssDNA producing a higher number of distinct reads when sequenced in depth ( Figure 3, Table 1) and lower levels of clonality ( Figure S.4, Table S.1). Similarly, read length distributions between 30-60 bp were significantly different between library types (paired t-test = Figure 3. Frequency of length distributions for uniquely mapped reads for all samples in in library comparison (A) and estimated library complexity curves (B and C). Curves were generated using preseq (Daley and Smith 2013). The y-axis shows the estimated amount of reads generated after one lane (180 million reads) or 20 lanes of sequencing (x-axis). Shadowed areas represent confidence intervals. 7.2736, df = 4, p-value = 0.0009) (Figure 3) with a greater proportion of smaller fragments being recovered with ssDNA ( Figure 3). Ávila-Arcos et al. (2015) demonstrated how the success of capture-enrichment protocols on aDNA correlates inversely with average size of endogenous DNA molecules, an observation consistent with the documented relationship (using modern DNA) of increased capture efficiency with length of target molecules (Clark et al. 2011;Querfurth et al. 2012). Mitogenome target capture on our libraries replicates this observation, with significantly higher proportions of mtDNA reads recovered for the dsDNA than for the ssDNA library method (paired ttest = 4.4673, df = 4, p-value = 0.0055) (Table S.4).
In summary, while some aspects of our results were consistent with previous observations, we noticeably failed to reproduce observed benefits relating to endogenous DNA levels. To explore this further, we took advantage of two similar datasets published by Bennett et al. (2014) and Wales et al. (2015), which focused on the performance of the two library methods on teeth and bone samples (Figure 4). Although their material spanned a wider range of ages and species than ours, notably their samples contained either much lower, or much higher, endogenous DNA content (<1% or >90%) than those studied here. Addition of their data to our own, enables characterization of the relative performance of the two methods across a wide range of endogenous DNA levels, and highlights the trend observed by Wales et al. (2015) that for teeth and bone samples at least a negative correlation exists between the ratio of endogenous DNA improvement of ssDNA over dsDNA and initial endogenous DNA content. In other words, while the ssDNA method generally yields more endogenous DNA from a given sample, this advantage is more apparent on samples with lower initial endogenous DNA content. Samples with high endogenous content, such as those from more recent time periods or those preserved under favourable conditions (e.g. permafrost), tend to yield DNA libraries of similar quality regardless of the library preparation method.

Conclusion
Although our dataset is relatively small, several clear conclusions can be drawn that we hope will be useful for others when deciding upon extraction and library build strategies. Firstly, although we did not find any clear difference with regards to endogenous DNA recovery, the ability of the SIS method to deal with larger levels of input material has the expected advantage of maximizing endogenous template molecules, thus increasing library complexity and enabling deeper sequencing. We note however, that Dabney et al. (2013) introduced a modification to the silica column method that facilitates handling of larger starting volumes through incorporation of reservoirs onto the spin columns. Recent studies by Gamba et al. (2016) and Barlow et al. (2016) observed that the Dabney method consistently recovered more DNA than the other methods tested, increasing also the molecular complexity of the extract. However, the combination of different extraction methods with different library construction protocols can have different outcomes. For example, Barlow et al. (2016) suggest that the combination of the SIS method with dsDNA libraries generates the largest amount of data.
Secondly, while the endogenous DNA enrichment potential of Gansauge and Meyer's (2013) ssDNA method is excellent when applied to poor quality samples (1-3% endogenous DNA) or for use in very high-depth sequencing, the relative benefit of this method rapidly decreases as material quality increases and/or amount of sequencing undertaken decreases. This observation is pertinent given both the significantly higher economic and labour requirements of ssDNA. New evidence presented by Barlow et al. (2016) suggests that in order to maximize library complexity, a combination of the Dabney method and ssDNA library protocol yields the best results. While many of our results are consistent with those of previous studies, we observe clear exceptions. We believe that discrepancies most likely relate to the fact that at least some of the previous studies have been restricted to either single samples (Meyer et al. 2012) or material that is characterised by extremely low (<1%) endogenous DNA contents (Bennett et al. 2014)quite distinct from the relatively good and intermediate quality material tested here and in recent studies by Barlow et al. (2016) and Gamba et al. (2016).
In agreement with previous methodological comparisons, we believe that, depending on the scope of the study and the priorities regarding the data generation, a careful initial assessment of the available methodologies must be taken by researchers in order to decide which combination of methodologies best suits their aims. Most of all, our findings highlight the danger of generalizing on the efficacy of a particular extraction or library preparation method from a limited number of samples. This in turn underscores the need for pilot studies prior to subjecting irreplaceable material to newly-developed extraction or library preparation methods.