Neandertal cannibalism and Neandertal bones used as tools in Northern Europe

Almost 150 years after the first identification of Neandertal skeletal material, the cognitive and symbolic abilities of these populations remain a subject of intense debate. We present 99 new Neandertal remains from the Troisième caverne of Goyet (Belgium) dated to 40,500–45,500 calBP. The remains were identified through a multidisciplinary study that combines morphometrics, taphonomy, stable isotopes, radiocarbon dating and genetic analyses. The Goyet Neandertal bones show distinctive anthropogenic modifications, which provides clear evidence for butchery activities as well as four bones having been used for retouching stone tools. In addition to being the first site to have yielded multiple Neandertal bones used as retouchers, Goyet not only provides the first unambiguous evidence of Neandertal cannibalism in Northern Europe, but also highlights considerable diversity in mortuary behaviour among the region’s late Neandertal population in the period immediately preceding their disappearance.

Supplementary Fig. S2. The 21 drawers of fragmentary, "indeterminate" fauna from Dupont's excavations at the Troisième caverne of Goyet that were systematically sorted in order to identify any overlooked human remains.
The comparative taphonomic analysis of the Goyet Neandertal remains was conducted on the faunal remains identified in drawers Q53, Q55, Q375, and Q376 (see Supplementary Table S5).

Supplementary Fig. S3. Left and right tibia pieces and lower left lateral incisor used to estimate the Neandertal MNI for Goyet.
Anterior (left) and posterior (right) view of each tibia piece; lingual (left) and distal (right) view of the lower left lateral incisor. Tibias I and II and Q375-2 are from the left side; all of the other tibia pieces are from the right side. Scale = 3 cm for the bone pieces and 1 cm for the tooth. Left: individual specimens; right: refit piece (from left to right: in medial, anterior, lateral and posterior views). Note that all of the specimens were found mixed with fauna from E. Dupont excavations and a small yellow label indicating their stratigraphic origin was glued to each at the beginning of the 20th century; red traces were also drawn on likely faunal fragments to delimit impact notches and retoucher areas. Scale = 5 cm.   Numbers at the main branch nodes represent bootstrap values after 1,000 iterations. Supplementary Table S1. Human remains from the Troisième caverne of Goyet identified as Neandertal with indication of the "fauna-bearing level" (FBL) and analyses performed. See Supplementary Fig. S8 for the placement of the anthropogenic modifications on the bones.  (41) were also performed. The PR is calculated as MNE*100/(MNImax*NEind) with MNImax being the highest MNI for the whole sample and NEind the number of elements per individual.

Supplementary Notes
Supplementary Note S1. The Troisième caverne of Goyet and its regional context The Goyet caves are located in Mozet, Belgium, some 20 km from the well-known site of Spy ( Supplementary Fig. S1). The Troisième caverne of Goyet (50°26'44"N, 5°00'48"E) is part of a   Fig. S1). In Germany, the remains of at least three individuals identified at the Neandertal type-site have been attributed to this group along with a partial parietal bone from Warendorf (7). While several other German sites have yielded Middle Palaeolithic human remains, their age is still too uncertain to be considered here (7). This is also the case with an isolated frontal bone fragment recovered from the North Sea (8). Further east, a small number of isolated teeth have been recently identified at Stajnia Cave (Poland; 9). In Belgium, the Trou de l'Abîme at Couvin and Walou Cave at Trooz have each yielded an isolated tooth, whereas the remains of two adults and one juvenile were discovered at Spy (10).

Supplementary Note S2. Assessment of the Goyet Middle Palaeolithic material and overview of the Late Mousterian from the Mosan Basin
The lack of field data from Dupont's excavations at the Troisième caverne makes it impossible to determine whether the typo-technologically Middle Palaeolithic material corresponds to one or several occupations. The latter possibility was suggested in the only, albeit incomplete study of this material by Ulrix-Closset (11), basing her conclusion on the diversity of surface alterations and several typological arguments.
Ulrix-Closset noted that Levallois technology is poorly represented, with débitage involving primarily "spherical cores." Various scrapers forms are the most abundant retouched tool type, followed by 45 Mousterian points and the occasional limace. She also stressed the presence of bifacial scrapers alongside 28, often 'atypical' bifacial tools, three of which resemble foliate pieces (11). The numerous denticulates and raclettes present in the assemblage are more appropriately interpreted as edge-damaged artefacts, highlighting the significant post-depositional reworking of the deposits. Ulrix-Closset (11) (14). The material from Trou de l'Abîme at Couvin is made on a fine-grained, non-local flint that was heavily reduced.
Scrapers are the most well-represented tool type and occur alongside small bifacial pieces (15).  Table S1). Moreover, the latter shows a very different aspect from the other human remains labelled 2861, which are fresh and partially eroded, pointing to a different taphonomic history. We believe that maxilla 2861-1 may have been misplaced with the 2861 material after the excavations but before the remains numbered 2878 were inventoried and labelled at the end of the 19th century.

