Pleistocene sediment DNA reveals hominin and faunal turnovers at Denisova Cave

Denisova Cave in southern Siberia is the type locality of the Denisovans, an archaic hominin group who were related to Neanderthals1–4. The dozen hominin remains recovered from the deposits also include Neanderthals5,6 and the child of a Neanderthal and a Denisovan7, which suggests that Denisova Cave was a contact zone between these archaic hominins. However, uncertainties persist about the order in which these groups appeared at the site, the timing and environmental context of hominin occupation, and the association of particular hominin groups with archaeological assemblages5,8–11. Here we report the analysis of DNA from 728 sediment samples that were collected in a grid-like manner from layers dating to the Pleistocene epoch. We retrieved ancient faunal and hominin mitochondrial (mt)DNA from 685 and 175 samples, respectively. The earliest evidence for hominin mtDNA is of Denisovans, and is associated with early Middle Palaeolithic stone tools that were deposited approximately 250,000 to 170,000 years ago; Neanderthal mtDNA first appears towards the end of this period. We detect a turnover in the mtDNA of Denisovans that coincides with changes in the composition of faunal mtDNA, and evidence that Denisovans and Neanderthals occupied the site repeatedly—possibly until, or after, the onset of the Initial Upper Palaeolithic at least 45,000 years ago, when modern human mtDNA is first recorded in the sediments.

The marginal log likelihoods for trees based on the mtDNA protein-coding genes generated from testing different clock and tree models with a path sampling approach…….……………………. 37 Supplementary Table 10| The estimated tip dates and divergence times (in years) reported from the Tracer program from BEAST2 for the protein-coding region………………………………………………………….…………………….. 37 Supplementary Figure 16| The phylogenetic mtDNA tree determined from Bayesian analysis with Beast2 using protein-coding genes from the mtDNA genome……………………………………………………………………………….  2 , together with provisional descriptions of the stratigraphy and sedimentology of the South Chamber profile sampled for optical dating 1 . Further excavations in South Chamber have revealed various stratigraphic complications, however, so analysis and interpretations of its Pleistocene stratigraphy, chronology, archaeology and environmental records are ongoing. We provide below a synopsis for each of the three chambers and highlight additional information of relevance to this study.

Main Chamber
The  The southeast profile sampled for sedimentary ancient DNA (aDNA) analysis did not reveal all layers (Supplementary Figure 1); layers 18-15, 13, 11.5 and 11.1 were not exposed at the time of sample collection. A total of 274 individual sediment samples were collected for aDNA analysis, including samples associated with the early Middle Palaeolithic, eMP (layer 22, n = 94; layer 21, n = 6; layer 20, n = 52), the middle Middle Palaeolithic, mMP (layer 19, n = 28; layer 14, n = 11; layer 12, n = 25), the Initial Upper Palaeolithic, IUP (layer 11, n = 47) and the Upper Palaeolithic, UP (layer 9, n = 11). The archaeological associations, environmental contexts and optical ages of these layers are summarised in Jacobs et al. (2019: Figure 3 and Extended Data Table 1) 1 . Figure 2b) include the yellowish, peak-shaped layer 22 near the base of the sequence, the reddish layer 12.3 higher up the sequence, and the distinct, almost horizontal contact between orangey layer 12.1 and greyish layer 11.4 that marks the boundary between the mMP and the IUP, which started 48-45 ka at Denisova Cave 2 .

