Comparative performance of the BGISEQ-500 vs Illumina HiSeq2500 sequencing platforms for palaeogenomic sequencing

1As with many other disciplines, the advent of next-generation sequencing (NGS) platforms has revolutionized ancient DNA (aDNA) research. During the era of Sanger sequencing, the datasets within most studies were restricted to short lengths of mitochondrial DNA (mtDNA) or nuclear DNA (nuDNA), and at most, if one used multiplexing techniques, one could aim for mitochondrial genomes (mitogenomes). However, thanks to NGS techniques, with the right sample and sufficient funds, today practitioners are able to aim for relatively complete ancient nuclear genomes (hereafter referred to as palaeogenomes), even at the population level. While there are now a range of NGS technologies available to choose from, those favoured by the aDNA field are suited to the characteristically short DNA molecules that dominate aDNA extracts; thus long-read technologies such as the PacBio (Pacific Biosciences, CA, USA) and Minion (Oxford Nanopore Technologies, Oxford, UK) are not widely used. A range of technologies have been explored in the aDNA context, including the Roche/454 series, SOLID-4, the now discontinued Helicos, and the Ion Torrent series. The undisputed workhorses, however, are the platforms within the Illumina series, principally due to a combination of factors that include sequencing cost per unit date, low sequencing error rate (ca 0.1%), and simply the number of machines available upon which to sequence. Thus in recent years, the focus has been placed on the development and optimization of methods in order to increase data quality and reduce overall cost. Steps taken have included both tailoring library constructions and amplification methods toward the damaged endogenous aDNA, e.g., through exploiting blunt end ligations, removing steps associated with DNA loss, enzyme choice, or even focusing on direct ligation to single-stranded DNA, as well as improvement in the informatic tools that are used to process the Fastq data generated.

2Today, therefore, Illumina-based sequencing has formed the basis of the overwhelming majority of palaeogenomic studies, including (but not limited to) draft genomes of humans and related hominids, animals, plants and even pathogens, population genomic datasets, metagenomic studies, and even insights into ancient transcriptomes and epigenomes. For recent reviews, see.

3Despite this remarkable progress, palaeogenomics still faces one significant limitation—the overall data generation cost. The per-base cost of Illumina-based NGS sequencing is falling thanks to improvements relating to flow-cell cluster density and the generation of longer reads (although for most aDNA this latter point is rarely beneficial). As such, today a modern human 3 gigabase haploid genome can be sequenced to ×30 coverage for as little as USD$1000, and possibly even as low as USD$100. Palaeogenomicists, however, are not so fortunate, given that much (if not, in many cases, the majority) of the DNA in most ancient samples is derived from exogenous contaminants such as microbes. While some methodological improvements such as optimized choice of tissue sources, extraction methods, and various forms of enrichment help improve the endogenous DNA content, costs can still be many fold that for modern DNA data. Thus while attractive to many, the application of palaeogenomics has been largely restricted to the most well-funded research teams and spectacular research questions.

4While Illumina has dominated the palaeogenomic sequencing market, in 2016, a new platform emerged that may offer considerable potential to the field—the combinatorial probe-anchor synthesis–based BGISEQ-500. The underlying technology combines DNA nanoball (DNB) nanoarrays with polymerase-based stepwise sequencing, and its use has recently been validated as comparative in performance to the Illumina platforms when sequencing small noncoding RNAs as well as resequencing modern human DNA. The BGISEQ-500 has several features that suggest it will be it attractive to aDNA users. First, its sequencing read-length capacity (currently up to either single read [SR] or paired end [PE] 100 bp) falls within lengths that are acceptable to most palaeogenomicists. Second, its high throughput—a single 2 channel flow cell can produce at least 500 million single-end reads per channel (thus up to at least 2 billion PE reads per flow cell) in only a few days. Third, at least the initial stages of the library construction method underlying the BGISEQ-500 are sufficiently close to the methods currently used for Illumina palaeogenomic sequencing, and thus can be easily modified based upon some of the abovementioned previous aDNA-related developments if needed. To fully explore this platform's potential for aDNA, we therefore undertook a direct performance comparison against Illumina technology by building libraries and sequencing 8 historic and ancient DNA extracts. To both keep the underlying variables as similar as possible and to exploit a recent (Illumina-based) methodological development that (i) simplifies library construction and minimizes hands-on time and economic cost and (ii) performs at least as well as the Meyer and Kircher blunt end ligation method that many palaeogenomicists favour, we did not use the original BGISEQ-500 library method, but rather developed a new protocol based on our recently developed blunt end single tube (BEST) method. We subsequently undertook a range of bioinformatic analyses aimed at exploring whether the resulting sequence datasets (i.e., Illumina vs BGISEQ-500) exhibited significant differences with regards to a number of parameters that are currently believed important for aDNA studies.

