A refined proposal for the origin of dogs: the case study of Gnirshöhle, a Magdalenian cave site

Dogs are known to be the oldest animals domesticated by humans. Although many studies have examined wolf domestication, the geographic and temporal origin of this process is still being debated. To address this issue, our study sheds new light on the early stages of wolf domestication during the Magdalenian period (16–14 ka cal BP) in the Hegau Jura region (Southwestern Germany and Switzerland). By combining morphology, genetics, and isotopes, our multidisciplinary approach helps to evaluate alternate processes driving the early phases of domestication. The isotope analysis uncovered a restricted, low δ15N protein diet for all analyzed Gnirshöhle specimens, while morphological examinations and phylogenetic relationships did not unequivocally assign them to one or the other canid lineage. Intriguingly, the newly generated mitochondrial canid genomes span the entire genetic diversity of modern dogs and wolves. Such high mitochondrial diversity could imply that Magdalenian people tamed and reared animals originating from different wolf lineages. We discuss our results in light of three ecological hypotheses and conclude that both domestication and the existence of a specialized wolf ecomorph are highly probable. However, due to their proximity to humans and a restricted diet, we propose domestication as the most likely scenario explaining the patterns observed herein.

In contrast, two pits (32 and 34) differed from that pattern with a corpus of 500 and 1,000 remains, respectively 46,47 . Aurochs (Bos primigenius) and roe deer (Capreolus capreolus) are the prominent species represented, followed by red deer (Cervus elaphus) and wild boar (Sus scrofa). Carnivores are less frequent with the occurrence of fox (Vulpes vulpes), wolf (Canis lupus), marten (Martes martes), and wild cat (Felis silvestris). Some bird bones belonged to the Anatidae family. No fish bones have been found despite the presence of a river close to the site.
In Pit 41, a portion of an axial skeleton of a canid, composed of eleven vertebrae in articulation, was found under a small sandstone slab. One cervical vertebra from a canid, with excellent preservation, has been found in Pit 57, under stones used for a post hole. Some bones of undetermined canids were also collected from Pit 34.

Frankfurt (Roman and Early Medieval)
A dog skull (inventory number α19496) comes from the Roman NIDA, modern day Frankfurt am Main -Heddernheim. The upper left M 1 was sampled from this particular dog remain.
Contextually, the dog's skull was found in the filling of a cistern in area 39 at a depth of 155-355 cm, which was placed in the IIB -III period 48 and is dated to around the end of the 3rd century AD.
Within the Franconian burial ground in Frankfurt am Main -Nieder-Erlenbach, two dog skeletons, 1986.03.001.001 and 1986.03.001.002 (NER 13/Grave 21), were retrieved from dog burials dating to the Merovingian period. A total of seven dog burials have been found within the burial ground 49 . Based on the finds and features, NER 13/Grave 21 is dated to the Young Merovingian period (JM) II, and possibly JM III 49,50 , i.e. from 640-670 AD and 670/80-720 AD, respectively 51,52 . For sampling, an upper right M 1 was used from both dogs.
Supplementary Note 2: Morphological concerns of specimens from Gnirshöhle and Hohle Fels

Susanne C. Münzel
The morphological variability of canids is considerably large. Morphological traits and metrics alone are not the determining factors when considering the ongoing debate for the domestication process [53][54][55][56][57] . A general consensus is held concerning the shortening of the snout, while the teeth remain large, resulting in tooth crowding, and in a subsequent step, the size of the teeth is reduced. However, tooth crowding alone is not a criterion for domestication either, as it has also been observed in wolves 57,58 . Conversely this means conclusions from measurements of single teeth are not appropriate to classify wolf or dog.

