The earliest Denisovans and their cultural adaptation.

Since the initial identification of the Denisovans a decade ago, only a handful of their physical remains have been discovered. Here we analysed ~3,800 non-diagnostic bone fragments using collagen peptide mass fingerprinting to locate new hominin remains from Denisova Cave (Siberia, Russia). We identified five new hominin bones, four of which contained sufficient DNA for mitochondrial analysis. Three carry mitochondrial DNA of the Denisovan type and one was found to carry mtDNA of the Neanderthal type. The former come from the same archaeological layer near the base of the cave's sequence and are the oldest securely dated evidence of Denisovans at 200 ka (thousand years ago) (205-192 ka at 68.2% or 217-187 ka at 95% probability). The stratigraphic context in which they were located contains a wealth of archaeological material in the form of lithics and faunal remains, allowing us to determine the material culture associated with these early hominins and explore their behavioural and environmental adaptations. The combination of bone collagen fingerprinting and genetic analyses has so far more-than-doubled the number of hominin bones at Denisova Cave and has expanded our understanding of Denisovan and Neanderthal interactions, as well as their archaeological signatures.

Introduction

The identification and analysis of Pleistocene hominin remains form the basis for unraveling the processes governing human evolution, interaction, and adaptation, yet discovery of new human fossils continues to present a significant hurdle. Recent developments in excavation practices and archaeological science cannot subvert an unavoidable problem; that human remains are rarely identified, especially in prehistoric contexts where formal burials were not observed. This is particularly true for the Denisovans, a sister population to the Neanderthals, whose discovery fundamentally changed our understanding of hominin diversity in Eurasia during the Late Pleistocene. The high-coverage nuclear genome of a Denisovan individual (Denisova 3) showed that they diverged from a common ancestor with Neanderthals between 440-390 ka (thousand years ago). The identification of Denisovan ancestry in indigenous peoples of Australia and Papua New Guinea and in East and South East Asians has led to the inference that modern humans met and admixed with at least two distinct populations of Denisovans. This raises the possibility that Denisovans may have been widespread across continental Asia, island southeast Asia and near Oceania.

So far, only five small and highly fragmented fossils, all discovered at Denisova Cave (Russian Altai, Siberia, Russia), have been identified as Denisovans on the basis of DNA analyses. These include worn and incomplete molars (Denisova 2, Denisova 4, and Denisova 8), partial phalanges (Denisova 3), and small bone chips (Denisova 11). Only one (Denisova 3) has yielded enough DNA for whole genome sequencing . Poor DNA preservation and modern contamination has thus far impeded nuclear genome analyses of the other specimens. Outside Denisova Cave, a mandible from Baishiya Cave (Xiahe, China) was tentatively assigned to Denisovans on the basis of proteomic evidence and sediment DNA extracted from the site further confirmed the presence of Denisovans at the site.

Advances in proteomic research, in particular the increasingly common application of peptide mass fingerprinting (or ZooMS; Zooarchaeology by Mass Spectrometry), has been shown to be an efficient way for determining hominin presence at archaeological sites through the taxonomic identification of bone based on collagen characterization. In vertebrates, it is commonly used to assign genus or family-level identifications and, in some instances, species-specific determinations are possible. The highly time- and cost-efficient nature of ZooMS, its reproducibility, and the long-term preservation of collagen compared with other biomolecules, including DNA, make it an invaluable screening tool for the identification of fragmentary, morphologically non-diagnostic bones. ZooMS has been used to successfully identify hominin remains in large assemblages of bones including Denisova 11, a female individual with a Neanderthal mother and a Denisovan father.

Here we present a high-throughput application of peptide mass fingerprinting to unidentified bones from Denisova Cave. Located in the northwest Altai mountains, Denisova Cave preserves the longest archaeological sequence in northern Eurasia dating from the Middle Pleistocene to the Holocene. The cave contains a rich stratigraphic record, most notable for its Middle and Upper Palaeolithic cultural, faunal and fossil remains. It is the only site where the presence of Denisovans and Neanderthals has been determined on the basis of DNA recovered from both fossils and cave sediments in several layers throughout the sequence. In addition, the presence of early modern humans was recently confirmed at the site on the basis of mtDNA recovered from sediments .The combination of good biomolecular preservation, rich archaeological assemblages, and the presence of multiple hominin groups makes Denisova Cave one of the most informative archaeological sites for Pleistocene Eurasia.

