aDNA lab work

aDNA work was performed in the dedicated facilities of the Institute for Archaeological Sciences in Tübingen, Germany. The HST femur was first irradiated with ultraviolet light on the selected sampling area and then drilled with a dentist drill along the cortical bone. A total of 130 mg of bone powder went into the DNA extraction following an established protocol22. DNA was eluted in 100 μl of TET and 20% of the extract (GX35) was used to build a double-stranded genetic library (GA87)23. The total copies in the resulting library were measured with quantitative PCR (qPCR) (7.53 × 109 copies). They were split into four 100 μl indexing PCR reactions with 10 cycles where an individual index pair (8 bp each) was assigned to create a unique double-indexed library41. The total copies were measured again via qPCR (2.60 × 1011 copies) and the reaction efficiency was calculated by dividing the number of total molecules after indexing by the number of total molecules before indexing PCR. An aliquot of two-fifth of the indexed library was split into two reactions that were amplified for seven cycles each with AccuPrime Pfx DNA polymerase. The PCR products were purified over a single MinElute spin column and the concentration after amplification was quantified to 286 ng μl−1 on an Agilent 2,100 Bioanalyser DNA 1,000 chip. Extraction and library negative controls were carried along the workflow and treated equally.

The amplified library was enriched for mtDNA using modern human baits as reported by Maricic et al.24. This protocol has been previously used to successfully capture complete Neanderthal mtDNA genomes4. Four hundred nanograms of the amplified library were pooled with the same amount of four other libraries for a total of 2,000 ng and captured with 500 ng of mtDNA probes. After purification, the isolated molecules were quantified with qPCR (4.84 × 106 copies) and reamplified for 20 additional cycles as described above. The captured pool as well as the uncaptured GA87 library was quantified with Agilent 2,100 Bioanalyser DNA 1,000 chip, diluted to 10 nM and sequenced with other equimolar libraries on an Illumina HiSeq2500 Rapid run via 2 × 100+8+8 cycles and on an Illumina NextSeq500 run via 2 × 75+8+8 cycles.

Sequence processing and mtDNA consensus reconstruction

Sequenced molecules were converted from bcl to fastq files and reads containing the defining library indexes were binned in an individual folder. The EAGER pipeline was used for all subsequent data processing42. Initially adapter and index sequences were trimmed off. Only merged reads where forward and reverse reads overlapped by at least 10 bp were retained. Shotgun sequences above 30 bp were aligned to the complete human genome (hg19) with Burrow–Wheeler Aligner (parameters −n 0.01 and seeding off) to calculate the percentage of human DNA. Duplicates and reads with mapping quality below 30 were discarded to estimate damage patterns and average fragment length (Supplementary Table 3). From the total of ∼3 Ma paired-end reads sequenced after mtDNA capture, 89.31% were successfully merged and fragments below 30 bp length were further discarded for mapping. The resultant ∼1.3 Ma merged reads were aligned to four reference mtDNA sequences: the RSRS25, the rCRS27, the Neanderthal Feldhofer 1 sequence26 and the RNRS originally proposed in Behar et al.25 and later updated when the more basal Altai mtDNA was published5. The rCRS and Feldhofer 1 references are two derived mtDNA sequences on the modern human and Neanderthal branch, respectively. Instead, RSRS and RNRS represent the MRCA mtDNA for modern humans and Neanderthals, respectively. Reads were mapped using Burrow–Wheeler Aligner43 with identical parameters (−n 5 and seeding off) for all four references, in combination with a tool able to consider the circularity of mtDNA as part of EAGER. The percentage of target DNA was calculated by dividing the total number of input reads by the reads mapping to each mtDNA reference. Duplicates with the same start and end coordinates were removed and the duplication factor was measured by dividing the total number of reads mapping before by the total number of reads mapping after duplicate removal. All fragments with map quality below 30 were removed to estimate the average mtDNA coverage. The resulting molecules were also used to calculate average fragment length and deamination patterns44 (Supplementary Fig. 3). Statistics for each processing step of the four reference sequences are reported in Supplementary Table 2.