5DNA was extracted from 8 historic and ancient canid samples, chosen so as to represent a range of materials that are currently interesting to the palaeogenomics community (Table 1)—in particular with regards to the fragment sizes of the surviving DNA and the range of endogenous DNA content within them. Two of the samples are preserved hides of wolves (Canis lupus), between 91 and 148 years old, that are believed to contain relatively pure (free of enzymatic inhibitors), although heavily fragmented, DNA (a presumed side effect of the tanning process). The remaining samples are naturally preserved wolf, dog (Canis familiaris), or undetermined large canid remains dated between roughly 600 and 14 000 years old.

6Sample CN 214 was acquired by and registered in the collections of the Natural History Museum of Denmark (NHMD) in 1869. According to museum records, the specimen was shot in Uummannaq, West Greenland, prior to 1869. CN 1921 is a wolf that was shot in Rosenvinge Bugt, East Greenland, in 1925, and then subsequently placed in the NHMD collections. MGUH VP 3332 is a find belonging to the Greenland National Museum (GNM), specifically a bone sample found on the surface 2 m above sea level in 1979 in Vølvedal Peary Land, North Greenland. The specimen has been directly dated to 6785 ± 100 14C years BP (Ua-1346, calibrated age: 7620 years BP). NKA 1950 × 2906 is a tooth sample excavated at the Greenlandic Norse GUS site (Gården Under Sandet/The Farm Beneath the Sand) and is placed in the GNM collections. The site was settled by the Greenlandic Norse and inhabited between ca. 1000 to 600 years BP. ZMK 350/1982 was excavated from the Saqqaq cultural Paleo-Eskimo site Qajâ and is placed in the GNM collections. Although the site in general has been dated to between 3600 and 2700 years BP, this particular sample is from the earliest occupation layers. Last, samples FRC, L, and M1 are tissue samples from an extremely well preserved mummified canid found in the permafrost near the village Tumat in the Sakha Republic, Siberia, Russia. The specimen belongs to the collections of the Mammoth Museum in Yakutsk (Russia) and has been directly dated to 12 223 ± 34 14C years BP (ETH-73412, calibrated age ranging from 12 297 BC to 12 047 years BC, with 95.4% likelihood calibrated to ca. 14 122 years BP); calibration was made using OxCal v. 4.2.4..

7Following DNA extraction, 2 aliquots of each extract were constructed in the Illumina and BGISEQ-500 libraries, respectively, using identical amounts of starting material (16.3 μl, ∼5–50 ng DNA input sample dependent), and then subsequently sequenced to enable bioinformatic comparisons on the data.