Non-diagnostic bone fragments, an important untapped source of potential human fossils, represent 95% of the bones excavated at Denisova Cave. We applied ZooMS to 3,791 bone fragments from the East Chamber, one of the three explored galleries of the cave. The fragments were specifically chosen for their lack of diagnostic features which precluded macroscopic identification. The analysed bones came from each of the archaeological layers of the East Chamber, specifically layers 9, 11, 12, 13, 14, and 15. We also analysed a small number of bones from layer 17 which contains no archaeological evidence for hominin occupation (layers 10 and 16 are composed of culturally sterile deposits; Table S1). The majority of analysed bones were excavated from layers 14 and 15 from which no hominin bones were previously found, although layer 15, the lowermost archaeological layer of the East Chamber, has previously yielded Denisovan sediment DNA. From each bone, a chip of approximately 20 mg was removed and, following established ZooMS protocols, collagen was extracted and analyzed using a MALDI-TOF mass spectrometer in order to carry out taxonomic identification (Materials and Methods). The vast majority of the analyzed bones were assigned to large herbivores (Bos/Bison, Equidae and Cervidae) and carnivores, in reasonable agreement with fauna previously identified at the site through morphological analysis (Fig. S1).

Results

ZooMS

Five bone fragments (Fig. 1a) generated peptide mass fingerprints with characteristic markers corresponding to the Hominidae (Fig. 1b, Table S2, Dataset S1). Four of them come from layer 15 (DC7277, DC7795, DC8591, and DC8846) and one from layer 12 (DC4969). Given that no great apes are known from the region these bones almost certainly belong to humans. Human fossils identified using ZooMS now account for the majority of the hominin bones discovered at Denisova Cave (nine of the 17 fossils; 52%).

Figure 1

The five new human fossils from Denisova Cave identified using ZooMS analysis

a. Photographs of the new hominin bones. The two larger surfaces (front and back) of each bone are shown. b. MALDI-TOF mass spectra for the newly identified Hominidae bones, Denisova 17 (DC4969), Denisova 18 (DC7277), Denisova 19 (DC8846), Denisova 20 (DC7795), Denisova 21 (DC8591). Labelled peaks highlight the peptides used as markers to identify these bones as Hominidae.

microCT Analysis

To digitally preserve the morphology of the bone fragments, four of the five new specimens were scanned with a microCT system (BrukerTM SkyScan 2211 X-ray Nanotomograph). We used image spatial resolutions ranging between 0.020 and 0.023 mm, following the recommendations of Immel et al. to avoid degrading effects of X-ray irradiation on ancient DNA (Materials and Methods). 3D surfaces of the fossil bones were extracted from the microCT scans (Fig. S2; Dataset S2).

mtDNA Analysis

Since peptide mass fingerprinting cannot be used for a more specific taxonomic assignation than Hominidae, we used DNA analysis to identify the groups these five bones belonged to on the basis of mitochondrial (mt)DNA sequences. Extraction, sequencing and authentication of ancient hominin DNA from each bone followed published procedures (Materials and Methods). Using a mtDNA enrichment approach we isolated sufficient ancient hominin DNA and reconstructed the mitochondrial genomes of four of the five specimens; Denisova 17 (DC4969), Denisova 19 (DC8846), Denisova 20 (DC7795) and Denisova 21 (DC8591) (Table S3-S4). These were sequenced to an average coverage of 2,695-fold, 15-fold, 31-fold and 28-fold, respectively. Pairwise differences and phylogenetic analyses showed that the mtDNA of Denisova 17 falls within the diversity of Neanderthal mtDNAs, while the mtDNAs of Denisova 19, Denisova 20 and Denisova 21 fall within the diversity of Denisovan mtDNAs (Fig. 2, Table S5-S8). Denisova 18 contains too few ancient DNA fragments to securely associate its mtDNA with a hominin group (Materials and Methods and Supplementary Information).

Figure 2

Mitochondrial (mt) DNA maximum parsimony phylogenetic trees for the newly identified hominin bones.