Gnirshöhle Mandible GN-999
Thus, we combined the length of the toothrow (ALP1M3) with the length of the M1 (CLM1) to measure the relation between snout length and length of the lower carnassial (M1) ( Figure 2).
Unfortunately, these two measurements for the mandible were only available for the Kesslerloch wolves (see Pleistocene wolves in Figure 2); the dog from Kesslerloch is represented by a maxilla. For the other dog from Bonn-Oberkassel, the length of the tooth row in the mandible is considerably shorter, because the length of the alveoli was measured from P1 to M2, although the alveolus of M3 is visible 59 . Furthermore, this mandible bears some pathologies, such as an agenesis (pathologically not erupted) of P2 and P3 contributing to the shortening of the jaw 60 .
Given this pathology, this dog mandible from Bonn-Oberkassel was not comparable with Gnirshöhle (GN-999).

Hohle Fels Maxilla HF-530
When Edgard Camerós 63 brought the wolf maxilla (HF-530) from Hohle Fels from the uppermost Gravettian layer (AH IIB) as a possible domesticate into discussion, he triggered an avalanche of debate about 'Palaeo-dogs' in the pre-LGM of the Swabian Jura 60,64 . Camerós' argument for the maxilla HF-530 as a possible domesticate was the prominent protocon on the upper P 4 resembling that of the dog maxilla from Kesslerloch 63 , but the debate about morphological traits and the discussion about teeth and mandible sizes in recent years demonstrates that with single traits and morphometrics alone, the classification is not possible if a wolf is domesticated or not [53][54][55][56][57]

Gnirshöhle
Faunal remains of Gnirshöhle had been previously dated in the 1970s; however, at that time, large sample quantities of bone were required for one radiocarbon measurement, e.g. 90g, including bones of different species, such as five reindeer, one horse and one canid bone fragment. The old dates of primarily game species range between 12-13 ka BP (15-16.5 ka cal BP) 10 . Since the old dates have a large standard deviation and were a mixture of bones from different species, it was especially important to obtain new dates for the canids.
We used five collagen samples from canids for the new radiocarbon dating, all of which were also extracted for isotopic and genetic analysis (Table S1). Preparation for the radiocarbon dating followed the protocols of Hajdas and colleagues 65,66 and dating was performed at the Laboratory of Ion Beam Physics (ETH Zurich) by measuring the 14 C/ 12 C ratio using the MICADAS accelerator mass spectrometry (AMS). Finally, the 14 C dates were calibrated by using  (Table S1). One important outcome of the radiocarbon dating is that both chambers of the cave, Gnirshöhle I and II, are contemporary and thus, the canid remains were deposited during the time range of Magdalenian occupation in the Bruder Valley between 16 and 15 ka cal BP 10 ( Figure S1 and Table S1). This period still belongs into the Pleniglacial stage GS-2a (NGRIP-oxygen isotopes). The previously obtained dates seem to be slightly older, but the standard deviation, +/-300 years, is much larger than for the new dates 10 . The newly dated canids from Gnirshöhle are ca. 1,000 to 1,500 years older than the dogs from Kesslerloch or Bonn-Oberkassel ( Figure S1).

Hohle Fels
Six wolf remains from Hohle Fels were dated, which were essential for the genetic diversity calculations. Two of the specimens (HF-530 and HF-1250.2) date to the Magdalenian and Late Palaeolithic periods, respectively (Table S1 and Figure S1). The dates of the other four canids (HF-912, HF-1965, HF-1035, and HF-1390, Table S1) align with the Gravettian and Aurignacian periods in Hohle Fels.