Supplementary Note S5. The Neandertal assemblage from the Troisième caverne
Identification. The Neandertal remains from the Troisième caverne were isolated from the rest of the human sample on the basis of their morphometric characteristics combined with their taphonomic aspect, isotopic ratios, radiocarbon dating and genetic analysis (see Table 1 and Supplementary Tables S1 and S3). -The most complete tibia of the Neandertal sample, left tibia I, cannot be the antimere of tibias III, V or VI as its mtDNA sequence differs from those obtained for these bones nor can it be associated with left tibia II as they both preserve the area of the tibial tuberosity in addition to carrying different mtDNA sequences. Unfortunately, the mtDNA of tibia IV is not sufficiently preserved to determine whether it carried the same mtDNA sequence as that of tibia I. The morphometric characteristics and taphonomic aspect of the two bones are, however, compatible with their being corresponding antimeres, and we cannot exclude that they belong to the same individual. with femur II whose morphometric and taphonomic characteristics are compatible with it belonging to the same individual, hence we tentatively associated them. Femur I shares an identical mtDNA sequence with right tibias III and V, and its morphometric and taphonomic characteristics are compatible with it belonging to the same individual as one of these tibias. Finally, we tentatively associated left and right tibias I and IV based on morphometric and taphonomic similarities (see above).

Supplementary Note S6. Radiocarbon dating of the Goyet Neandertals
Radiocarbon dates are reported in years before present (BP) following the convention proposed by Mook and van der Plicht (20). These dates require calibration to obtain calendar ages. The presently recommended calibration curve is IntCal13 (21)

Supplementary Note S7. Palaeogenetic analyses of the Goyet Neandertals
Each bone fragment selected for palaeogenetic analysis was first irradiated with UV light in order to reduce surface DNA contamination. A dental drill was used to remove a thin layer of bone surface and to sample inside the bone. An aliquot of between 30 mg and 120 mg of bone powder was utilized in the DNA extraction following an optimized protocol to retrieve typical short ancient DNA (aDNA) fragments (23). 10-20 µl out of 100 µl of extract were transformed into a sequencing library using a double stranded library preparation protocol (24) and indexed with an individual double index combination (25). Different amplification cycles were used for each indexed sequencing library in order to avoid heteroduplex products formation. Mitochondrial DNA (mtDNA) was subsequently enriched through a bead-capture protocol that uses modern human mtDNA fragments as baits (26). The enriched libraries were re-amplified, quantified on a DNA 1000 chip (Agilent), pooled in equal concentration with other samples, and sequenced on an Illumina HiSeq 2500 rapid run (2x101+8+8 cycles).
The sequenced reads were quality filtered and merged using established protocols (27). Only merged reads with a length above 30 base pairs (bp) were mapped against a Neandertal mitochondrial reference sequence using BWA (28) in combination with an in-house developed circular mapping tool (29).
After removal of identical reads appearing more than once, sequences with a mapping quality score lower than 30 were excluded using the SAMtools software package (30). Only unique sequences securely placed within the mtDNA (Supplementary Table S7) were used to reconstruct the mitochondrial consensus sequence of each sample for positions with at least 5-fold coverage using the custom iterative assembler MIA (31).
In order to authenticate taxonomic assignment, reads of the three low coverage samples (C5-1, Q55-4 and Q119-2) were mapped against the modern human mitochondrial reference sequence (rCRS).
The mapping results are in the same range as the values obtained using the Neandertal reference (Supplementary Table S8) therefore excluding a reference bias in the taxonomic assignment. Moreover, potential contamination with modern human DNA was assessed with a contamination estimation software that considers positions where the Neandertal mtDNA reference sequence differs from at least 99 % of 311 worldwide modern human mtDNAs (31). For each of these diagnostic positions, we compared the number of sequences matching the Neandertal reference better (clean fragments) than modern human mtDNAs (polluting fragments) to calculate contamination (Supplementary Table S7). Moreover, the characteristic damage pattern of aDNA was calculated as the percentage of reads showing C to T or G to A misincorporations respectively at the 5' and 3' ends of DNA fragments (Supplementary Table S7 and Supplementary Fig. S17) using a program first used in Briggs et al. (32).
The phylogenetic placement of the seven newly generated complete or almost complete mitochondrial sequences (i.e. at least 98 % complete) was assessed by comparing them to modern human, Neandertal and Denisovan mitochondrial genomes. The MUSCLE software (33) was used to align mtDNA consensus sequences of the Goyet specimens, 54 modern humans belonging to different worldwide language groups (34), eight Neandertals (31,(35)(36)(37), and one Denisovan individual (38). A maximum parsimony tree and a maximum likelihood tree with complete deletions (16,110 positions considered) and 1,000 iterations as bootstrap support were built using MEGA 5.2 (39) and refined with FigTree (ref. 40; Fig. 2 and Supplementary Fig. S4).
Subsequently, reads with nucleotide misincorporations (postmortem damage or PMD score ≥ 3) indicating authentic ancient origin (Supplementary Table S7) were selected with PMD tools (37). The percentage of filtered fragments with damaged termini increased up to 71 % whereas contamination decreased for all samples (Supplementary Table S7). Only the filtered reads of the Goyet samples represented in Fig. 2 were used to generate new mitochondrial consensus sequences (for positions covered at least 5 times) that were aligned to the same Neandertal mtDNA reference and assembled in additional maximum parsimony and maximum likelihood trees (Supplementary Figs. S18 and S19) using the same parameters mentioned above (10,234 positions considered). This confirmed the phylogenetic placement of the original reconstructed mitochondrial consensus sequences within the Neandertal mtDNA diversity and validated the intragroup matrilineal relationships.