a. Neanderthal mtDNA parsimony phylogeny including Denisova 17 and 26 previously published Neanderthal mtDNAs. The mtDNA of Denisova 17 differs by 27, 26 and 28 substitutions from the mtDNA of Denisova 11, Denisova 5 and Denisova 15. The topology of the tree corresponds to the 50% majority-rule consensus of 16 equally parsimonious trees rooted using the highly divergent mtDNA of the Hohlenstein Stadel Neanderthal (Supplementary Information, Fig. S5). b. Denisovan mtDNA parsimony phylogeny including Denisova 19, 20 and 21 and four previously published Denisovan mtDNAs. The tree is one of the two most parsimonious trees (Supplementary Information, Fig. S6) rooted using the mtDNA of the Middle Pleistocene Hominin from Sima de los Huesos. The trees are drawn to scale, with branch lengths computed in number of substitutions.

Discussion

The presence of Neanderthals in the Altai was originally identified in Okladnikov Cave, a site located 50 km to the north of Denisova Cave, on the basis of mtDNA evidence. Further archaeological and genetic data suggests that Neanderthals were in Siberia on several separate occasions. They appear at Denisova Cave (layer 12, East Chamber) at least ~150-130 ka. Five Neanderthal fossils have been found in the East Chamber so far, of which three are from layer 12 (Denisova 9, 11, 17) and two from the overlying layer 11.4 (Denisova 5, 15) (Fig. 3a). A single sediment sample from layer 14 of the East Chamber yielded Neanderthal DNA, but further work is required to replicate and confirm this signal. We estimated the molecular age of the mtDNA of the newly identified Neanderthal (Denisova 17) to ~134 ka (95% HPD: 94 – 177 ka) using Bayesian dating as implemented in BEAST v.1.10.4 and the mtDNA of 12 radiocarbon dated Neanderthal individuals as calibration points (Table S7). Phylogeny inferences show that the mtDNA of Denisova 17 is more distantly related to the mtDNAs of the two other Neanderthals from Denisova Cave, Denisova 5 and Denisova 15, who are more closely related to one another (Fig. 2a) (Supplementary Information, Fig. S4 and Table S5). In contrast, Denisova 11 mtDNA is more closely related to the mtDNAs of Neanderthals from western Eurasia and to other Siberian Neanderthals, such as those from Okladnikov Cave and Chagyrskaya Cave (Fig. 2a). Gene flow between Neanderthals and Denisovans provides additional indirect evidence of earlier interactions between the two groups. Analysis of the genome of a female Denisovan individual (Denisova 2), for example, has revealed that she had Neanderthal ancestry deriving from an introgression ~1500 years before she lived, as early as 250-200 ka. Two other Denisovans from higher up the stratigraphic sequence (Denisova 8 and 3) also show Neanderthal introgression from two different Neanderthal populations. Although it is not possible to tell where these interbreeding events occurred, they provide evidence for potential co-habitation and frequent interactions between the two hominin groups from >200 ka (Denisova 2) until their disappearance from the Altai around 50 ka (Denisova 3). Neanderthal presence, while more pronounced during the Last Interglacial at Denisova Cave (MIS5) (Fig. 3b), is discontinuous in the Altai region and may reflect occasional eastward migration of Neanderthal groups across large tracts of Eurasia. Since no gene flow from Denisovans to late European Neanderthals has been identified so far, these interactions seem most likely to have occurred in northeastern Eurasia. The Altai, in particular, appears to be an overlapping zone for both Denisovan and Neanderthal groups for over 150,000 years, witnessing and possibly facilitating population admixture as well as sustaining distinct hominin populations over this long period.

Figure 3

Stratigraphic and chronological relationship of the newly identified fossils from the East Chamber.

a. Stratigraphy of the East Chamber of Denisova Cave. The position of hominin fossils (circles) and sediment DNA (trowel) is shown. The newly identified hominin fossils, Denisova 17, Denisova 18, Denisova 19, Denisova 20, and Denisova 21 are shown next to the relevant stratigraphic layers they were excavated from. To the left, the colored bars and the numerals represent the age range in thousand years before present (ka BP) of each dated layer based on modelled optical ages. b. Updated age ranges of all hominins from Denisova Cave, including the newly identified specimens. The ages for Denisova 17, 19, 20, 21 derive from the Bayesian statistical treatment of optical ages. The ages for all other fossils derive from Bayesian modeling described in Douka et al. which includes both genetic, radiocarbon and optical ages. The marine-oxygen isotope curve compiled from benthic δ O records is shown to the top; the Last and Penultimate Interglacials (Marine Isotope Stages 5 and 7, respectively) that were the warmest parts of the last 300 ka, are highlighted.

