GoKawiilDNA extraction and library preparation
All laboratory work was conducted in dedicated ancient DNA clean room facilities, at the Max Planck Institute for Evolutionary Anthropology, the Royal Belgian Institute of Natural Sciences or the University of Tübingen. Bone powder (1–57.1 mg; Supplementary Table 1.1) was obtained from 35 skeletal remains using a sterile dentistry drill. More than one-third of the powdered samples underwent a pretreatment with either 0.5% hypochlorite solution or a phosphate buffer, in an attempt to reduce present-day contamination13, before DNA extraction using a silica-based protocol optimized for recovery of short, degraded DNA fragments14. Single-stranded sequencing libraries were prepared using between 5 μl and 30 μl of the extracts15, for most libraries using automated liquid handling (https://doi.org/10.17504/protocols.io.kqdg32bdpv25/v1)61. A subset of the earliest libraries were prepared using a modified single-stranded DNA protocol with partial uracil-specific excision reagent (USER) treatment, reducing uracil content while preserving terminal deamination signals62. To enable multiplex sequencing, each library was indexed with a pair of unique 7 or 8 base pair (bp) barcodes63,64. We quantified library yields by quantitative PCR or digital droplet PCR15 (https://doi.org/10.17504/protocols.io.bp2l6xwd5lqe/v2)65, and quantified library preparation efficiency using spiked-in control oligonucleotides66. Negative controls were included in all extraction, library preparation and sequencing steps to monitor for contamination.
Shotgun sequencing
To screen for endogenous DNA preservation, libraries were shallowly sequenced on Illumina platforms (Supplementary Information section 1.2). Sequencing data processing included adapter trimming and read merging using leeHom67 (v.1.2.18, https://bioinf.eva.mpg.de/leehom/), followed by mapping to the human reference genome (hg19/GRCh37) with BWA68 (v.0.5.10-evan.9-1-g44db244, https://github.com/mpieva/network-aware-bwa) using parameters optimized for ancient DNA69. PCR duplicates were removed using bam-rmdup (https://github.com/mpieva/biohazard-tools/) and reads were filtered for a minimum length of 35 bp and a mapping quality of 25. DNA authenticity was evaluated by determining the frequency of deamination-induced C-to-T substitutions at the ends of reads, requiring more than 10% substitutions at terminal bases as evidence for the presence of authentic ancient DNA. Contamination levels were assessed by comparing deamination rates before and after conditioning on damage at one end70 (Supplementary Table 1.1).
Mitochondrial capture
Libraries were enriched for mtDNA through two rounds of in-solution hybridization capture (Supplementary Information section 10). We used probes tiling the revised Cambridge reference sequence of the human mtDNA genome (rCRS71) for all samples except A57813, which was captured with an array containing mitochondrial genomes from 241 mammalian species in addition to the rCRS reference72. After sequencing on Illumina MiSeq or HiSeq platforms, reads were processed with leeHom67 and mapped using BWA39 (v.0.5.10) with ancient DNA parameters69 to both rCRS and a circularized Vindija 33.16 mitochondrial reference. After duplicate removal and quality filtering (mapQ > 25, length > 35 bp), authenticity was confirmed by requiring at least 10% C-to-T damage. Contamination from modern human DNA was conservatively estimated by counting fragments with modern human-specific alleles at lineage-diagnostic positions in the rCRS alignment73. We compared the effect of using all sequences versus only those with putative terminal deamination, as well as two different thresholds for contamination (10% and 25%), and found no significant differences in the resulting mitochondrial consensus genomes. Consensus calling was performed with different parameters, requiring a consensus support of either 66% or 80% and a minimum coverage of either 3× or 4×, all while requiring a base quality score of 20 and masking transitions at the three last bases at the ends of the reads to mitigate errors due to damage-induced substitutions.
We constructed a dated phylogeny with a multiple alignment, generated using mafft with 1,000 iterations74 (v.7.453), of the new Neanderthal mtDNA genomes, as well as previously published mtDNA genomes from 23 Neanderthals1,2,4,8,11,21,23,34,75,76,77,78,79,80, 55 present-day humans71,81, 10 ancient modern humans16,70,82,83,84,85,86, the Sima de los Huesos hominin70 and 4 Denisovans84,87,88,89 using BEAST2 (v.2.1.3)33. We determined that TrN as the substitution model90 described our data best, and set the proportion of invariant sites of 0.8. We ran this analysis multiple times, using all sequences or only those in the coding portion of the mitochondria, as well as testing different models with different clocks and trees, and using different priors for the tip ages of Couvin G8-0083 and Thorin. Moreover, we also explored the viability of combining non-overlapping radiocarbon dates for GN2, and the effect that using different dates from the different specimens would have on the posterior dates. These analyses are described in more detail in Supplementary Information section 10.
Y chromosome capture
Three out of the four genetically male individuals had enough sequences mapping to the non-recombining portion of the Y chromosome91 for downstream analyses (Supplementary Information section 11). Owing to the low coverage, we did not construct a de novo phylogeny. Instead, the new data of Goyet Q305-1, Fonds-de-Forêt 1 and Trou Magrite 2422-36 were used to place the new individuals into an existing high-quality phylogeny of archaic Y chromosomes1 to determine their affiliations with previously sequenced male Neanderthals.
GN1 high-coverage genome
... continue reading