Key stratigraphic features (Extended Data
Neanderthal DNA has previously been extracted from three of the sediment samples collected for optical dating from the southeast profile of the 1984 excavations (Jacobs et al., 2019: Extended Data Figure 1a,

East Chamber
The Pleistocene deposits in East Chamber comprise layers 17-9 and their subdivisions.   3 for four samples collected from the southeast profile of the 2013 excavations. Optical dating of 37 sediment samples 1 suggests that the deposits accumulated episodically from before 284 ± 32 ka to after 38 ± 9 ka (end dates for deposition of layers 17.1 and 11.1, respectively).
A total of 134 sediment samples for aDNA analysis were collected from the southeast profile exposed after excavation in 2015 (Extended Data Figure 2c; Supplementary Figure 2a). This includes samples associated with the eMP (layer 15, n = 4; layer 14, n = 35), the mMP (layer 13, n = 14; layer 12, n = 33; layer 11.4, n = 8; layer 11.3, n = 11) and the IUP (layer 11.2, n = 22; layer 11.1, n = 7). A further 118 samples were collected from the northwest profile exposed after excavation in 2016 (Supplementary Figure 2b), including samples from layers 15 (n = 8), 14 (n = 33), 13 (n = 36) and 12 (n = 41). Samples were not collected from layers 17, 16 or 9 in either profile. The archaeological associations, environmental contexts and optical ages of these layers are summarised in Jacobs et al. (2019: Figure 3 and Extended Data Table 1) Key stratigraphic features (Extended Data Figure 2c,d) include the yellowish, culturally sterile layer 17 at the base of the sequence, the blackish colour of layers 15 and 14 that distinguishes them from the overlying, yellowish layer 13, and the distinct charcoal band separating the wedge of pale cream-coloured layer 11.4 from the overlying, reddish layer 11.3. Single-grain equivalent dose (De) distributions suggest that deposits associated with layers 17-11.3 are largely intact, with minimal evidence for post-depositional mixing between layers 1 . Post-depositional disturbance (bioturbation) of some parts of layers 11.2, 11.1 and 9 is evident from the macro-stratigraphy and De distributions, but undisturbed areas are still present. Sediment samples collected for optical dating from the same profile as that sampled for sedimentary aDNA analysis have De distributions consistent with layers 11.4 (DCE16-3), 11.3 (DCE16-2) and 11.1 (DCE16-1) being intact at these locations (Jacobs et  Denisovan DNA was extracted previously from one of the sediment samples collected for optical dating (DCE12-12: 191 ± 26 ka) from layer 15 in the 2012 excavation profile 1,4 , providing the earliest evidence for hominin presence in East Chamber. Mitochondrial DNA was also extracted from two Denisovan teeth: Denisova 8 (from the interface between layers 12 and 11.4) and Denisova 3 (from layer 11.2). Neanderthal DNA was extracted from two of the sediment samples collected from layer 11.4 for optical dating: DCE14-13 and DCE14-15 (104 ± 12 and 123 ± 22 ka, respectively). The latter sample was erroneously indicated as having originated from layer 14 4 and was later correctly assigned to layer 11.4 1 ; it was displayed incorrectly as being from layer 14, however, in Jacobs et al.  1 . The bone fragment of a Neanderthal-Denisovan offspring (Denisova 11) was found in layer 12.3 and has a modelled age of 118-79 ka (95.4% HPD interval) 2 .

South Chamber
Excavations in South Chamber are currently in progress and further work to clarify the stratigraphy and chronology is ongoing. Information about environmental context, including the faunal assemblage, is not available at present and archaeological interpretations are preliminary. The provisional stratigraphic assignments given in Jacobs et al. (2019: Supplementary Table 18) 1 for the profile sampled for optical dating are under revision, and only two of the layers tentatively identified by them-layers 11 and 22, using their numbering system-are considered intact.
Ninety-one sediment samples were collected from the southeast (upper) profile (Supplementary Figure  3a,b). The Pleistocene sediments can be divided into rocky deposits preserved closest to the cave walls, and white-speckled, reddish sediments deposited in between. The rocky deposits have been assigned to layer 11 (Supplementary Figure 3b) and are associated with the IUP; 21 sediment samples were collected from this layer. The intermediate, reddish sediments have been significantly affected by post-depositional phosphatization, and we refer to them as 'phosphate deformation deposits' (pdd) (Supplementary Figure  3a,b). Extensive burrowing near the top of these deposits has also been observed in places. Tentative layer assignments were made in the field, based on associated archaeological finds, and should be considered provisional: pdd-9 (associated with the UP, n = 13 samples) immediately beneath the Holocene deposits, and pdd-12 (associated with the mMP, n = 57 samples).
A further 111 sediment samples were collected from the southeast (lower) profile (Supplementary Figure  3c,d). These sediments immediately underlie pdd-12 in the upper profile and are associated with mMP assemblages. A distinct stratigraphic boundary is evident between the yellowish deposits of layer 22, from which 25 samples were collected, and the overlying brownish sediments, from which 86 samples were collected. None of these deposits have been visibly impacted by post-depositional phosphatization, but excavations further into these deposits since 2017 have revealed complex post-depositional deformation that prohibits meaningful stratigraphic subdivisions; we refer to these deposits as 'deformed MP' (dMP) (Supplementary Figure 3d).
Owing to uncertainties in the stratigraphic attribution of the pdd (upper profile) and dMP deposits (lower profile), the corresponding samples cannot be interpreted in terms of stratigraphic or chronological patterns in the same way as can be done for Main and East Chambers. Accordingly, the aDNA data from South Chamber are not presented alongside those from Main and East Chambers, but they nonetheless carry valuable information about DNA preservation. None of the aDNA samples from Main and East Chambers were collected from parts of the deposit affected by post-depositional phosphatization, thus providing an opportunity to compare DNA preservation between chambers and between the upper and lower profiles in South Chamber. The aDNA data from South Chamber also provide preliminary insights into the presence or absence of hominin and other mammalian taxa, so a summary of these data is presented in Extended Data Figure 3.
To place the aDNA data for South Chamber in relative stratigraphic and chronological order, we have given tentative layer assignments to samples from the upper profile (pdd-9, layer 11 and pdd-12) and lower profile (dMP and layer 22). The latter samples are assumed to be stratigraphically lower-and therefore older-than those from the upper profile, but some dMP samples may overlap in time with those from pdd-12. The optical ages obtained for layers 22 and 11 indicate the approximate time span of sediment accumulation, with deposition of these layers estimated to have ended 269 ± 97 ka and started after 47 ± 8 ka, respectively 1 .
No Neanderthal remains have been found in South Chamber, but a Denisovan tooth (Denisova 4) was recovered from layer 11 in 2000. Mitochondrial and nuclear DNA have been extracted from this tooth, which has a modelled age of 84-55 ka (95.4% HPD interval) 2 .

SUPPLEMENTARY SECTION 3: ANCIENT DNA PRESERVATION IN SEDIMENTS
The ability to perform high-density sampling of sediment at a single site also enables the examination of potential trends in the preservation and degradation of ancient DNA in relation to both time and different sediment characteristics, while controlling for temperature fluctuation and geographical location. Previous studies relating DNA preservation and sediment are mainly focused on binding properties of DNA to different minerals from laboratory prepared samples 4,6,7 . Studies containing samples collected in the field have been focused on the presence or absence of ancient DNA and do not include in-depth evaluations of DNA preservation and degradation under different conditions or correlations to sediment characteristics [8][9][10] . Evaluations of DNA degradation from ancient skeletal remains have shown that deamination and fragment size, which are the key parameters frequently used to describe ancient DNA degradation, do not correlate very strongly with sample age in studies involving material from multiple archaeological sites 11,12 13 . Stronger correlations between age and DNA degradation have been observed when analysing material of different ages from a single site 14 . The data set produced in this study thus provides ideal conditions for studying the effect of age, as well as various physical and chemical properties of the samples, on DNA preservation in sediment.
In this section, we use the following four parameters to characterize DNA preservation in sediments from Denisova Cave: The frequency at which cytosines in the reference genome are substituted by thymines (C-to-T substitutions) at the terminal positions of sequence alignments provides a proxy for the level of deamination in a library 15 . The reference in each family that resulted in the most identified fragments was used for each library (Supplementary Data File 1). Initial analyses showed that deamination rates can vary across sequences from different families, even in the same sample, which may be the result of mapping biases (e.g., differences in evolutionary distance between the sequenced fragments and the reference genome), real differences in the properties of DNA molecules, or both. We therefore provide deamination rates separately for the mammalian families that are most abundantly represented throughout the stratigraphy (Bovidae, Canidae, Hyaenidae and Ursidae). We also restricted our analyses to the 5' ends of sequences, which are less prone to ligation biases in single-stranded library preparation than 3' ends 16 .
(ii) Average fragment length Average fragment length was computed using all unique mammalian mtDNA sequences from families identified as ancient in each library. Sequences shorter than 35 bp are not used in this calculation, as they were not subjected to taxonomic identification.
(iii) Total number of mammalian mtDNA fragments This measure reflects the number of unique mammalian mtDNA sequences that were assigned to ancient taxa, normalized by the amount of material (milligrams of sediment) used for DNA extraction. It is important to note that this value is not independent of sequence depth. For some libraries, not enough sequences were generated to ensure the sequencing of all or most unique mtDNA fragments (see Supplementary Data File 1 for duplication rates).
(iv) Inhibition The co-extraction of inhibitory substances may reduce the efficiency of library preparation, which was determined using the conversion rate of a spiked-in control oligonucleotide 17 .

Dependence of DNA preservation on stratigraphic depth (age)
We tested whether there is a correlation between the relative age of the sediments, as reflected by the layer from which samples were collected, and each of three distinct measures of DNA degradation, using Spearman's correlation test (data determined not to be normally distributed using a Shapiro-Wilk normality test). First, in order to determine whether deamination rates increase with the relative age of the sediments, we computed the frequency of C-to-T substitutions by layer in the depositional sequence, from top to bottom for each chamber. In all three chambers, and for each ancient mammalian family tested, a significant positive correlation (maximum p-value = 2.7e-08, minimum n = 133) was found between increasing depth in the stratigraphy as represented by layer and deamination rates (Extended Data Figures  3d and 4a Table 1). In summary, significant correlations were found between all parameters of DNA degradation studied here and depth of the sampling locations in the stratigraphy for all three chambers.

Dependence of DNA preservation on sediment pH
It is known that depurination and other mechanisms contributing to DNA decay proceed faster at low pH 18 .
To determine whether the pH of sediment affects the preservation of ancient mammalian DNA at Denisova Cave, the pH was measured for a subset of sediment samples from each layer and chamber. For this purpose, at least 1 g of sediment was aliquoted into a 2 mL tube. The sample was then weighed on a weigh boat and left to air dry at room temperature for about 1.5 hours. The weight after drying was recorded. The dried sediment was then homogenized using a mortar and pestle. Each sediment sample was split in half (based on its weight), transferred into a fresh 2 mL tube and resuspended through vortexing (2-3 minutes) in a volume of water corresponding to approximately twice the sample weight (e.g., 1 ml water for 500 mg sediment). The samples were left to settle for at least 60 minutes and were then centrifuged in a table-top centrifuge (2 minutes at maximum speed). The supernatant was transferred to 5 mL tubes and pH measurements were performed using a Mettler Toledo Seven Easy pH meter. Three measurements were taken for each original sample: one for each half of the sample, and the third after pooling the supernatants from both halves. For some sediment subsamples, the volumes of supernatant remaining after the settling phase were below the specifications of the pH meter, so only the pH values resulting from the pooled supernatants were used for analysis. The associated data can be found in Supplementary Data File 2.
For the vast majority of samples, we obtained pH values in a narrow range between 7.5 and 8 (Supplementary Figure 4). Only four out of 58 samples produced neutral or slightly acidic pH values (between 6 and 7). Of those, only one yielded small amounts of ancient mammalian mtDNA. A significant correlation was found between the amount of DNA recovered and the pH of the sample (Spearman's correlation test: rho value = 0.39, p-value = 0.003, n = 58). When controlling for time, as represented by layer (excluding samples that cover multiple layers), this correlation remains significant (multiple regression analysis: p-value = 0.013, n = 46). We also tested whether the inferred level of inhibition in library preparation correlates with the pH of the samples, to investigate whether low or high pH values may lead to the co-extraction of inhibitory substances, but found no correlation (Supplementary Figure 5) (Spearman's correlation test: rho value = 0.017, p-value = 0.90, n = 57). Taken together, these results indicate that low pH negatively impacts DNA preservation in sediments, in agreement with a previous study on lake sediments 19 .

Dependence of DNA preservation on post-depositional phosphatization
The presence of phosphates can impact the pH of sediments and may interfere with DNA binding to minerals. We tested this by comparing the presence of post-depositional phosphatization in the Pleistocene layers of each chamber (Supplementary Table 2) 1,3,20 to the DNA characteristics described above: the deamination frequency, average DNA fragment size, number of fragments assigned to ancient taxa per mg of sediment and the inferred inhibition level ( Supplementary Figures 6 and 7). Many of the sedimentary layers at Denisova Cave contain localized zones of phosphatization (e.g., phosphatic rinds around limestone clasts 3 ), but areas of extensive phosphatization in the Pleistocene deposits are restricted to parts of the uppermost layers in Main Chamber (layer 9) 1 , East Chamber (layers 9 and 11.1) 20 and South Chamber (pdd-9 and -12, where 'pdd' denotes phosphate deformation deposits; Extended Data Figure 3a) 1 .
There was no significant difference in the amount of ancient DNA recovered or inferred inhibition observed when comparing layers with and without extensive phosphatization when using a Wilcoxon test (data not normally distributed based on a Shapiro-Wilk test) and correcting for multiple testing 21 (amount of ancient DNA: p-value = 0.15, n = 664; inhibition: p-value = 0.7, n = 707, respectively). However, the DNA recovered from these layers had longer average fragment sizes and lower deamination rates (average size: p-value = 3.7E-07, n = 664; 5' deamination: p-value = <2E-16, n = 554 (Canidae), p-value = <2E-16, n = 587 (Bovidae), p-value = <2E-16, n = 473 (Hyaenidae), p-value = 3.5E-15, n = 493 (Ursidae); Wilcoxon test corrected for multiple testing 21 ). These trends remain significant when controlling for time, as represented by layer, with a type III anova test (average size: p-value = 3.36E-04, n = 572; 5' deamination: p-value = 8.07E-06, n = 399 (Canidae), p-value = 2.89E-06, n = 410 (Bovidae), p-value = 3.12E-05, n = 362 (Hyaenidae), p-value = 5.1E-12, n = 385 (Ursidae)).The longer fragment sizes may be due to competition between phosphate in the DNA backbone and free phosphate in the sediment for binding to mineral surfaces, or because these are younger layers. The presence of phosphate is linked to the pH values of sediments, as the heavily phosphatized layers produced the lowest measured pH values from all layers across all three chambers (Supplementary Figure 4). Phosphates can acidify neutral or basic soils, so post-depositional phosphatization may be the driving factor behind the observed decrease in pH. However, the Pleistocene layers with the most extensive phosphatization are also located closest to the interface with the overlying Holocene deposits, so the pH of the uppermost Pleistocene sediments may also have been affected by diagenetic changes associated with the microbial degradation of bat guano in the lowermost Holocene deposits 20 .

DNA preservation and clast size
The size range of clastic material comprising each layer could relate to the surface area on sediment particles available for the binding and preservation of DNA. Most of the Pleistocene layers in Denisova Cave are composed of poorly sorted silt and sand grains washed or blown into the cave, or reworked from preexisting cave sediments, interstratified with variable amounts of angular clasts of spalled limestone (gravel, cobble and boulder in size) 1 . We examined the average DNA fragment size, number of fragments assigned to ancient taxa per mg of sediment, inferred inhibition and deamination for samples collected from various layers in Main and East Chambers (Supplementary Table 2 and Supplementary Figures 8 and 9). DNA was recovered from layers composed predominantly of fine-grained sediments (e.g., layers 9 and 21 in Main Chamber) and from layers containing larger proportions of coarser clasts (e.g., parts of layer 11 in both Main and East Chambers), which indicates that DNA can be preserved in deposits containing clasts of a variety of sizes.

DNA preservation and sediment colour
The colour of sediment may indicate the presence of different types of organics (e.g., charcoal, bone and coprolite fragments), metals and other substances that could affect DNA preservation. We compared the colour of each layer 1 to the deamination rate, average fragment size, ancient mammalian mtDNA fragments per milligram sediment, and inferred inhibition for samples collected from Main and East Chambers ( Supplementary Figures 10 and 11). While no significant trends were found in relation to inferred deamination rates, blackened sediments had significantly shorter and less ancient DNA fragments than other sediment colours, except for rusty ochre or yellow (Supplementary Tables 5 and 6). Blackened sediments showed decreased library preparation efficiency in comparison to all other sediment colours, except for grey (Supplementary Table 6). It was possible to recover DNA from layers of all colours, but those black or yellow in colour produced the lowest yields. However, each sediment colour was observed in no more than three layers (Supplementary Table 1) and typically deposited over a limited time span. Blackened, yellow and rusty ochre sediments all originate from the lowest layers in Main and East Chambers, so these results are likely driven by the impacts of time. Further data from additional sites are needed, therefore, to ascertain whether the differences observed here are attributable solely to sediment colour or whether they result from uneven sampling of the stratigraphy.  Each point represents one sample from Main, East or South Chamber for which pH was determined. Library preparation efficiency was determined by comparing the number of spike-in control molecules recovered in the relevant library to that in the library negative controls. Samples with an efficiency lower than 0.5 are considered to show inhibition. The points are coloured depending on if evidence for ancient DNA preservation was found in the sample. Note that technical variation in qPCR measurements may lead to library preparation efficiencies greater than 1.  The different clastic sizes were defined as coarse and fine (n=76), coarse with boulders (n=7), coarse (n=64), and fine (n=109). Each point represents a single sample and is coloured depending on if that sample was identified as containing ancient DNA. The box plots show the distribution of the observed average fragment lengths, number of fragments assigned to ancient taxa per mg of sediment, and library preparation efficiency following the standard Tukey representation (box limits represent the lower and upper quartiles, whiskers represent 1.5 times the interquartile range, outliers are represented by black dots). The "NA" assignment indicates that there was no published information for this layer (n=408). Note that technical variation in qPCR measurements may lead to library preparation efficiencies greater than 1. Supplementary Figure 10 | Inferred 5' deamination rates of ancient sequences assigned to Bovidae (n=435), Canidae (n=442), Hyaenidae (n=334) and Ursidae (n=393) across the layers in East and Main Chambers stratified by the recorded colour of sediments (black/brown, n=16; brown, n=131; brown/red, n=72; grey, n=39; greyish brown, n=132; rusty ochre, n=64; yellow, n=29). Each point (centre of error bar) represents the average observed 5' C to T substitution frequency in a library from a specific layer in the relevant chamber. The grey bars represent the 95% binomial confidence intervals based on the number of alignments starting at a C in the reference genome. The box plots show the distribution of the observed 5' C to T substitution frequencies following the standard Tukey representation (box limits represent the lower and upper quartiles, whiskers represent 1.5 times the interquartile range, outliers are represented by black dots). The "NA" assignment indicates that there was no published information for this layer.

Supplementary
Supplementary Figure 11 | Average fragments size, the number of ancient mtDNA fragments recovered per milligram sediment and library preparation efficiency in samples from East and Main Chambers stratified by the recorded colour of sediments (black/brown, n=16; brown, n=131; brown/red, n=72; grey, n=39; greyish brown, n=132; rusty ochre, n=64; yellow, n=29). Each point represents a single sample and is coloured depending on if that sample was identified as containing ancient DNA. The box plots show the distribution of the observed average fragment lengths, number of fragments assigned to ancient taxa per mg of sediment, and library preparation efficiency following the standard Tukey representation (box limits represent the lower and upper quartiles, whiskers represent 1.5 times the interquartile range, outliers are represented by black dots). The "NA" assignment indicates that there was no published information for this layer (n=164). Note that technical variation in qPCR measurements may lead to library preparation efficiencies greater than 1.

Identification of diagnostic positions
In order to identify different groups of hominins throughout the stratigraphy, we determined diagnostic positions in the hominin mtDNA tree, following a strategy described earlier 22 . We first created a multiple sequence alignment, using MAFFT 23 , of 19 Neanderthal mtDNA genomes [24][25][26][27][28][29][30][31] , including that of the Hohlenstein-Stadel (HST) Neanderthal (which was retained as a separate branch due to its divergence from other, 'typical', Neanderthal mtDNA genomes) 32 , a Middle Pleistocene hominin from Sima de los Huesos 33 , 4 Denisovans 5,34-36 , a chimpanzee 37 , a set of 311 present-day modern humans from a wide geographical distribution 26 and the revised Cambridge Reference Sequence (rCRS) 38 (see Supplementary Data File 3 for a list of genomes used). We then determined positions where all mtDNA genomes representing one branch in the tree show a base difference to all other mtDNA genomes. In order to maximize the number of diagnostic positions available without reducing the accuracy of lineage assignment, two parameters were explored when determining diagnostic positions. First, we included or excluded the chimpanzee, and second, we required either at least 99% or 100% of the 311 humans to share the same state in order to call a diagnostic position. The combination of these parameters led to four sets of diagnostic positions which were tested further (Supplementary Figure 12).
Each set of diagnostic positions was tested in two different ways. First, we determined the support for each branch using previously published sequence data from a sediment sample from Denisova Cave, which had been shown to contain Neanderthal DNA 4 , as well as unpublished data from an ancient modern human bone. Sequences from hominin mtDNA fragments were isolated using the analysis pipeline described in Methods and the support for each branch was examined using fragments overlapping the diagnostic positions. Only mtDNA fragments carrying C-to-T substitutions at the 5' and/or 3' end (putatively deaminated fragments) were included in the analysis in order to deplete DNA fragments originating from present-day human contamination. With all four sets of diagnostic positions, we observed high support (greater than 96%) for the typical Neanderthal and modern human branches, respectively, and consistently less than 10% for all other branches, for most even less than 2% ( Supplementary Figures 13 and 14). Based on the above test, we concluded that all four sets of diagnostic positions produced valid and nearly indistinguishable results. We therefore chose the set with the largest number of diagnostic positions for further analyses, i.e. the set that was generated without the chimpanzee and for which 99% of the modern human mtDNA genomes were required to agree on one base.
Second, the homogeneity of coverage along the mitochondrial genome was evaluated by merging all putatively deaminated sequences generated in this study from sediment samples that were identified as containing ancient hominin DNA and determining the coverage for each position in the mtDNA genome. To evaluate the possibility of increased coverage at more conserved areas of the genome due to capture bias, we compared the observed coverage at each position to the phyloP conservation score [39][40][41] . A significant correlation between coverage and the phyloP score (Spearman's correlation test: rho = 0.0244, p-value = 1.7e-3) was found, but when we plotted a histogram of coverage at all diagnostic positions (determined without chimpanzee and requiring at least 99% of humans to be identical at the position) we observed only four positions with a coverage greater than 2 standard deviations from the mean. These diagnostic positions were removed to reduce imbalance in the contribution of individual positions in subsequent analyses (Supplementary Figure 15).

Assignment of mtDNA fragments to known hominin mitochondrial groups
Libraries that were identified as containing ancient hominin DNA were then evaluated for support of known hominin mitochondrial groups (human, Denisovan, Sima de los Huesos, HST-like Neanderthal, other Neanderthals). This was done by determining the proportion of all hominin mtDNA fragments sharing the group-specific state at positions diagnostic for the respective group. The support for each group was computed twice: (i) using all fragments overlapping diagnostic positions and (ii) using fragments carrying a C-to-T substitution within the first or last three positions in their alignment to the rCRS (putatively deaminated fragments). In order to assign mtDNA fragments from a library to one or more groups, the following criteria had to be met: (i) Fragments had to support the group-specific state at three or more diagnostic positions, and (ii) significantly more than 10% of the fragments had to support the group-specific state based on 95% binomial confidence intervals. For the identification of modern humans, only deaminated fragments were used in order to minimize the impact of present-day contamination with modern human DNA.
In the initial screening, which included only one library per sediment sample, 124 of the 168 samples (74%) that showed evidence for the preservation of ancient hominin mtDNA could be assigned to specific hominin groups (Supplementary Data File 1). In order to understand why the lineage assignments were not resolved for all samples, the 44 libraries identified as containing ancient hominin DNA, but not assigned to a specific lineage were examined further. Thirty-two of these libraries contained high amounts of modern human contamination (≥ 50% based on the support for the modern human branch when using all fragments, i.e., without restricting to putatively deaminated fragments only), low numbers of hominin mtDNA fragments (≤40 deaminated fragments), and/or low-level support for multiple lineages (>5%, but non-significantly more than 10%), all of which likely hindered the ability to make a lineage assignment.
To test the reproducibility of identifying ancient hominin DNA within the same sediment sub-samples, we repeated the DNA extraction, library preparation, hybridization capture and sequencing on aliquots of lysates that originally produced libraries that tested negative (n = 19) and positive (n = 61) for ancient hominin DNA. Of the 61 samples that initially tested positive, 28 (46%) were positive in the second experiment (Extended Data Figure 6). However, our ability to replicate results varied strongly with the number of putatively deaminated DNA fragments that were recovered in the first screening: 16 out of 17 (94%) of the samples that initially yielded more than 100 fragments tested positive again (the sample that failed replication produced 105 fragments in the first screening), whereas positive results could be replicated for only 2 out of 31 (6%) of samples that showed less than 50 fragments in the first screening. Lineage assignments were consistent across experiments for libraries from the same lysates where these could be made, with the exception of a single lysate from the deformed Middle Palaeolithic portion of South Chamber (S149) that showed Denisovan support in the first and both Denisovan and Neanderthal support in the second library. Of the 19 lysates that initially tested negative, 5 (26%) turned positive in the second screening. These results indicate that the hominin DNA content of many samples is close to the detection limit, where our ability to detect ancient hominin DNA is hampered by stochasticity in the sampling of molecules, small fluctuations in the efficiency of sample preparation and variations in the amount of present-day human contamination that is introduced during laboratory work. It should also be noted that not all libraries from the second lysate aliquot were sequenced deeply enough (duplication rate > 3) to recover most of the unique fragments present in the library, which likely further impaired our ability to replicate positive results from the first screening.
After the sequencing of additional libraries produced to test for reproducibility, sequenced DNA fragments from each original lysate were merged for the remaining analyses (see Supplementary Data File 1 for more information on which samples had data merged). Based on these data, 142 of the 174 samples (82%) identified as containing ancient hominin mtDNA were assigned to a specific mitochondrial group. Of these, 45 were assigned to non-HST-like Neanderthals, 66 to Denisovans, 2 to both, 35 to modern humans and 2 to both modern humans and non-HST-like Neanderthals. No samples were assigned to HST-like Neanderthals or the Sima de los Huesos lineage.
The homogeneity of DNA preservation within samples was also tested by taking between two and seven additional sub-samples from 42 of the sediment samples that had been shown to contain ancient hominin DNA. All of the 7 samples that had more than 100 deaminated fragments in the initial library yielded hominin DNA in at least 50% of the sub-samples. In contrast, this was the case for only 10 out of 20 samples (50%) that yielded less than 50 deaminated fragments in the first library that was prepared. The lineage assignments were not always consistent for sub-samples of the same sediment samples, especially for those with low amounts of deaminated fragments in the original library (Extended Data Figure 6c). These observations indicate that even though some samples tend to produce richer libraries upon repeated subsampling than others, preservation of ancient hominin DNA varies substantially within samples.
When combining all data, 244 sediment sub-samples (out of a total of 868) were found to contain ancient hominin DNA, of which 194 (80%) could be assigned to a lineage. From these 194 sediment sub-samples, 79 were assigned to non-HST-like Neanderthals, 87 to Denisovans, 6 to both Neanderthals and Denisovans, 2 to both modern humans and Neanderthals, 36 to modern humans, and none to HST-like Neanderthals or the Sima de los Huesos lineage (see Figure 1 and Extended Data Figure 3 for an overview of lineage support across the layers of all three chambers).

SUPPLEMENTARY SECTION 5: DETERMINING THE NUMBER OF SEQUENCE VARIANTS
Previous work on hominin mtDNA from sediments from Denisova Cave has shown that in most cases multiple individuals contributed to the DNA that was isolated from sediment. However, in one instance a single sequence variant was recovered, likely originating from only one individual 4 . We thus aimed to determine whether additional such cases can be detected in the much larger data set generated in the present study. We repeated the processing of hominin mtDNA sequences from sub-samples identified as containing either Neanderthal and Denisovan DNA using the Vindija 33.19 Neanderthal mtDNA genome (NC 011137.1) and the Denisova 3 mtDNA genome (NC_013993.1), respectively, as reference for mapping instead of the rCRS in order to minimize the loss of fragments due to sequence divergence to the reference. We estimated present-day human contamination based on the support for the modern human lineage in all fragments (Supplementary Section 4). For sub-samples with less than 5% estimated present-day human contamination all mtDNA fragments were used, for sub-samples with more than 5% estimated present-day human contamination putatively deaminated fragments were used. We then computed the mtDNA coverage from each sub-sample. Only sub-samples with at least 3x coverage (n = 7) were used for estimating the number of sequence variants (Supplementary Table 7).
We then used a previously published maximum likelihood method 4 that uses the consistency of the observed bases at each position to estimate the likelihood for a given library to contain one, two or three sequence variants. Hominin DNA fragments from three samples were estimated to originate from single contributors. In the remaining four samples, two sequence variants are most likely present (Supplementary Table 7). In all cases, it was estimated that the major sequence variant makes up at least 90% of the sequences generated from each sample. We therefore used all seven libraries for consensus calling (Supplementary Section 6) and reconstructing phylogenetic relationships with previously published mtDNA genomes (Supplementary Section 7). Table 7 | Estimated coverage and number of sequence variants present based on a maximum likelihood analysis for libraries with at least 3-fold coverage. The yellow highlights indicate when deaminated or all fragments were used for consensus calling. The support for the human lineage was used as a proxy for evaluating contamination. Data source indicates if data from a single library or merged data from the same lysate (sub-sample) were used.

SUPPLEMENTARY SECTION 6: RECONSTRUCTING ARCHAIC HUMAN MTDNA CONSENSUS SEQUENCES
To perform phylogenetic analyses, we attempted to reconstruct consensus genome sequences from the seven sub-samples identified as containing DNA predominantly from a single individual (Supplementary Section 5). For sub-samples containing Neanderthal or Denisovan DNA, alignments to the Denisova 3 or Vindija 33.19 mtDNA genomes, respectively, were used for consensus calling, while for samples containing ancient human mtDNA we used the alignment to the rCRS. Consensus calls were made at positions covered by at least 3 fragments if 60% of the fragments supported the majority base. In order to minimize the impact of deamination, T's on the terminal three bases of each fragment were disregarded. This strategy yielded mtDNA consensus sequences with between 217 and 6745 missing positions due to low coverage and up to 11 missing positions per sample due to low support of the majority base (Supplementary Table  8).
As the mtDNA capture probes were designed based on the human mtDNA genome, they may cause a bias towards modern human contamination in regions of the mtDNA genome where the divergence between modern and archaic humans is highest (the D-loop). Undetected contamination with mammalian mtDNA, on the other hand, is likely to be highest in the most conserved regions of the mtDNA genome (the 12S and 16S rRNA genes). We therefore decided to only use the 13 protein-coding genes of the mtDNA genome for further analysis to minimize the potential impact of erroneous consensus calls on phylogenetic reconstructions.
After limiting our analysis to the 13 protein-coding genes, we examined the missing positions that were due to low support of the majority base. We found that some of them were positions where some fragments and the reference genome carried a C and other fragments a T. This indicates that these positions likely remained unresolved due to deamination that occurred outside of the terminal three bases of DNA fragments. We manually corrected positions in the consensus sequences that showed low consensus support due to deamination in the 4 th to 6 th position at fragments ends (5 instances), resulting in consensus sequences with between 13 and 3889 missing calls due to low coverage and between up to 8 missing calls due to low consensus support for the approximately 11.3 kb of protein-coding genes in the mtDNA genome (Supplementary Table 8

Constructing a Neighbour-Joining tree
To place the four consensus mtDNA genome sequences we reconstructed within the variation of hominin mtDNA genomes, we generated a multiple sequence alignment containing the mtDNA genome sequences of 55 present-day (from 26 and the rCRS) and 10 ancient modern humans 43-51 , 24 Neanderthals 2,24-32 , 4 Denisovans 5,34-36 , and a chimpanzee 52 using MAFFT 23 . We then removed all non-protein-coding regions and positions containing gaps or missing data in one or more sequences. Trees were inferred using the Neighbour-Joining method in MEGAX 53 , with 500 bootstrap replicates to estimate the support for each node. This analysis was carried out separately for the three Neanderthal and the one Denisovan mtDNA genomes reconstructed from sediment sub-samples (Extended Data Figure 6a,b).
All three newly reconstructed Neanderthal mtDNA genomes cluster together. This cluster falls outside of the Altai-like Neanderthal clade, and outside of the branches leading to other Neanderthal mtDNA genomes recovered from ancient individuals living in the Altai (Chagryskaya 8, Denisova 11 and Okladnikov 2), further emphasizing the Neanderthal mitochondrial diversity in the region during the late Pleistocene. The newly reconstructed Denisovan mtDNA genome falls outside of both the Denisova 2 and 8 and Denisova 3 and 4 clades, implying that there is more diversity in Denisovan mtDNA genomes than previously known within Denisova Cave.

Branch shortening estimates
One of the sub-samples (from sample M65, layer 19 of Main Chamber) produced a Neanderthal consensus mtDNA genome sequence of a sufficient quality to date it by molecular methods. We therefore constructed a phylogenetic tree using BEAST2 54 and the multiple sequence alignment described in the previous section but excluding other consensus sequences from sediment samples. We identified the best fitting clock and tree model for the analysis by using a path sampling approach from the MODEL_SELECTION package 55-57 in BEAST2 54 . For each model combination, 40 path steps were used, each with a chain length of 25,000,000 iterations, parameter alpha of 0.3, pre-burn-in of 75,000 iterations and an 80% burn-in of the whole chain. A mutation rate 1.57 X 10E-8 was used as the analysis was restricted to the protein-coding region. For both the relaxed log normal and strict clock models, a normal distribution was used with the mean set to the mutation rate mentioned above and a sigma of 1.00E-10 45 . For all models, the substitution model Tamura-Nei 1993 (TN93) 58 was used as this was estimated to be the best model for archaic hominins 32 . Previously published radiocarbon dates for ancient modern humans and Neanderthals 2,24-32,59 were used to calibrate the tree. All modern samples were set to present day (date = 0). The ages of Neanderthals of unknown age, including the sediment sample, were constrained to a range of 30,000 to 200,000 years, with the exclusion of Sima de los Huesos which was constrained to 200,000 to 780,000 years. As has been done elsewhere 32 , Denisova 3 was restricted to the range of 30,000 to 100,000 years based on previous molecular dating of its nuclear genome 60,61 . For the other Denisovans (Denisova 2, 4 and 8), ranges of 30,000 to 300,000 years were used. For each individual, we used a uniform prior over the allowed range of dates. Neanderthals, modern humans, and Denisovans were constrained to monophyletic groups and the time to their most recent common ancestor (TMRCA) was estimated for each group.
Both the strict and relaxed clocks with a constant population size were found to be significantly better than using a Bayesian skyline population model (Bayes factors 56 7.8 and 6.08 for strict and relaxed clocks, respectively). No significant difference was found between the strict and relaxed clocks with a constant population (Bayes factor 0.76, Supplementary Table 9), therefore the combination of a strict clock and a constant population size was used as it is the simplest model. Three Markov chain Monte Carlo (MCMC) runs of 75,000,000 iterations, with a pre burn-in of 10,000,000 iterations and sampling every 2,000 trees were then performed. The log and tree files of the runs were then merged using logcombiner2 from BEAST2 54 . The merged tree file was subsequently annotated using the program treeannotater from Supplementary Figure 16 | The phylogenetic mtDNA tree determined from Bayesian analysis with Beast2 using approximately 11.3 kb of protein-coding genes from the mitochondrial genome. Each node shows the corresponding posterior probability of the branch and the x-axis represents the time in years from the present. The chimpanzee branch used to root the tree is not shown.

SUPPLEMENTARY SECTION 8: MTDNA HAPLOTYPE IDENTIFICATION FROM SPARSE DATA USING KALLISTO
As described in Supplementary Section 5, many sediment libraries contain low coverage of the mtDNA genome, or appear to contain mtDNA from multiple haplotypes, limiting the ability to construct consensus sequences. In this study, DNA from 244 sub-samples could be assigned to either the modern human, Neanderthal or Denisovan lineage via diagnostic sites (Supplementary Section 4); of which only 7 were suitable for constructing a consensus sequence of the mtDNA genome (Supplementary Section 4). Here, we use a method based on the program kallisto 62 to assign a library (or its major component in case of mixtures) to specific branches in the hominin mtDNA tree, extending the categorization of these samples beyond the identification of broad hominin groups.
We first constructed a multiple sequence alignment (MSA) for the full length of 95 mtDNA genomes, comprised of: 25 Neanderthals (including Denisova 11), one Neanderthal consensus sequence inferred from Denisova cave sediment sample M65 (Supplementary Section 6), 4 Denisovans, Sima de los Huesos, 54 modern humans and 10 ancient modern humans, as described in Supplementary Section 6 (full list of genomes in Supplementary Data File 3). To aid in analysis, we have partitioned the tree into sets of "major groups" that are closely related to each other (Supplementary Figure 17, node colours; Supplementary Figure 18, red boxes). Each of these groups is composed of mtDNA genomes which form a clade, or diverge from the larger tree at around the same time. When constructing kallisto references out of tip mtDNA genomes, we first remove all columns from the MSA in which any genome has an unknown base (i.e., an "N").
We next simulated 5,000 ancient DNA reads from each of the Sima, Neanderthal or Denisovan mtDNA genomes, as described in Supplementary Section 9. In order to simulate DNA from the entire mtDNA genome -that is, to avoid removing N columns -we inferred ancestral states across the full phylogeny (Supplementary Figure 17) with the software treetime 63 , using a TN93 mutation model and a mutation rate of 1 x 10 -10 bases per year. These ancestral states then resolve any unknown bases in the tip genomes.
We then tested the ability of kallisto to correctly place simulated DNA reads. Using the full set of mtDNA genomes as a reference set, we find that the correct genome is always the most abundant hit, receiving approximately 80-95% of the total estimated abundance (Supplementary Figure 18). The only exception is when multiple genomes are identical, in which case the abundance is equally distributed between them (e.g., Goyet individuals).
A more realistic scenario is one where the correct genome is now known -i.e., because the aDNA originates from a new sediment or skeletal sample, and is a novel mtDNA haplotype. We therefore repeated this analysis, each time dropping the source genome from the set of reference genomes. In almost all cases, the reference genome with the highest abundance is closely related to the source genome, and falls within the same major group (Supplementary Figure 19). The primary exceptions are HST and Sima de los Huesos, which have no closely related genomes -in these cases, the abundance falls on a distantly related genome. Some mtDNA references seem to "attract" abundance from distantly related genomes (e.g. Chagyrskaya 8 often has low levels of abundance -up to 21% of all non-modern human abundance), and thus may be more challenging to classify.

Probabilistic genome assignment with ancestralized mitochondrial genomes
We next used the method from Vernot et al. 64 to assign each sample to the mtDNA phylogeny (Supplementary Figure 17). In short, we simulated DNA from ancestralized genomes spanning the phylogeny, and for each reference we calculated abundance thresholds at which there is a >90%, 95% or 99% probability that the DNA originated from a closely related genome. Here, "closely related" is defined as on the same branch, or within 50,000 years of the root of the branch. These thresholds vary based on the reference (Supplementary Figure 20), with some references (e.g. Chagyrskaya 8, discussed above) requiring larger proportions of the total abundance for there to be considered strong evidence of a closely related genome. These thresholds vary depending on the number of reads used in the analysis, with larger amounts of reads generally requiring lower abundance to reach given certainty. We therefore calculated thresholds for each reference using 100 We then calculated the number of ancient hominin reads for each sub-sample as (1 -modern human contamination) * (number of unique reads). The modern human contamination proportion was calculated as described in Supplementary Section 5, using diagnostic positions. Finally, we ran kallisto on all sediment sub-samples (Supplementary Figure 21) with >100 ancient hominin reads. For each sample we calculated the proportion of non-modern human abundance assigned to each reference, under the assumption that any modern human signal largely originates from contamination. To convert these proportions in to probabilities, we applied the thresholds as calculated above. For samples with 250-999 ancient hominin reads, we conservatively use the 250 read thresholds. The results are presented in Supplementary Figure  22 and shown within the stratigraphy in Main Figure 1 for East and Main Chambers and Extended Figure 3 for South Chamber. Normalized abundance for each of 10 reference genomes (columns), for simulated ancient DNA from ancestralized genomes spanning the mitochondrial tree. Each point is one simulated mtDNA genome, and its abundance for that particular target. Red vertical line at 50,000 years divergence -points to the left of this line are from ancestralized genomes that are less than 50,000 years diverged from the reference genome branch, and are considered "closely related". For each reference, three abundance thresholds X were calculated (red dotted lines), such that 90%, 95% and 99% (bottom to top) of all simulated sequences with at least X abundance are more closely related than 50,000 years. Thresholds shown for 100-1000 simulated reads (rows). The sub-sample of SP6720 (M65) from which the consensus haplotype was inferred is denoted with a red box.

SUPPLEMENTARY SECTION 9: IDENTIFICATION OF PREVIOUSLY UNKNOWN HOMININ MTDNA LINEAGES
Even though diagnostic positions have to be defined using known mtDNA sequences (which are also required for analyses based on kallisto), archaic hominin mtDNA lineages that are yet unknown can in principle also be detected by identifying sub-samples with relatively high numbers of hominin deaminated fragments that could not be assigned to any of the known lineages. To investigate whether there is evidence for the presence of such unknown lineages in our data, we plotted the distribution of the number of deaminated fragments in sub-samples that were identified as containing ancient hominin mtDNA, but which were not assigned to any lineage (n = 51). All except two sub-samples contained less than 100 putatively deaminated fragments, compatible with that their assignment to a hominin lineage was prevented by limited statistical power. However, two sub-samples (lysate IDs: Lys565 and Lys10333), both originating from sample M76, contained 246 and 770 deaminated fragments (996 and 2646 unique hominid sequences), respectively, more than any of the other unassigned sub-samples (Supplementary Figure 23).
A closer examination of the support of diagnostic positions observed for these two sub-samples revealed significantly more than 10% support (based on 95% binomial confidence intervals) for the branch shared by all Neanderthals (i.e., the shared Neanderthal-Hohlenstein-Stadel branch, subsequently referred to as 'N-HST' branch; 31-33% support), but not for the HST branch (0.5-1.4% support) or the typical Neanderthal branch (8-13% support) (Supplementary Figure 24). Support for the N-HST branch increases to between 39% and 45% when using deaminated fragments only, while support for the HST branch and typical Neanderthal branch remains low (0% and 13-18%, respectively), suggesting the presence of Neanderthal mtDNA that shares only few derived sites with known Neanderthal mtDNA genomes. In addition to Neanderthal mtDNA, both sub-samples show significantly more than 10% support for the modern human branch (40-54%), which reduces to insignificant levels after filtering for deaminated fragments (3-14%), indicating contamination with present-day human DNA. Both sub-samples also show low but insignificant levels of support for the Denisovan branch, including after filtering for deaminated fragments (6-9% and 6-15%, respectively).
To investigate if the observed difference in branch support between the internal (N-HST) and external (N) Neanderthal branches is unique to these sub-samples, we examined all sub-samples presenting significant support for the internal (N-HST) branch and containing at least 100 deaminated fragments and tested whether there is a significant difference in support between the two branches. In order to maximize the amount of available data this comparison was performed using all fragments. Of the 37 sub-samples where this test could be performed, only the aforementioned sub-samples originating from sample M76 showed a significant difference of support between the two branches (Fisher exact test with a BH correction 65 for multiple testing, Lys565 p-value = 1.5E-2 and Lys10333 p-value = 1.2E-03). Using the percent support for the modern human mtDNA lineage as an estimation of contamination, these 37 sub-samples contained a range of present-day human contamination from 5% to 88% with a mean of 29%. The sub-samples from M76 contained an estimated 40 and 54% of human mtDNA contamination, indicating that this signal is not due to modern human contamination.
In order to test if the presence of both Neanderthal and Denisovan mtDNA in one sample (as likely is the case for M76) could lead to an imbalance in support between the internal (N-HST) and external (N) Neanderthal branches, we also examined 10 sub-samples that contained mtDNA from both groups. These sub-samples contained between 27 and 403 deaminated fragments (114 and 2,150 unique hominid sequences) and 0-52% present-day human mtDNA contamination. Based on a fisher exact test no libraries showed a significant difference in support between the two branches, indicating the presence of both Neanderthal and Denisovan mtDNA in a sample does not result in an imbalance of support between the external and internal branches.
We next investigated whether the increased support for the N-HST branch displayed by sample M76 is consistent with a mitochondrial genome sequence falling basal to the known mtDNA diversity of known Neanderthals. To do this, we simulated six ancestralized mtDNA sequences based on the Altai Neanderthal (Denisova 5) mtDNA genome, going back to six time points between 40,000 and 200,000 years before the Altai Neanderthal lived (approximately 130 ka). This was done using a tree of published mitochondrial genomes constructed from Bayesian analysis with Beast2 54 and the program treetime 63 , following the procedure outlined in Supplementary Section 8. Briefly, the Tamura-Nei 1993 (TN93) 58 mutation model was used with a mutation rate of 1e-10 mutations per base pair per year. The tree built using Beast2 estimated the time between the tip point of the Altai Neanderthal mtDNA genome to the shared node between all non-HST Neanderthals at 39,700 years and to the shared node between all Neanderthals, including HST, at 145,000 years (Supplementary Figure 25). The output of this program was a single fasta file which was then used to generate different numbers of simulated reads (100, 300, 500, 1000 and 2,000) mimicking ancient DNA fragments, using a read length distribution observed in a sample previously enriched for human mitochondrial DNA via hybridization capture (limits of 30 bp and 100 bp) and deamination rates as used in ref 61 . The simulated fragments were then processed as described in the Methods section.
The log ratio of support between the N and NHST branch was calculated for each set of ancestralized fragments and compared to the log ratio of the combined fragments from both sub-samples Lys565 and Lys10333 (Supplementary Figure 26). The difference in support between the shared N-HST branch and N branch for the sub-sample is consistent with a genome diverging from the N branch 100,000-125,000 years before Altai lived, i.e., between 20,000 and 45,000 years after the divergence of HST from other Neanderthals (255-230 ka). Sample M76 is from layer 20 of Main Chamber, which is dated to 250-170 ka 1

Identification of mammalian taxa
Of the 728 samples test, 685 contained traces of ancient DNA pertaining to at least one mammalian family. In our initial analysis, 22 families were identified. Whereas 12 of the identified ancient taxa were present in over 100 libraries, some taxa (n = 10) were detected in less than 10 libraries. These families included Eupleridae (Malagasy carnivores), which are not expected for the region and time period. We visually inspected the sequences assigned to the 10 least abundant families and found that they often concentrated in short, highly conserved parts of the mtDNA genome. This contradicts the expectation that recovered DNA fragments would originate from random parts of the mtDNA genome, and thus be distributed randomly across the reference genome when assigned correctly to a taxon of origin. We therefore excluded assignments to families for each library that were based on alignments to less than 500 bp of the mtDNA reference genome. Using this filter, we removed assignments to the Malagasy carnivores. Summary statistics for each sample and their associated negative controls are provided in Supplementary Data File 3. The relative locations of samples identified as containing ancient mammalian mtDNA, the number of DNA fragments recovered as well as their taxonomic assignments are shown in Extended Data Figures 8 and 9 for Main and East Chamber and Supplementary Figure 27 for South Chamber.

Comparison to zooarchaeological records
In depth evaluations of the zooarchaeological record have been published throughout the excavations of Denisova Cave for Main and East Chambers 1,66-71 . In comparing the proportion of mtDNA fragments assigned to taxa identified as ancient to the published zooarchaeological records we observed that small mammals (e.g. squirrels (Sciuridae), rodents (Spalacidae), etc.), are poorly represented by the genetic data and are often entirely absent (Extended Data Figure 7a,b). In contrast, all large mammals represented by skeletal remains were detected through genetic analysis, even though there are differences in the proportions of skeletal remains and DNA fragments recovered (Extended Data Figure 7c,d). One library contained ancient DNA from Camelidae, which has not been identified in the skeletal record, but is known to have been present in the region during the middle and late Pleistocene 72 . Comparisons on a per-layer basis, which were performed only for large mammals (excluding hominids) that were observed in more than one library, largely show a good congruence between the genetic and zooarchaeological data. For example, as one moves towards lower layers in the stratigraphy, the concentration of bovids decreases in Main and East Chambers in both types of data whereas the concentration of ursids and canids increases (Extended Data Figure 7c

SUPPLEMENTARY SECTION 11: FAUNAL IDENTIFICATIONS ON THE SPECIES LEVEL
To investigate potential faunal transitions based on genetic data, it is necessary to assign DNA sequences to specific species. However, analyses beyond the family level rely on the availability of complete mitochondrial genome sequences that cover the genetic diversity of extant and extinct species within a family, as well as the grouping of mtDNA sequences from the relevant species in monophyletic clades. Based on these considerations, we selected for further analysis three mammalian families for which abundant mtDNA sequence data were recovered from multiple layers throughout the stratigraphy of Denisova Cave: ursid (bears), elephantid (elephants), and hyaenids (hyaenas) (Supplementary Figures 28 to  30, respectively).
In an attempt to assign sequences to specific species, we first identified lineage-specific diagnostic positions in the mtDNA genome, following the approach in ref. 22. This was achieved by creating multiple sequence alignments for species within ursids (bears), elephantids (elephants), and hyaenids (hyaenas) using a selection of published complete mitochondrial genomes from each family. At least two different mitochondrial genomes were used per species whenever possible, with additional genomes used for the parts of the species trees in which the remains from Denisova Cave are expected to fall (brown and cave bears for ursids and woolly mammoth for elephantids, all spotted and cave hyaenas for hyaenids) 1, [66][67][68][69] (Supplementary Data File 3). These genome sequences were then aligned to each other using MAFFT 23 . To visualize the phylogenetic relationships between the genome sequences used for identifying diagnostic positions, we used Mega X 53 to build a Neighbour-Joining tree. Missing positions in any of the genome sequences were removed and evolutionary distances were measured by counting the number of differences 73 . A bootstrap test with 500 replicates was used and the support noted for each branch 74 ( Supplementary Figures 31 to 33).
Diagnostic positions were determined by identifying positions where all mtDNA genomes in one group showed a base difference to all other groups. We note that the probe design for the mammalian mitochondrial capture consisted of 5-bp tiling across the mitochondrial genomes of 242 mammals 75 . This could lead to a bias for the increased recovery of fragments from more conserved portions of the mitochondrial genome, which may in turn negatively impact the ability to identify species based on the recovered DNA fragments. To evaluate this possibility, we compared the observed coverage across the mtDNA genome to phyloP conservation scores [39][40][41] . This was carried out by merging all sequences attributed to ursid or elephantid from all samples were these families were identified as ancient (for a total of 860,488 and 336,506 sequences, respectively) and mapping them to the reference mtDNA genome available for the phyloP value determination (panda bear and African elephant, respectively). No reference with phlyoP values was available for the hyaenid family so this analysis was not performed for these sequences. The correlation between the coverage and relative conservation at each position was then tested with a Spearman's correlation test, resulting in a positive correlation for both families (bears: rho = 0.167, p-value <2.2e-16; elephants: rho = 0.178, p-value <2.2e-16 for all sequences). The sequences for each family were then mapped to mtDNA reference genomes for species identified at Denisova Cave based on skeletal remains (cave bear for bears, NC_011112.1; woolly mammoth for elephants, NC_007596.2; spotted hyaena for hyaenas, NC_020670.1). We then visualized the variation in coverage across positions with histograms, which showed that a small proportion of positions had more than twice the mean observed coverage ( Supplementary Figures 34 to 36). To avoid the dominance of a small number of highly covered diagnostic positions in subsequent analyses, the diagnostic positions showing a sequence coverage of more than twice the mean were removed from each set of positions ( Supplementary Figures 37 to 39).
Following the identification of diagnostic positions, the number of DNA fragments that support each mitochondrial group were determined separately for each sample containing ancient ursid and elephant DNA ( Supplementary Figures 40 and 41) using alignments to the cave bear mtDNA genome (NC_011112.1) for bears, the woolly mammoth mtDNA genome (NC_007596.2) for elephants, and the spotted hyaena (NC_020670.1) for hyaenas respectively. Support for a specific species was considered significant if at least 3 diagnostic positions were covered by DNA fragments from a given sample and if significantly more than 10% (using 95% confidence intervals from a binomial test) of DNA fragments covering the diagnostic positions shared the state diagnostic for the respective species.
From the 506 samples that were identified as containing ancient bear mtDNA, 413 (82%) were assigned to at least one bear mtDNA group (Supplementary Figure 40 for Main and East Chambers). From these 413 samples 174 contained cave bear mtDNA, 299 contained brown bear mtDNA and sixty of these samples contained both cave and brown bear mtDNA. Based on optical dating of layers 1 the lowest layers of Main Chamber are older than the lowest layers sampled in East Chamber. These layers predominantly contain cave bear DNA. Moving upward in the stratigraphy for Main Chamber and in the lower layers of East Chamber both brown and cave bears were identified. In the upper layers, predominantly brown bears were identified (Extended Data Figure 4). The change to predominantly brown bears when moving up the stratigraphy is consistent with the skeletal records 1 , although it should be noted that we do not have mtDNA references for all the bear species identified in Denisova Cave.
Of the 408 samples that were identified as containing ancient elephant mtDNA 367 (90%) were assigned to at least one elephant mtDNA group (Supplementary Figure 41). All 367 of these samples contained significant support for woolly mammoth mtDNA. Fifty-one of these samples also contained significant support for the African and straight-tusked elephant mtDNA clade. Three of these samples also contained evidence for Asian elephant mtDNA, however only unique four diagnostic positions out of the 167 possible diagnostic positions for the Asian elephant mtDNA were covered in these samples. The sequences covering these diagnostic positions were also blasted 76 and showed similarities to both woolly mammoth and Asian elephant. As woolly mammoths and Asian elephants are sister clades it's possible that a lack of published mtDNA genomes of elephants from the relevant time periods and geography precludes the resolution of these samples. Due to the low level of support for the Asian elephant in each of the three samples (<30%) and few diagnostic positions covered, these assignments were removed from further analyses. While woolly mammoth mtDNA was identified throughout the stratigraphy, the presence of African elephant/straight-tusked elephant mtDNA (presumably reflecting straight-tusked elephant) was most concentrated between approximately 150 and 60 ka.
For the hyaenas, all 390 samples that could be assigned significantly to a mtDNA group were assigned to cave and spotted hyaena. To differentiate between the different mtDNA haplogroups within the cave and spotted hyaenas we identified new diagnostic positions for haplogroups A, B, C and D with the striped hyaena as an outgroup (Supplementary Figure 42). We defined diagnostic positions as described above and we removed positions that had above two standard deviations from the mean coverage of all merged ancient hyaena reads. The final number of diagnostic positions used for the identification of different mtDNA groups within the bear, elephant, and hyaena families is shown in Supplementary Figure 43.
With this second set of diagnostic positions for hyaenas, we were able to assign 363 sample of the 490 samples that contained ancient hyaena mtDNA to a mtDNA haplogroup (Supplementary Figure 44). Two hundred and twenty-two of the samples contained significant support for haplogroup A, a clade containing cave and spotted hyaenas from Africa and Eurasia. Twelve samples contained significant support for haplogroup B, a European cave hyaena clade. One hundred and fifty-eight samples contained support for haplogroup D, a clade of Asian cave hyaenas. Eight samples contained support for both haplogroups A and B, twenty-on for both haplogroups A and D, and one sample contained support for both haplogroups B and D. Haplogroup A is dominant throughout much of the stratigraphy, except for the period between 200 and 100 ka when haplogroup D was dominant.
This analysis demonstrates the ability to use sediment DNA to understand transitions between mtDNA groups within different mammalian families. It also emphasizes the importance of having relevant mtDNA reference genomes in the portions of the tree that one is interested in. We are also able to identify that there is a change in mtDNA groups for bears, elephants and hyaenas between approximately 200 and 100 ka (Extended Data Figure 10).