The specimens with the Denisovan mtDNAs (19, 20 and 21) all come from layer 15 of the East Chamber. The mitochondrial sequences of Denisova 19 and 21 are identical indicating that they may belong to the same individual or be maternal relatives. They differ by 4 substitutions from Denisova 20’s mtDNA. In phylogenetic trees, the mtDNAs of the newly identified Denisovans form a clade with the mtDNAs of Denisova 2 (layer 22.1, Main Chamber) and Denisova 8 (layer 11.4, Easter Chamber) from which they differ by 20 and 30 substitutions, respectively (Fig. 2b and Fig. S3 and S6). Parsimony analyses are consistent with Denisova 19, 20 and 21 being of similar age or slightly older than Denisova 2, and substantially older than Denisova 8, and Denisova 3 and Denisova 4 (layer 11.2, East Chamber, and layer 11.1, South Chamber, respectively).

The mtDNA age estimates for the newly-identified fossils (Supplementary Information, Table S6) and their relationship to Denisova 2 agree with the overall stratigraphic context and previous attempts to cross-correlate the three Chambers of Denisova Cave on the basis of absolute dates, archaeological sequence and hominin groups. Previously, the earliest Denisovan (Denisova 2) was estimated to date to 122-194 ka using a Bayesian approach incorporating optical, genetic and stratigraphic data (Fig. 3b), or as early as 280 ka on the basis of optical ages only. That specimen was discovered in 1984 in the Main Chamber and its contextual integrity has been questioned, whereas the new fossils reported here were excavated in 2012-13 from a secure context. Layer 15 is the oldest archaeological layer of the East Chamber and is estimated to date to ~200 ka (205-192 ka at 68.2% probability, or 217-187 ka at 95.4% probability) based on Bayesian modelling of existing optical ages (Fig. 3a). Using these date estimates as calibration points in a Bayesian statistical framework, we inferred a divergence date for the mtDNAs of the three new and the four previously published Denisovans to ~229 (95% HPD 206 - 252 ka; Table S8) during the Interglacial period MIS 7. Both the mtDNA age estimates and the established chronology for layer 15 render Denisova 19, 20, 21, or their maternal relatives, the oldest Denisovans currently documented (Supplementary Information, Fig. S3).

The presence of individuals carrying Denisovan mtDNA in the lowermost archaeological layer 15 of the East Chamber offers us an opportunity to consider the wider archaeological and subsistence context of this group of hominins. So far this has not been possible because previous Denisovan fossils either derived from layers impoverished in archaeological material or from layers where Neanderthal cohabitation could not be excluded. Denisova 19, 20 and 21 date to the Penultimate Interglacial (MIS 7) (Fig. 3b); a warm climatic period with comparable conditions to today that would have rendered the Altai a favorable location for hominin expansion and intensified occupation. During this phase, a mosaic of landscapes can be detected in the vicinity of the cave, including both broad-leaved forests and open steppe landscapes. Both traditional zooarchaeological and ZooMS analyses revealed that the inhabitants of the cave targeted a variety of taxa living in these environments, including interglacial forest and forest-steppe species, such as roe deer (Capreolus pygargus), Siberian red deer (Cervus elaphus), and giant deer (Megaloceros giganteus), as well as species typical of more open country, such as horse (Equus ovodovi and Equus ferus), bison (Bison priscus), woolly rhinoceros (Coelodonta antiquitatis), and Mongolian gazelle (Gazella guttursza) (Fig. S1). Frequent anthropogenic impacts on bones, including splitting, burning and butchery cut-marks, confirm that these species were procured regularly. Humans appear not to have been the only occupants of Denisova Cave during this period, however. About a quarter of the macroscopically-identified faunal assemblage from layer 15 comprised carnivore remains, predominantly Canis lupus and Cuon alpinus. This high proportion of carnivore taxa suggests that humans may have been actively competing with these predators over resources and perhaps the cave itself.