Consensus reconstruction was performed in a two-step approach. First, schmutzi28 was used to infer the endogenous sequence. An internal program of the software package (contDeam) was first run to calculate the endogenous deamination rate and a contamination prior. To each nucleotide a base likelihood value was assigned incorporating damage, base quality and mapping quality information in a Bayesian framework28. The endogenous consensus was then determined by endoCaller after the first iteration of the program. No cutoff to the nucleotide posterior probability was selected resulting in base called even in positions covered with only one fragment. This produced a consensus sequence with three unassigned positions.

Second, the four consensus sequences, one from each reference, were visually compared in Geneious 8.1.7 (http://www.geneious.com)45. A multiple genome alignment was produced and each of the 19 inconsistent positions between the four consensus was evaluated. We imported the bam files in Geneious and for each read covering those positions we inspected if they also overlapped neighbouring confidently assigned SNPs (for example, called in all four consensus). Fragments containing such SNPs were considered as endogenous, whereas reads containing the alternative allele were considered as contaminants. In every case, the consensus sequence reconstructed after mapping against the RNRS reference was found to exhibit the endogenous base. This confirms that mapping against a reference sequence that is phylogenetically closer to the consensus sequence increases mapping accuracy (see Supplementary Note 3). Using the same criterion described above we then manually screened the RNRS mapped consensus and edited the following positions according to rCRS coordinates. Two miss-mapped insertions were removed (pos. 247delT and 16184delA), two uncertain positions with low coverage were edited (A189G and A4296N) and two regions covered with only contaminant reads were masked (pos. 203–214Ns and 5486–5508Ns). We additionally masked the known troublesome regions of poly-C (pos. 303–315) and poly-AC (518–524) stretches46. Combining the two approaches resulted in a total of 59 unassigned positions in the final consensus sequence that was used for phylogenetic (Supplementary Fig. 5) and mtDNA diversity analyses. We then generated a more conservative consensus by setting a coverage cutoff to twofold. The resulting mtDNA sequence exhibits 81Ns, but none of the 22 additional unassigned bases overlapped with polymorphic positions within the known Neanderthal mtDNA diversity. Therefore, the tree topology and mtDNA diversity within Neanderthals was not affected.

Contamination with modern human mtDNA

We followed three different approaches to estimate modern human contamination levels in the isolated mtDNA. The first method is implemented in contDeam28 and it relies on deamination patterns. This program works on two assumptions: that modern human DNA contamination presents no damage and that the damage at one end of a molecule is independent of the one at the other end. Reads with deamination at the 5′ end are selected and the deamination rate is measured at the 3′ end and vice versa. The calculated value is supposed to represent the true damage signal of the endogenous mtDNA fragments. A contamination estimate is then computed as the percentage of undamaged reads necessary to shift the damage rate from the endogenous value to the one initially calculated on all fragments. We obtained an estimate ranging from 9.5% to 11.5% (for all four references combined) that could be an underestimation of the real contamination level if the contaminant DNA is also damaged. However, simulations have shown that this effect is marginal if the deamination rate of the endogenous DNA is over 50% at the molecule termini28, like observed for HST mtDNA fragments (Supplementary Fig. 3).

The second method makes use of the probabilistic iterative method implemented in schmutzi28. The program was run with the fallowing parameters: ‘-- notusepredC -uselength’. The contamination estimate is performed with the tool mtCont using sites where the endogenous sequence differs from a non-redundant data set of 256 Eurasian mtDNAs. For this method, we used the mtDNA reads mapping to the rCRS reference according to which base frequencies of the comparative data set are calculated. While contDeam measured contamination rate on a fragment level, mtCont provided an estimate at a base level of 9–11%.

The third method is based on the diagnostic positions where the reconstructed HST consensus differs from present-day worldwide mtDNAs. All polymorphic positions with a frequency above 1% in a data set with 311 worldwide mtDNAs are not considered. We then identified 123 positions where HST has a different base compared to more than 99% of the 311 mtDNAs. Of those, we restricted the analysis to only transversions (positions 2,831, 6,265, 7,105, 9,328, 9,354, 11,457, 13,761, 13,878, 14,457, 14,925, 16,138). Of the total 262 reads covering the 11 positions, 239 reads presented the endogenous base, while 23 reads the contaminant variant. This resulted in a contamination rate of 8.8% (confidence interval 95%, 5.4–12.2%). The last approach30 provides a direct measure for the proportion of contaminant fragments and overlaps with the two previous estimates.

Phylogenetic analyses

To further explore the maternal relationships of the HST mtDNA with other archaic and modern human mtDNAs we compared the phylogenetic placement of the HST consensus sequence with and without D-loop with 17 Neanderthal, 54 modern human47, three Denisovan1,2 and Sima de los Huesos3 mtDNAs, plus a chimpanzee mtDNA (GenBank: X93335.1) to root the tree. Two maximum parsimony trees with the 77 mtDNAs and 1,000 iterations each were built, one including the whole molecule and 97% partial deletion (16,536 positions) and one with the coding region only and 98% partial deletion (15,417 positions) (Fig. 1c and Supplementary Fig. 5). The topology of both trees is consistent with HST diverging from the Neanderthal branch more basally than any other sequence and presenting a short phylogenetic branch length.

We also tested in Modelgenerator48 the same multiple genome alignment with only the coding region, including missing sites but not gaps. The substitution model that best fits the data (AIC1) was GTR with invariant sites and γ-distributed correction for rate heterogeneity. These parameters were selected in MrBayes49, used to build a Bayesian phylogenetic tree (Supplementary Fig. 6). Fifty millions iterations of the Markov chain Monte Carlo (MCMC) were run with 10,000 sampling interval. From the total trees, the first 10% were removed as burn-in and a summarized tree was generated. All major branches show posterior support of 1 and confirmed the maximum parsimony trees topology.

We finally explored the diversity of the Neanderthal HVRI including four additional sequences for which only HVRI was available in GenBank (Valdegoba JQ670672, Scladina DQ464008, Teshik-Tash EU078679, Monti Lessini DQ836132). We aligned them to the HVRI of 17 complete Neanderthal mtDNAs (excluding Denisova 11 because of several unassigned positions in the HVRI), three Denisovan mtDNAs and the rCRS using MUSCLE50. We then built a maximum parsimony phylogeny in MEGA6 (ref. 51) with complete deletion (105 positions) and 1,000 bootstrap iterations (Supplementary Fig. 10).

Mitochondrial DNA diversity

The pairwise nucleotide distance among Neanderthals with (n=18) and without HST (n=17), Denisovan (n=3) and modern human (n=311) mtDNAs was calculated in MEGA6. For this analysis, we used the complete mtDNAs sequences and the number of differences between them was counted with pairwise deletion where all unassigned positions were removed for each sequence pair. We plotted the pairwise nucleotide distance against their frequencies for each of the four data sets (two Neanderthals, Denisovan and modern human) in Fig. 2. We also reported the average distance (and minimal–maximal values) of 311 modern humans, three Denisovans, Sima de los Huesos and 17 Neanderthals to the HST complete mtDNA (Supplementary Table 4). The lowest distance is with Neanderthals, followed by modern humans, Sima de los Huesos and Denisovans in agreement with the phylogenetic assignment. However, the nucleotide distances between HST and other Neanderthals are the largest observed among Neanderthals (89–111 nucleotides). These values are higher than between Denisova 3–Denisova 4 and Denisova 8 and around the uppermost edge among 311 worldwide mtDNAs (Fig. 2).

We further measured the mtDNA diversities of the enlarged Neanderthal mtDNA data set with the Watterson’s estimator, as reported in ref. 2. We first prepared a multiple genome alignment of Neanderthal mtDNAs both including HST (18 sequences) and excluding HST (17 sequences) using MUSCLE. Then, the number of segregating sites (K) was estimated with DNA Sequence Polymorphism (DnaSP) v.5.10.01 (ref. 52). Finally, θ w was calculated as follows: K/a n /16,595, where a n is Σ i=1 n−11/i to take in consideration the number of mtDNA sequences in each data set. Adding HST to the 17 Neanderthal mtDNAs, the number of segregating sites almost doubled (from 78 to 145), whereas θ w increased from 1.37 × 10−3 to 2.50 × 10−3 (Supplementary Table 4). The latter value is closer to the mtDNA diversity estimated within three Denisovan mtDNAs (3.46 × 10−3)2.

BEAST analyses

We used the software package BEAST v.1.8.1 (ref. 33) to both estimate the divergence times between and within modern and archaic humans as well as to track the changes in the maternal effective population size (N e ) of Neanderthal mtDNAs through time.

For the skyline analyses, we first created a multiple genome alignment with only the mtDNA coding region of 18 Neanderthal mtDNAs and the rCRS as outgroup. We then removed from the alignment all columns where at least one mtDNA presented a gap or missing data, resulting in 15,345 positions. We run Modelgenerator v.85 (ref. 48) on our data set to identify Tamura-Nei 93 with a fixed fraction of invariable sites as the best-supported model. We set a fixed mutation rate (1.57 × 10−8 μ site−1 year−1)10 calculated for the coding region of modern humans with ancient mtDNAs as calibration points (Supplementary Note 4). As tree prior we selected the Bayesian skyline coalescent with 10 as group number and piecewise linear as the skyline model. We tested both a strict clock and an uncorrelated lognormal-distributed relaxed clock. For both models three MCMC runs with 50,000,000 iterations were run, with 10,000 sampling frequency. We discarded 10% of the states from each run as chain burn-in and then combined the three independent runs for both models using LogCombiner v.1.8.1 (included in the BEAST package), resulting in a total of 135 millions iterations. The two models were compared with a marginal likelihood estimation using path sampling (PS) and stepping-stone sampling (SS)34. The skyline tree prior in combination with a strict variation among tree branches performed better according to PS, while lognormal-distributed relaxed clock was best supported according to SS (Supplementary Table 7). The strict clock provided higher effective sample size (ESS) values because of earlier chain convergence; therefore, it was the chosen model to reconstruct a skyline plot for the 18 Neanderthal mtDNAs. We used Tracer v.1.6 selecting linear change as Bayesian skyline variant and a default of 100 as number of bins. In Supplementary Fig. 7 we report the mean N e (black line) and the 95% HPD interval (purple lines) of the Neanderthal mtDNAs in logarithmic value on the y axis and the temporal range from 350 to 32 ka on the x axis. We observe a N e reduction until around 42 ka followed by a rapid and short growth inversion, predating the Neanderthal disappearance.

For the dating analyses, we instead used a data set composed of 18 Neanderthal, 54 modern humans and 1 Denisovan mtDNAs as outgroup. As described above, we removed the D-loop from the alignment and further excluded all positions containing gaps and missing data for a total of 15,334 positions. The best-supported model for this data set was again Tamura–Nei 93 with invariable sites as tested in Modelgenerator v.85. The same fixed mutation rate was selected and tip dates were indicated for the eight dated Neanderthal mtDNAs (Supplementary Table 6). The date for all 54 mtDNA was kept as zero, while a range between 30 and 500 ka (initial value 50 ka) was given for all the undated Neanderthals. We tested two models of rate variation within branches: a strict clock and an uncorrelated lognormal-distributed relaxed clock. We also investigated two different tree priors: a Bayesian Skyline coalescent and a constant population size. As before, for each of the four model combinations MCMC was run three times with 50,000,000 iterations, sampling frequency 10,000 and 10% burn-in. The resulting 135,000,000 iterations per model were combined using LogCombiner v1.8.1. Best-supported model assessment was again implemented with marginal likelihood estimation using PS and SS. The strict molecular clock with the skyline tree prior provided higher likelihoods than the three other tested models (Supplementary Table 7).

Modelling the mitochondrial replacement

Under neutrality and assuming the Eurasian Neanderthal effective population size (N e ) to be relatively small (that is, <10,000 N e units), we calculated the mean time period necessary for an introgressing mtDNA lineage below 20% in frequency to reach fixation, when conditioning for that (Supplementary Table 8). This was computed using the following formula from Kimura and Ohta53, where N is the Eurasian Neanderthal N e and p is the proportion of the introgressing mtDNA lineage: T (p)=−2N (1−p) ln(1−p)/p that was readapted from 4N for autosomes to 2N for mtDNA in ref. 35. Generations were converted into years assuming a generation time of 29 years (Supplementary Table 8). We also calculated the likelihood of a complete replacement as a function of the Neanderthal effective population size and the admixture rate from a branch basal to modern humans (Supplementary Note 5).

Data availability

The HST mtDNA consensus sequence reported in this paper is available in GenBank with the accession code KY751400. All other data are available on request to the corresponding authors.