8We initially generated between 1.35 × 107 and 5.94 × 107 reads per Illumina library, and 2.32 × 107 - 3.39 × 108 reads per BGISEQ-500 library (Table 2; Supplemental Table S1). The dataset supporting the results of this article is available in the ERDA and GigaDB repositories (see Availability of supporting data). Following normalization of the data for read length and depth (Table 2), we found no statistically significant difference between the 2 datasets with regards to the % endogenous nuclear DNA and average length of endogenous DNA, several of the most important parameters for palaeogenomicists, given their fundamental role in affecting the overall financial cost of a study (Table 3). In contrast, there was a statistically significant platform-dependent difference in the % reads mapping to the mitochondrial genome, with fewer reads mapping in the BGISEQ-500 libraries. However, closer inspection of the data indicates that the total numbers of mtDNA reads are extremely low for 3 of the samples (Supplemental Table S1) and that for the remainder the numbers are extremely similar. As such, we do not believe there to be much significance behind this observation.

9With regards to sequence accuracy, although double-strand (δD) sequence damage rates as estimated using MapDamage2.0 showed no statistically significant difference, a small, yet statistically significant difference was observed for δS, the single-strand damage parameter (lower rate for BGISEQ-500) (Tables 2 and 3). Furthermore, we also observed a small, yet significant difference in the background rate of differences from the reference genome (MapDamage2.0 θ), with slightly higher values observed in the BGISEQ-500 platform (Tables 2 and 3). We hypothesize that both differences may be explained by the fact that, while the initial steps of the library build methodologies are similar, a greater number of polymerase chain reaction (PCR) cycles was used to amplify the Illumina libraries (Supplemental Table S3). This had a clear effect on overall library complexity as while there was no statistically significant difference with regards to library clonality levels or the % reads retained after initial filtering (Tables 2 and 3), when we used preseq to extrapolate on the library complexity, we observed that in all but 1 case, the BGISEQ-500 platform provided richer libraries (Fig. 1). Alternatively, we hypothesize that an alternative explanation for the observed differences in δS and θ might relate to the relatively low genome coverage that we have for each sample. As such, each sample was sequenced over different parts of the genome, which in turn may lead to small biases in the error profiles. Ultimately, however, we feel that full resolution of the differences will require the generation of extensive extra data, and thus more will be learnt in future studies that use the BGISEQ-500.

10We subsequently explored 2 further parameters that relate to whether there are method-specific biases with regards to which part of the genome is sequenced: k-mer frequency and GC content. k-mer content was consistent between methods for most of the samples, each sample pair clustered together. However, samples P83 and 1921 were exceptions to this pattern, with each method yielding slightly different k-mer distributions (Fig. 2). We note that the k-mer content of sample P83 is very similar to sample M1, which makes accurate clustering more challenging. The differences for sample 1921 are more difficult to explain, however, although 1 obvious point is that this is the sole BGISEQ-500 library to exhibit lower complexity than its Illumina pair, although it is not clear if/how this may affect the results.

11GC content was also largely consistent between methods. At a global level, we found no statistically significant difference in the average GC content (Tables 2 and 3), and in more refined analyses, we observed the fragment count for the same windows to be well correlated between BGISEQ-500- and Illumina-derived reads, both of which are correlated with GC content (Figs 3 and 4). We find high genome-wide coefficients of determination for samples 1921,214,L and M1, while these values are lower for samples FRC, P83, and P84 (see Table 4; the sample P79 was excluded from this analysis because of insufficient data). We believe these differences are most likely attributable to the overall endogenous DNA quality in the samples rather than the platforms’ technical performance as there is a trend of samples with lower endogenous DNA content having poorer correlations.

12Our final analysis explored copy number variation (CNV) levels although, as mentioned above, the low genomic coverage of the data makes CNV analyses challenging. Nevertheless, the r2 values for the comparisons that pass our quality control range from 0.35 to 0.96 (Table 5). Furthermore, the observation of particular DNA extractions with excellent concordance values despite the nature of our experiment make it tempting to speculate that indeed both technologies are viable for high-quality CNV calls. For example, using 36-mers and accounting for all possible placements of a 36-mer, the sample M1 has a coverage of above ×1 on both platforms. Ultimately, however, it is not possible to discern from the present data whether the observed variation in CNV calls in the samples is due to differences in the sequencing platforms or the nature of the libraries; thus these results should be taken as preliminary, pending future validation.