Archaeologically, layer 15 (and layer 14) of the East Chamber contain the highest frequency of stone artifacts in the entire sequence of the cave, with more than 3,000 pieces per 1 m2. The lithic assemblage comprises discoidal, Levallois, and parallel cores to produce flakes using primary reduction techniques. Scrapers are the dominant tool type, including those shaped by steep Quina-type retouch, as well as spur-like, denticulate, and notched forms (Fig. 4). Large ventrally thinned and basally truncated flakes, or truncated-faceted flakes, are typical pieces (SI Section 4). A small number of blades with a longitudinal dorsal scar pattern are also present. Analyses of organic residues collected from a retouched flake from layer 15 revealed saturated and unsaturated fatty acids and, alongside the absence of bone and plant micro-residues, its proposed use was for animal skin processing activities, such as scraping, cutting, and/or sawing.

Figure 4

Denisovan lithic tools from the lowermost archaeological layers (15 and 14) of the East Chamber. 1, 2, 3, 7, side scrapers; 4 ventral thinned flake; 5, denticulate tool; 6, notched tool; 8, core.

Based on their techno-typological characteristics and chrono-stratigraphic position, the lithic assemblage of layers 14 and 15 of the East Chamber is attributed to an Early Middle Palaeolithic stone tool industry that has no direct counterparts in North and Central Asia. If we were to look further afield, the closest parallel is the Acheulo-Yabrudian cultural complex (AYCC) from the Near East. The AYCC has been identified at several cave (mostly) and open-air sites such as Tabun, Qesem, Hayonim, Misliya dating to between 400/350-250 ka. This is a period that marks the transition from the Lower to Middle Palaeolithic, and is linked to major transformations in hominid adaptive and cognitive abilities and major technological and subsistence innovations.39 These include, among others, the habitual use of fire and the systematic hunting and butchering of medium-size ungulates, such as fallow deer. Techno-typological similarities between the AYCC with Denisova Cave layers 14 and 15 of the East Chamber include comparable forms of ventrally thinned and basally truncated flakes, and the presence of Quina scrapers, denticulate and notched tools (examples in Fig. 4). There are no bifacial tools in the Denisova assemblage; bifaces are a typical element of the Acheulean variant of the AYCC, but are rare or absent in the other two facies of the complex. Yet, since there are no intermediate occurrences of similar traditions between the Levant and the Altai, and no hominin remains that could be directly linked to Denisovans outside the Altai and the Tibetan Plateau, further work is required to resolve issues surrounding Denisovan cultural adaptations and innovations. A focused attempt to characterize the lithic component of the earliest Denisovan layers is currently underway and will allow further understanding of the evolution of the Denisovan toolkit through time.

The distribution of Denisovan DNA in present-day humans suggests that Denisovans were widely dispersed occupying large tracts of Pleistocene Asia, and that there was spatial and temporal structure in their population. The Denisovan DNA introgressed in present-day humans from Siberia, East Asia, and indigenous Americans shares the highest similarity with the high quality genome of Denisova 3. However, the mtDNAs of the three older Denisovans we identify here, Denisova 19, 20, and 21, belong to a different mtDNA lineage to that of Denisova 3. Characterization of the nuclear DNA of these individuals is required to determine whether these early Denisovans are more closely related to the Denisovans that admixed with the ancestors of present-day humans living in Island Southeast Asia and New Guinea. Deciphering the relationship of various Denisovan groups is necessary to further understand how their distribution across Central, East and Southeast Asia may reflect the variability in the material culture observed in these regions during the Pleistocene. The challenges encountered by Denisovans while living in extremely diverse and changing environments, from the Altai mountains to the high-altitudes of the Tibetan Plateau, and possibly from north China to Island Southeast Asia, would have required adaptation in novel ways to survive. The application of state-of-the-art biomolecular approaches, such as palaeoproteomics and DNA analyses, to bone fossils and sediments holds great potential in identifying new hominins dating back to the Middle Pleistocene in Central Asia and elsewhere, and provides an opportunity to calibrate past demographic and dispersal events, while also linking them to the development of specific techno-complexes and cultural traditions.

Materials and Methods

Zooarchaeology by Mass Spectrometry (ZooMS)