Some concerns with regard to the chronology in Hohle Fels
The maxilla HF-530 originates from the uppermost Gravettian layer (AH IIB), which is in direct contact to the Magdalenian layer. Therefore, pre-and post-LGM layers are sometimes difficult to distinguish in Hohle Fels. The assignment of a possible domesticate by Camerós and colleagues 63 initiated a discussion about 'Palaeo-dogs' in the pre-LGM of the Swabian Jura.
Later, Prassack and colleagues 64 listed this specimen as a Gravettian dog and in the study of Janssens and colleagues 60 , it has even been cited as a domesticate from the Aurignacian. In the light of this discussion, the sample was dated as part of this study and is of Magdalenian provenance. However, we cannot confirm the previous assignment as a dog.
Another chronological problem occurred for sample HF-912, a tibia shaft fragment found in layer GH 3AD/ AH IIAD assigned to the Magdalenian. However, during the excavation an admixture of some older fauna was recognized by a darker patina of the bones, while the typical Magdalenian fauna show a fresher 'yellow-' or 'honey-colored' patina. Sample HF-912 has a darker patina and most likely belongs to the pre-LGM period 29 , which was corroborated by radiocarbon dating (31.5 -31.2 ka cal BP, Table S1). This is also supported by the BEAST analysis ( Figure 4, Figure S6, Table S6).
The samples HF-1250.1 and HF-1250.2, a 2nd metacarpal and a carpal bone, respectively, originated from a collapsed profile in one of the uppermost layers. The carpal bone (HF-1250.2) was dated to 13.3 -13.1 ka cal BP. From an archaeozoological perspective, the two elements probably belong to one individual. This assumption is supported by the identical mitochondrial DNA sequence as well as by isotopic analysis. Accordingly, both samples were treated as one individual for isospace, and especially for all genetic downstream analysis, the sequencing data of both samples were merged. Both bone fragments could be assigned to the Late Palaeolithic period (14.2 -11.6 ka cal BP), which is slightly younger than the Magdalenian period 10,70-72 (ca. 17.5 -14.2 ka cal BP).

Detailed methods for elemental and isotopic analysis
For the isotopic analysis, we cut small samples (0.3 -0.7 g) from the bones using a Saeshin Forte 200 alpha micro-circular saw. After successive cleaning in purified Millipore water and acetone, we manually ground the samples to a powder (grain size less than 0.7 mm). To evaluate if collagen was still present in the samples, we performed a CN elemental analysis of the bone samples following the protocol of Bocherens and colleagues 73 . This analysis was performed at the Hydrogeochemisty department (University of Tübingen) using the Vario EL elemental analyzer. Sulfanilic acid from Merck was used as the international standard.
Collagen extraction followed the protocol of Bocherens and colleagues 74 , and was performed in the laboratory of the Biogeology working group (University of Tübingen). Depending on the nitrogen content (%Nbone) of each sample, as measured by the CNS analysis, we used 120 mg The laboratory used the international laboratory standard IAEA 600 (caffeine) as well as two inhouse reference materials (modern collagen of camel and elk). An analytical error below 0.2‰ (1σ) was determined for δ 13 C and δ 15 N in all the repeated analyses. The reproducibility error for the amounts of C and N was lower than 2%.
Following the protocols of DeNiro 75 and Ambrose 76 , we only used collagen samples with a carbon-to-nitrogen-ratio (C:Ncoll) between 2.9 and 3.6 and a nitrogen percentage higher than 5% for palaeoecological interpretations.

Detailed methods of niche modeling and dietary reconstruction
To reconstruct the niches of the canids, we first applied a multivariate cluster analysis to the δ 13 C and δ 15 N values in JMP 14. As a result, we obtained three different clusters that were not dependent on species determination; instead, they depended exclusively on the individual trophic source reflected by carbon and nitrogen isotopes. We then used SIBER (Stable Isotope Bayesian Ellipses in R) to calibrate the niches out of the clusters 77 . With this R package, it was possible to reconstruct the complete niches (= convex hull 78 ) and the core niches (= standard ellipse area 77 ) of the canid communities. The complete niche includes all individuals of a niche and is quite sensitive to sample size. In contrast, the core niche depicts the center of a niche that is calculated using the maximum likelihood estimation and explains 40% of the variability 77 .
According to Jackson and colleagues 77 , this is more reliable for analyzing small sample sizes and recommended for niche interpretations. Since Magdalenian sites in Central Europe were ecologically similar to those in Southwestern Germany and Northern Switzerland 79 , we included 49 published δ 13 C and δ 15 N values from bone collagen of large and small herbivores from other parts in Central Europe to obtain a comprehensive isotopic dataset (Table S5). TEF values correspond to the difference between the stable isotope ratios of the consumer (predator collagen) and its diet (prey collagen) and are the result of the discrimination of stable isotopes due to the behavior and physiology of the consumer [82][83][84] . For this study, we used TEF values (Δ 13 C = 1.1 ± 1.1‰; Δ 15 N = 3.2 ± 1.8‰) from a study on modern fox by Krajcarz and colleagues 83 .
To get a robust statistical analysis, we set the MCMC (Markov Chain Monte Carlo 81 ) chain length to 1,000,000 with a burn-in of 500,000 in 3 chains. We verified the model convergence with Gelman-Rubin and Geweke tests, which shows model convergence, if the values are near 1, but in most analyses, values below 1.1 are acceptable 85 . Additionally, the Geweke test compares the mean of the first part of each chain with the second part's mean, using a two-sided z-test. If both means are the same, the model is convergent 81 .