13Our study represents the first exploration of the applicability of the BGISEQ-500 as an alternative sequencing platform to the Illumina series for palaeogenomic sequencing, and in doing so we present a library build protocol to generate such data. Although our study is based on only 8 specimens, given that their ranges of endogenous DNA content (<1–75%) and normalized average endogenous DNA sequence lengths (ca 42–76 bp) are typical of many other ancient samples, we anticipate that our results will be indicative of the platform on such material in general. Overall, the results are extremely promising—the BGISEQ-500’s performance is comparable over all parameters tested, with the exception of the very slightly elevated error rate observed (although in contrast we observe higher library complexity and lower δS; thus overall we feel this will not represent a major concern to palaeogenomic studies). We do caution, however, that due to the small size of the dataset (both sample numbers and sequencing depth), at this point we are not able to offer any comment as to how this overall evidence of consistency may translate into downstream analyses involving whole genome summary statistics. Thus we strongly advocate that those who may be interested in using the BGISEQ-500 platform in genomic population data explore this point further. Furthermore, as additional datasets are generated, we look forward to the results of analyses that might wish to compare the relative performance of different sequence alignment and variant calling software on such data. Ultimately, however, we anticipate that our findings will stimulate considerable interest in the BGISEQ-500 platform by palaeogenomic research teams attempting to reconstruct ancient genomes and transcriptomes, and we look forward to future exploration of its potential across a wider range of ancient substrates.

14DNA was extracted using 1 of 3 different methods (designated A, B, C) (Table 1), as deemed appropriate for the choice of tissue. Methods A and C involved digestion in a proteinase K–containing buffer, following Gilbert et al., while method B involved digestion in a proteinase K–urea buffer, following Ersmark et al.. All samples were predigested at 56°C for 1 hour, after which the buffer was changed and then a second 12-hour digest was performed. Digests from method A used organic solvents (phenol: chloroform) and Qiagen MinElute columns (Qiagen, Hilden, Germany), following Carøe et al.. Digests from methods B and C were centrifuged at 6000 ×G for 1 minute, after which 500 μl supernatant was mixed 1:8 with a binding buffer as detailed in Allentoft et al., then centrifuged through Monarch DNA Cleanup Columns (New England Biolabs, MA, USA). DNA bound to the columns was washed with 800 μl buffer PE (Qiagen), then eluted using 2 washes in 17 μl buffer EB (Qiagen)—each with an incubation for 5 minutes at 37°C. Prior to library construction, small aliquots of each extract were analysed on an Agilent 2200 TapeStation HS chip (Agilent Technologies, Palo Alto, CA, USA) for fragment size estimation and molar concentration.

15Two aliquots of each extract were constructed in the Illumina and BGISEQ-500 libraries, respectively, using identical amounts of starting material (16.3 μl, ∼5–50 ng DNA input sample dependent) (Supplemental Table S2). Library blanks and index PCR blanks were also included to evaluate the potential contaminations during the library building process. Illumina libraries were constructed using a method based on the recently published single tube “BEST” protocol, largely following Carøe et al., although with some modifications (Supplemental File F1). To enable direct comparison of the sequencing methods, we chose not to use the conventional BGISEQ-500 library construction protocol. Rather, given the similarities between the initial processes of library construction between both methods (DNA end repair and adapter ligation), we modified the BEST protocol to be BGISEQ-500 compatible. Specifically, the standard Illumina-compatible adapters were replaced with BGISEQ-500-compatible adapters AD1 and AD2 (Supplemental Table S4). These adapters were synthesized as 2 pairs of complementary oligonucleotides (AD1_Long and AD1_Short, and AD2_Long and AD2_Short, respectively), then prepared into the final adapters, AD1 and AD2. Specifically, adapters were first diluted to 500 μM with ×1 TE buffer (10 mM Tris-HCl, 1 mM EDTA, pH 8.0, Sigma-Aldrich). Subsequently, an equimolar concentration of each pair of long and short adapters was mixed together and hybridized through incubation at 95°C for 1 minute, followed by a decrease in temperature with 0.1°C/s from 95°C to 12°C. After hybridization, adapters AD1 and AD2 were mixed and diluted at a concentration of 10 μM prior to their use in the library construction. We additionally designed BGISEQ-500-compatible library amplification primers for use in the library amplification steps that included 8 alternate sequencing indices in the reverse primers (Supplemental Table S4).