Analysis was carried out at the ZooMS facility of the Department of Archaeology at the Max Planck Institute for the Science of Human History, Jena, Germany. We followed established protocols. In brief, from each bone approximately 10-20 mg was removed using a circular diamond drill-bit. Samples were rinsed in ammonium bicarbonate overnight and incubated for 1h at 65°C. The supernatant was treated with 0.4 μg trypsin (Thermo Scientific Pierce™ Trypsin Protease) and allowed to digest at 37 °C for 18 h. The incubated samples were concentrated and desalted using C18 ZipTips (Thermo Scientific Pierce™ C18 Tips) and eluted in a final solution of 50μl of 50% acetonitrile and 0.1% TFA. 0.5 μl of the resulting solution was mixed with 0.5 μl of α - cyano-4-hydroxycinnamic acid solution (10 mg/mL in 50% acetonitrile (ACN) and 0.1% trifluoroacetic acid (TFA) and allowed to crystallize. The samples were analysed using a MALDI TOF/TOF (Bruker Autoflex Speed LRF) mass spectrometer. The resulting spectra were screened for diagnostic markers using flexAnalysis 3.4 (Bruker Daltonics) and mMass software. The spectra were compared against a reference library of known peptide markers.

MicroCT Scanning

Prior to sampling for aDNA analysis, the five bones identified as Hominidae using ZooMS, were scanned with image spatial resolutions ranging between 0.020 and 0.023 mm using the BrukerTM SkyScan 2211 X-ray Nanotomograph housed at MPI-SHH in Jena, Germany. Following the recommendations of Immel et al. for avoiding degrading effects of X-ray irradiation, we strictly limited the scan image spatial resolution to 0.020 mm although smaller voxel sizes could have been achieved. We used a 0.5 mm titanium filter to remove the lowest energies of the X-ray spectrum. All scans were done applying a source voltage of 110 kV and a source current of 170 μA. Using the ‘Isosurface’ module the Avizo™ 9.4 (Visualization Science Group), we extracted 3D surfaces of the fossil bones from the microCT scans.

Mitochondrial DNA Analysis

DNA extraction and library preparation

After the removal of approximately 1 mm of surface material using a sterile dentistry drill, multiple small samples of ~7 to ~21 mg of bone powder were obtained from each specimen. DNA was extracted from each sample (or a sub-sample thereof, not using more than 15 mg bone powder) with a method that uses silica-coated magnetic particles for the retrieval of short DNA molecules on an automated liquid handling platform.42 Due to the low quantities of material that were removed, the volume of lysis buffer was reduced to 300μl, of which 150 μl were used for DNA extraction.

For Denisova 18, 19, 20 and 21, additional sampling of 9 mg to 17 mg of bone powder was performed and the samples were pre-treated with 0.5% sodium hypochlorite (bleach) solution following the protocol developed by Korlevic et al. DNA was extracted from the bleach treated samples following the same silica-based protocol as for non-bleached samples. The entire DNA extracts were converted into single-stranded DNA libraries. Extraction and library negative controls were carried through all steps of the experiments. The libraries were amplified according to a double indexing scheme and purified as described in the aforementioned library preparation protocol.

A total of 29 single stranded DNA libraries were made for the 5 samples, including 9 for which extracts were pre-treated with bleach (Table S3). Using quantitative qPCR, we estimated the number of DNA molecules incorporated in each library between 1.6x1010 and 5.6x109 for non-bleached samples and between 1.8x109 and 1.5x108 for bleached samples, which is on average higher than for libraries prepared from negative controls. Unfortunately, we found that the bleach treatment greatly reduced the amount of the endogenous DNA, making it unsuitable for these samples.

Mitochondrial DNA captures and sequencing

Each amplified library was enriched for human mitochondrial DNA (mtDNA) in two consecutive rounds of hybridization capture with a probe set covering the full human mitochondrial genome. The enriched libraries were pooled and sequenced on a MiSeq (Illumina) in 76-cycle paired-end run.

Data processing and Mapping to a reference genome

Base calling was performed using Bustard (Illumina) and sequences that did not exactly match the expected index combinations were discarded. Adapter sequences were removed and overlapping paired-end reads were merged using leeHom with the parameter « --ancientdna ». Overlap-merged sequences were mapped to the human mitochondrial revised Cambridge reference sequence (rCRS) using Burrows-Wheeler Aligner (BWA) with parameter « -n 0.01 −o 2 −l 16500 ». PCR duplicates were collapsed into single sequences by consensus calling using bam-rmdup (https://github.com/mpieva/biohazard-tools). Sequences shorter than 35 bases or with a mapping quality lower than 25 were discarded. Initial investigations based on the sharing of derived sites carried by sequences covering positions in the mtDNA genome that are diagnostic of modern humans, Neanderthals and Denisovans showed that the mtDNA of Denisova 17 is of Neanderthal type and that Denisova 18, 19, 20 and 21 are of Denisovan type. In order to recover sequences that may be difficult to map due to their divergence to the human reference mtDNA genome, we re-aligned the raw sequences of the DNA libraries of Denisova 18, 19, 20 and 21 to the mtDNA sequences of the Denisova 32 and the Denisova 88 individuals. Because sequences aligned to Denisova 8 mtDNA were slightly more numerous than the sequences mapped to the Denisova 3 mtDNA, we used them for subsequent analyses.

Reconstruction of mtDNA genome sequence

Data from different libraries of the same specimen were merged, and only sequences with a length greater than 35 base pair (bp) and a mapping quality of at least 25 were retained to call the mtDNA consensus sequences. Because of the level of present-day human DNA contamination in Denisova 19, 20 and 21, we restricted the analysis to sequences showing evidence of cytosine (C) to thymine (T) mismatches to the reference genome at the 3 first or last bases (deaminated sequences). Using deaminated sequences only, we called a consensus at each position covered by at least three sequences where at least 66% of the fragments carry the same base. The state of positions covered by two or fewer deaminated sequences were called using the alignment of all sequences when supported by more than 5 sequences for libraries with contamination estimate < 5%. For Denisova 17, we used all sequences and called a consensus at each position covered by at least 5 sequences where at least 80% of the fragments carry the same base. For all consensus sequences, manual correction of the alignment was necessary to confidently call certain positions, especially over cytosine homopolymer stretches.

Pairwise, phylogenetic analyses and relative molecular age estimates

For the phylogenetic analyses, the reconstructed mtDNA sequences of Denisova 17, 19, 20 and 21 were aligned to the mtDNA genomes of 26 Neanderthals, four Denisovans, the Middle Pleistocene Hominin from Sima de los Huesos, six ancient modern humans, , six present-day humans and a chimpanzee using Clustal Omega (Table S9). Pairwise differences between mitochondrial genomes and Neighbor Joining phylogeny were inferred using MEGA X. Maximum Parsimony analyses were conducted in PAUP* using minimization of F-value (MINF) as character-state optimization, gaps in the sequences were treated as “missing”. Optimal trees were generated by heuristic search method using a Nearest-neighbor interchange (NNI) branch swapping algorithm and consensus tree was called by 50% majority-rule method of several equally parsimonious trees. Relative molecular age of the Denisovan and Neanderthal mtDNAs were estimated by inferring by parsimony the number of substitutions accumulated in each mtDNA sequences since the split from the mtDNA of the middle pleistocene hominin from Sima de los Huesos for Denisovans and since the split from the highly diverged mtDNA of the Hohlenstein Stadel individual for Neanderthals. We caution that back mutations and multiple substitutions occurring at the same position will not be accounted for and may affect our inferences.

Bayesian phylogeny, molecular dating

Molecular dates of specimens and divergences were estimated using the dataset above with an addition of 60 modern human mtDNAs (Table S9) in a Bayesian statistical framework as implemented in BEAST v.1.10.4. The Tamura-Nei substitution model (TrN+I+G4) was determined as the best fit for the data using Model Test-NG. We used marginal likelihood estimations (path sampling method) to determine the best fitting clock model and tree model. For each model combination we used a chain length of 75,000,000 iterations with an alpha of 0.3 for the Beta distribution and a burn-in representing 10% of the chain. We used a mutation rate of 2.53x10-8 (95% CI: 1.76x10-8 - 3.33x10-8) substitutions per site per year for the whole mtDNA genome and 1.57x10-8 (95% CI: 1.17x10-8 - 1.98x10-8) substitutions per site per year for the coding region. Following the scale of Kass and Raftery, a strict clock and Bayesian skyline tree model were supported over the other model combinations (log10 BF > 5.2). The molecular age of Denisova 17 mtDNA was determined using as calibration points, the uniform priors for the ages of the Neandertals as defined in Peyrégne et al., 2019 (Supplementary Information). We inferred the age of the TMRCA of Denisovan mtDNAs using as calibration points the 95% confidence intervals of archaeological date estimates for the Denisovan remains from Douka et al., 2019 (Supplementary Information). We then performed six MCMC runs of 75,000,000 iterations for the complete mitochondrial genome and the coding region. We sampled trees every 2500 iterations after a burn-in of 10% of the number of iterations and merged the runs using BEAST post-analysis program logcombiner.