With regards to sequence accuracy, although double-strand (δD) sequence damage rates as estimated using MapDamage2.0 [ 17 ] 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 [ 85 ] 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.

Figure 1:

Open in new tab Download slide

Library complexity estimated as the number of unique reads as a function of the total number of reads sequenced. These numbers are estimated and extrapolated using the program preseq [ 84 ]. The total number of reads sequenced for each library can be found in Table 2 and Supplemental Table S1. The solid lines are the estimates for the libraries sequenced on the Illumina HiSeq 2500 platform, while the dotted lines are the estimates for the libraries sequenced on the BGISEQ-500. Each of the 8 samples is represented by a different colour.

We 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.

Figure 2:

Open in new tab Download slide

Heatmap of k-mer counts across libraries. Libraries (columns) were hierarchically clustered based on Pearson correlation. Proportion of each of the 4096 6-mer (rows) are depicted using colours.

GC 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.

Figure 3:

Open in new tab Download slide

Top: median normalized fragment count (NFC) per 100 Kb windows, with 10 Kb offset for the sample 214 along scaffold_0. The solid line shows Illumina data, and the dotted line shows BGISEQ-500 data. Bottom: percentage GC calculated over the same the same windows as in the upper panel.

Figure 4:

Open in new tab Download slide

Median NFC of Illumina vs BGISEQ-500 for all samples in windows of 100 Kb, with an offset of 1 0Kb along scaffold_0. The color of each point corresponds to the windows’ GC content. For the high-quality samples (1921, 214, FRC, M1), a very good correlation of NFC between the 2 platforms can be observed. The fragment count seems to be correlated with GC content.

Our 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 r ² 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.

Potential Implications

Our 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.

Methods

DNA extraction

DNA 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. [ 62 ], while method B involved digestion in a proteinase K–urea buffer, following Ersmark et al. [ 86 ]. 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. [ 13 ]. 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. [ 42 ], 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.

Library construction

Two 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. [ 13 ], 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).

Following 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.

Illumina library PCR amplification and sequencing

Quantitative 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 MgCl ₂ (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 [ 12 ]), 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).

Post-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 MgCl ₂ (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.

BGISEQ-500 library PCR amplification

Initial 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.

BGISEQ-500 library circularization and sequencing

All 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.

Data analyses

The raw reads obtained from the HiSeq 2500 and BGISEQ-500 were analysed using FastQC (FastQC, RRID:SCR_014583 ) [ 87 ] 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.

Once the read qualities were verified using FastQC (FastQC, RRID:SCR_014583 ), we used the PALEOMIX pipeline (PALEOMIX, RRID:SCR_015057 ) [ 18 ] 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.

Adapter removal and trimming

The 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) [ 88 ] 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.

Mapping, indel realignment, and duplicate removal

The trimmed reads were mapped to the wolf reference genome [ 89 ] 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 ) [ 90 , 91 ] 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) [ 92 ].

DNA damage patterns

The DNA damage patterns and parameters were estimated using mapDamage (v.2.0.6; mapDamage, RRID:SCR_001240 ) [ 17 ] 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.

Clonality, endogenous DNA content, and library complexity estimation

The 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.

The complexity of each library was estimated and extrapolated using the library complexity extrapolation model in the program preseq [ 85 ], 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.

Mapping to the wolf mitochondrial genome

Since 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) [ 93 ] 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.

K-mer frequency

To 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 ) [ 94 , 95 ] and seqtk (v. 1.0) [ 96 ]. From these subsampled reads, we computed the 6-mer frequencies using jellyfish (Jellyfish, RRID:SCR_005491 ) [ 97 ].

Relative abundance vs GC content

The relationship between read abundance in a given genomic region and its GC content is well known and characterized for the Illumina platform [ 98 ]. 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.

CNV on low-coverage data

Fluctuations 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 [ 99 ]. 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.

We masked out any repeats in the reference assembly, as identified by both repeat masker (RepeatMasker, RRID:SCR_012954 ) [ 100 ] and tandem repeat finder [ 101 ]. 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 ) [ 102 ] 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 ) [ 99 ], 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.

Additional files

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

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

Supplemental Table S2: Sequence library identifiers.

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

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

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

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

Abbreviations

δ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.

Acknowledgements

The authors would like to acknowledge the assistance of the Danish National High-Throughput Sequencing Centre for assistance in Illumina data generation. The authors would like to acknowledge the ERC Consolidator Grant (681396 – Extinction Genomics), the Marie Skłodowska-Curie Actions (H2020-MSCA-ETN-643063 “Microwine”), Danish Council for Independent Research (4005–00107 Wine-ometrics), the Qimmeq project, funded by The Velux Foundations and Aage og Johanne Louis-Hansens Fond, China National GeneBank, and BGI-Shenzhen China for funding. We also gratefully acknowledge the Danish National Supercomputer for Life Sciences–Computerome (computerome.dtu.dk) for the computational resources to perform the sequence analyses. L.F.K.K is supported by an FPI fellowship associated with BFU2014–55090-P (FEDER), T.M.B. is supported by MINECO BFU2014–55090-P (FEDER) and BFU2015–6215-ERC, a U01 MH106874 grant, and the Secretaria d’Universitats i Recerca del Departament d’Economia i Coneixement de la Generalitat de Catalunya.

Availability of supporting data

Raw sequencing data are available from the SRA [PRJEB21089]. All other datasets supporting the results of this article are available in the ERDA repository [ 103 ] and GigaScience repository, GigaDB [ 104 ]. DNA extraction and library construction protocols presented here are also archived in protocols.io [ 105 ].

Competing interests

The authors declare that Hui Jiang, Chunyu Geng, Guojie Zhang, Wenwei Zhang, Shujin Fu, and Shanlin Liu are employees of BGI.

Author contributions

M.T.P.G., G.Z., and H.J. conceived the study with critical input from S.L., C.C., and S.S.T.M. C.C. adapted the BGISEQ-500 library construction method for aDNA. S.S.T.M. prepared the aDNA libraries. M.H.S.S. extracted the aDNA. C.G., W.Z., and S.F. performed BGISEQ-500 library circularization, ssDNA synthesis, and BGISEQ-500 sequencing. S.G., F.G.V., L.F.K.K., and T.M.B. analysed the data with assistance from S.L., T.S.P., and B.P. M.T.P.G. drafted the manuscript, with input from all authors.

References