16Following the final Bst fill-in step during library build, all libraries were mixed with 1:5 volume of PB binding buffer (Qiagen) and purified using Monarch® DNA Cleanup Columns, then washed with 750 μl buffer PE (Qiagen) and eluted in 40 μl buffer EB (Qiagen) after a 5-minute incubation at 37°C.

17Quantitative real-time PCR (qPCR) was used to estimate the required number of cycles during library index PCR. Each qPCR was performed in a 20 μl reaction volume using 1:20 dilution of purified library template, 0.2 mM dNTPs (Invitrogen), 0.04 U/μl AmpliTaq Gold DNA polymerase (Applied Biosystems, Foster City, CA, USA), 2.5 mM MgCl2 (Applied Biosystems), 1X GeneAmp® 10X PCR Buffer II (Applied Biosystems), 1 μl SYBR Green (Invitrogen, Carlsbad, CA, USA), 0.2 μM forward and reverse primers mixture (IS7 and IS8 primers), and 13.48 μl AccuGene molecular biology water (Lonza). qPCR cycling conditions were 95°C for 10 minutes, followed by 40 cycles of 95°C for 30 seconds, 60°C for 60 seconds, and 72°C for 60 seconds using the MX3005 qPCR machine (Agilent).

18Post-qPCR, library index amplifications were performed in 100 μl PCR reactions that contained 20 μl of purified library, 0.2 mM dNTPs (Invitrogen), 0.1 U/μl AmpliTaq Gold DNA polymerase (Applied Biosystems), 2.5 mM MgCl2 (Applied Biosystems), 1X GeneAmp® 10X PCR Buffer II (Applied Biosystems), 0.4 mg/ml BSA (New England Biolabs Inc), 0.2 μM of each forward (Illumina InPE 1.0 forward) and custom made reverse primers, and 51.2 μl AccuGene molecular biology water (Lonza, Basel, CH). PCR cycling conditions were: initial denaturation at 95°C for 12 minutes followed by 13 to 21 cycles of 95°C for 30 seconds, 60°C for 30 seconds, and 72°C for 40 seconds, and a final elongation step at 72°C for 5 minutes. Post-PCR, libraries were purified with QiaQuick columns (Qiagen) and eluted with 30 μl buffer EB (Qiagen) after an incubation for 10 minutes at 37°C. Small aliquots of this purified product were used for quantification and fragment size estimation using the High-Sensitivity DNA Assay for the Bioanalyzer 2100 (Agilent). Subsequently, there was a final purification using the AMpure XP system (Agentcourt, Beckman Counter, Indianapolis, IN, USA) with ×1.8 beads:library ratio, in order to remove any persisting primer dimers or other molecules with a fragment size of <100 bp. Last, libraries were pooled in equimolar concentrations (∼9.4 nM) and sequenced on the Illumina HiSeq platform in 80 bp single read mode by The Danish National High-Throughput DNA Sequencing Centre.