Results of the elemental and isotopic analysis, and convergence test for MixSIAR
The percentage of nitrogen in bone (%Nbone) was measured in eight samples. However, only six of them confirmed the favorable conditions of preservation (1.1 -3.1% Nbone), establishing quantitatively that collagen is well-preserved in the samples. Moreover, the atomic C:Ncoll ratios of all analyzed extracts (3.3 -3.4) confirmed that collagen preservation was appropriate for the interpretation of the isotopic analysis (Table S1, Table S5). Among the isotopic values, we found a clear difference between the newly analyzed fox from Bockstein and canids from Gnirshöhle.
The red fox from Bockstein has a δ 13 C value of −20.3‰ and a δ 15  BWA seedlength (-l) of 1,000 to effectively turn off seeding, BWA Max # Diff (-n) of 0.01 allowing less differences of reads to the reference sequence, and BWA Qualityfilter of 20, to discard reads with a mapping quality score lower than 20. After read mapping, the genome coverage was estimated by QualiMap 100 integrated within the EAGER pipeline. First, a full genotyping was performed using GATK 101  Additionally, the Roman specimen (F-α19496) was positioned basally to the two modern mitochondrial DNA sequences of the dog D clade ( Figure S5). Taken together, these three ancient mitochondrial DNA sequences, APC-19, APC20, and F-α19496, positioned close to the modern dog haplogroup D specimens, were part of a sub-branch encompassing 15 specimens ranging in the time from the Magdalenian period to modern times. Moreover, we observed that the Medieval dog specimens from Frankfurt, were placed within the dog A clade.
In the phylogenetic tree, the mitochondrial genome of HF-1250 is closely located to the Magdalenian mitochondrial genome of the GN-999 specimen ( Figure S5). The four Canadian mitochondrial genomes clustered together and were positioned closely and basally to a cluster of modern wolf mitochondrial DNA sequences originating from North America ( Figure S5).

Discussion of the Phylogenetic Results of the non-Magdalenian canids
While the main focus of this study concentrates on the interpretation of the new southwestern German mitochondrial genomes, it was also necessary to assemble a more comprehensive dataset with which those samples could be more precisely analyzed phylogenetically. This Lastly, contrary to previous findings 93 , and mainly due to the fact that we included dog mitochondrial genomes in our phylogenetic analyses, we were able to more precisely assign the two samples KSL-189 and Siberia5. With regard to the KSL-189, we can provide genetic evidence that this sample is more closely associated with dogs rather than with wolves , which matches previous morphological evidence 61 . The same holds true for Siberia5, a sample previously considered as a wolf but herein receiving genetic support to potentially resemble a dog.        Table S1: Archaeological context of the canid samples in chronological order (GH=geological horizon; AH=archaeological horizon; (1) bone comes from a collapsed profile, after dating it was assigned to the Late Palaeolithic; (2) the maxilla is from the uppermost Gravettian layer AH IIb, which is in direct contact with the Magdalenian; the new date corrected its affiliation into the Magdalenian; (3) these Umingmak samples were not directly dated, the time range is given by earlier radiocarbon dates 35