19Initial processing steps for the purified BGISEQ-500 libraries were largely similar to those used in the Illumina libraries, although with the following modifications. First, the libraries were qPCR-quantified using the CommonprimerBGI forward primer and 1 of the indexed reverse primers (Supplemental Table S4). Second, subsequent index PCR amplifications used 8 to 15 cycles (Supplemental Table S3) with CommonprimerBGI forward primer and the indexed reverse primers (Supplemental Table S4). Third, because several of the BGISEQ-500 libraries exhibited residual adapter dimers after the initial purification post-index PCR, each purified BGISEQ-500 library was split to 2 aliquots (∼12.5 μl each), and 1 of each aliquot was subject to an extra purification to remove any residual primer dimers (Supplemental Table S2). Each of these aliquots was sequenced independently. We note that several of the extrapurified libraries showed small improvements with regards to overall adapter dimer content in the generated sequence (Supplemental Table S5), and our initial impression is that this extra purification step may be worth undertaking if high levels of adapter dimers are found post-index PCR.

20All amplified libraries were subsequently sent to BGI for circularization and sequencing on the BGISEQ-500 platform. For circularization, PCR products with different barcodes were pooled together at equimolar concentration to yield a final amount of 80 ng. Pools contained both the samples relevant to this study as well as those from other projects (Supplemental Table S6). Each pool was subsequently heat-denatured, and the single-strand DNA was mixed with MGIEasyTM DNA Library Prep Kit V1 (PN:85–05533-00, BGI, Shenzhen, China), containing 5 μl splint oligo, 6 μl splint Buffer, 0.6 μl ligation Enhancer, and 0.2 μl ligation (Enzyme and NF water) to form a 60 μl reaction system, which was subsequently incubated at 37°C for 30 minutes. Last, 20 μl of each single-circle-library pool was used as input to prepare the DNB. Each pool was then sequenced on 1 lane, using 100SR chemistry with BGISEQ-500RS High-throughput sequencing kit (PN: 85–05238-01, BGI). Postsequencing, the data were automatically demultiplexed by index.

21The raw reads obtained from the HiSeq 2500 and BGISEQ-500 were analysed using FastQC (FastQC, RRID:SCR_014583) to compute the quality metrics of the reads, such as base sequence qualities, base sequence content, %GC, and sequence composition. With the exception of the analysis on the standard vs extrapurified BGISEQ-500 libraries (Supplemental Table S5), both BGISEQ-500 libraries from each extract were treated as a single dataset. We also compared the quality metrics of the reads from the same samples across the 2 platforms to ensure that the sequencing platform did not have a large impact on the quality metrics of the reads.

22Once the read qualities were verified using FastQC (FastQC, RRID:SCR_014583), we used the PALEOMIX pipeline (PALEOMIX, RRID:SCR_015057) to trim the adapter sequences, trim Ns and low-quality bases from the ends of reads, estimate ancient DNA damage, and finally map the trimmed reads to the reference genome. The individual steps of the pipeline are detailed below. We highlight that the values presented in Table 2 are normalized to account for sequencing read depth and length, while Supplemental Table S1 contains both the original and the normalized values.

23The first step of the initial processing of the reads involved trimming the adapter sequences from the ends of the reads. Since the samples consist of degraded DNA, many of the sequenced reads contain the platform-specific adapters at the 3΄ end of the reads. AdapterRemoval (v. 2.1.3) was used to trim the adapter sequences from the ends of the reads using the default mismatch rate of 1/3. In addition, bases with a quality score of less than 2 and unidentified bases (Ns) at the ends of reads were trimmed. Finally, only reads that were longer than 25 bases were retained for downstream analyses.

24The trimmed reads were mapped to the wolf reference genome using the mem algorithm in bwa (BWA, RRID:SCR_010910; v. 0.7.10), using the default settings for the mapping algorithm. The mapped reads were subsequently processed using the GATK (v. 3.3.0) indel realigner (GATK, RRID:SCR_001876) to fix the alignment issues arising from the presence of short indels at the beginnings and ends of reads. Since there are no catalogs of indel variations in the species included in this study, the realignment step was done using a set of indels within each sample. After the indel realignment step, the PCR duplicates were removed from the alignments using the MarkDuplicates program from Picard tools (Picard, RRID:SCR_006525; v. 1.128).

25The DNA damage patterns and parameters were estimated using mapDamage (v.2.0.6; mapDamage, RRID:SCR_001240) using a subsample of 100 000 reads from the set of mapped reads. The 3 main parameters estimated using mapDamage were θ, δD, and δS. δD and δS estimate the probability of cytosine deamination (driven by hydrolytic DNA damage) in a double- and single-stranded context, while θ estimates the background rate of difference between the reference and sample after accounting for DNA damage. Using these estimated parameters, the base qualities of putatively damaged bases were recalibrated to a lower score. The program was also used to compute the relative abundance of C→T changes at the 3΄ ends and A→G changes at the 5΄ ends of the reads and compare them across the 2 platforms.

26The clonality of each library was computed from the reads that were identified by the MarkDuplicates program during the duplicate identification and removal step. The clonality was computed as the ratio of the number of reads retained after duplicate removal and the number of reads retained after the adapter removal and trimming step. The endogenous content of the library was computed as the ratio of the number of reads mapping uniquely to the reference genome and the number of reads retained after adapter removal. Note that this is 1 possible definition of the endogenous content, here defined as the proportion of usable reads obtained from a library, and the numbers given in Table 2 and Supplemental Table S1 will allow you to compute the values for other definitions of endogenous content.

27The complexity of each library was estimated and extrapolated using the library complexity extrapolation model in the program preseq, which uses a nonparametric Bayesian Poisson model to estimate the gain in number of unique fragments when the library is sequenced deeper. Instead of using the aligned reads to estimate the library complexity, we used the counts of the number of duplicates in the bams generated by paleomix as input to preseq. The library complexity was estimated up to a maximum of a total of 10 billion reads sequenced per library.

28Since the draft de novo wolf genome does not contain information on scaffolds that are annotated as belonging to the mitochondria, we could not identify reads that mapped to the mitochondrial genome using the initial set of mapped reads. To overcome this problem, we downloaded a complete mitochondrial genome from NCBI (GenBank Accession: AM711902) and mapped the adapter trimmed reads to this complete mitochondrial genome. The same steps, including indel realignment and DNA damage–related recalibration of quality scores, were performed for the reads aligned to the mitochondria.

29To compare the sequence content of the reads obtained from the 2 sequencing platforms, we computed the k-mer frequencies in the reads from the same sample using the 2 technologies. Since the raw reads are enriched in adapter sequences and do not accurately reflect the sequence content of the underlying endogenous DNA molecules in the library, we restricted the k-mer analysis to reads that mapped to the genome after going through both adapter trimming and duplicate read removal. For each library, we sampled 100 000 reads from the reads mapped to the reference genome using samtools (v. 1.2; SAMTOOLS, RRID:SCR_002105) and seqtk (v. 1.0). From these subsampled reads, we computed the 6-mer frequencies using jellyfish (Jellyfish, RRID:SCR_005491).

30The relationship between read abundance in a given genomic region and its GC content is well known and characterized for the Illumina platform. For methods that depend upon depth of coverage or fragment count, such as measuring absolute copy number or expression levels, this bias needs to be taken into consideration and corrected for; otherwise, its magnitude might confound the signal in question. We therefore compared the GC content of the mapped endogenous DNA for the 2 platforms in several ways. First, the basic GC percentage was calculated from all endogenous reads. Second, we partitioned the reference genome into bins of 100 Kbps, with an offset of 10 Kbps, and calculated the GC percentage of each bin. We then mapped all datasets onto the reference and counted the number of mapped fragments in each bin. To account for differences in sequencing depth, we randomly subsampled mapped reads from the platform with the higher coverage to an equal amount of mapped bases of the platform with lower coverage, and then normalized the number of mappings by the median number of mappings for each extract.

31Fluctuations in depth of sequencing coverage can be used to generate personal genome-wide copy number maps of an individual as read depth is known to strongly correlate with copy number for several platforms. We sought to assess whether the same techniques might be applied to data generated on the BGISEQ-500. To this end, we generated individual genome-wide CN maps of all extracts and both platforms in varying window sizes, from 1 Kbp to 1 Mbp, to account for fluctuation in coverage, and checked concordance between them. It is worth noting that using ancient DNA libraries poses a particular challenge to this assessment as some inherent characteristics of this type of data (such as unequal degradation, fragmentation, or clonality during library preparation) make it difficult to pinpoint the source of variability between 2 call sets for a given extract, given a lack of concordance. Specifically, low effective coverage and poor DNA quality make high-resolution maps not feasible for many of the libraries used in this part of the project.

32We masked out any repeats in the reference assembly, as identified by both repeat masker (RepeatMasker, RRID:SCR_012954) and tandem repeat finder. Additionally, to identify repeats that have been potentially missed by the aforementioned algorithms, we chopped up the masked assembly into 36-mers with an offset of 5 bp. These were then mapped back onto the assembly using GEM (GEM, RRID:SCR_005339) with a maximum divergence set to 95% and retaining all possible mappings. All 36-mers with more than 20 placements along the genome were additionally masked out. We then generated nonoverlapping 36-mers from the production reads and mapped them onto the extensively masked reference assembly using GEM, allowing for a maximum divergence of 95% and retaining all possible placements. To call absolute copy number, the reference was portioned in windows of 1, 5, 10, 50, 100, and 1000 Kbps of nonoverlapping, nonrepetitive sequences with mrCanavar (mrCaNaVaR, RRID:SCR_003135), meaning that the genomic coordinates of the windows may span more than the window size if repeats are present within it. Importantly, as reads may not properly map at the boundaries of maskings, we introduced an additional padding of 36 bp. We then iteratively excluded all windows that represent outliers with respect to a normal distribution to identify a set of “control regions.” After correcting for GC content, the median depth of coverage in these control regions was used to normalize all windows and thus assign an absolute copy number to them. The concordance was calculated as the coefficient of determination of a linear model over corresponding to windows of the same extract between the 2 platforms. Additional quality control involved visually inspecting the normalized read depth distribution of the aforementioned control regions. In a good sample, this should be a symmetrical, bell-shaped curve centered at 2. We visually inspected all distributions and classified them as good, neutral, or bad, based on shape and symmetry. In addition to the aforementioned windows (called copy windows [CWs]), we also calculated normalized read depths the same size as CWs in terms of nonrepetitive sequence, with a fixed offset of the window size, but including repetitive sequence (called short windows [SWs]), and windows 5 times the size of a copy window (called long windows [LWs]), with an offset of 5 times the size of a copy window, but including repetitive sequence. As an additional quality control, the ratios of read depth of SW/CW should be around 1, and the ratios of read depth of LW/CW around 5, given proper sampling of the genome.

33Supplemental File F1: Improvements to original BEST library building protocol (see additional file).

34Supplemental Table S1: Full sequence data information (see additional file).

35Supplemental Table S2: Sequence library identifiers.

36Supplemental Table S3: The number of index PCR cycles used in each sample.

37Supplemental Table S4: The sequences of BGISEQ-500 adapters and index primers used in this study.

38Supplemental Table S5: Adapter dimer content of initial and extrapurified BGISEQ-500 libraries.

39Supplemental Table S6: Library pooling for BGISEQ-500 library circularization reactions.

40δD: MapDamage 2.0 double-strand DNA damage rate; δS: MapDamage 2.0 single-strand DNA damage rate; θ: MapDamage 2.0 DNA damage-corrected error rate; aDNA: ancient DNA; BEST: blunt end single tube; CN: copy number; CNV: copy number variation; CW: copy window; DNB: DNA nanoball; GNM: Greenland National Museum; LW: long window; NFC: normalized fragment count; NGS: next-generation sequencing; NHMD: Natural History Museum of Denmark; PE: paired end; SR: single read; SW: short window; YBP: years before present.