Palaeogenomics suggest domesticated camelid herding and wild camelid hunting in early pastoralist societies in the Atacama Desert

13 The domestication of camelids played a significant role in the development of Andean 14 societies. Although South American camelids can be osteometrically classified as large (putatively 15 Lama), or small (putatively Vicugna), it is difficult to differentiate between wild or domesticated. 16 Here we utilise palaeogenomic data to reveal the sex and ancestry of camelids from archaeological 17 sites from the Tulán ravine, a key area for understanding the transition from hunter-gatherer to 18 pastoral societies in the South-Central Andes. We inferred the ancestry of 49 individuals with 19 genome-wide coverages >0.001x and found evidence for both Vicugna and Lama. Investigations 20 into 26 individuals of >0.01x showed that all individuals, except one male llama, represented 21 ancestry not found in modern individuals. Similar male-to-female sex ratios suggest hunting rather 22 than herding. Moreover, while intergeneric admixture is widespread in modern domesticated 23 individuals, we find limited evidence of this in our dataset. Overall, our findings suggest that lost 24 early domesticate lineages and/or wild camelids, not the direct ancestors of contemporary 25 domesticates, were the primary source of camelids in both domestic and ritual contexts in the Tulán 26 ravine during the Early Formative period (3,360–2,370 cal. yr BP). 27 28

29 One of the most significant cultural developments of humanity was the domestication of 30 plants and animals 1. Domestication refers to the process by which humans selectively breed wild 31 animals and plants for specific traits that serve human needs, leading to genetic, morphological, and 32 behavioral changes over generations 2. As domesticated animals became central to human societies, 33 pastoralism emerged as a subsistence strategy, where humans raise and manage these animals for 34 resources such as food, fiber, and transport. However, despite its remarkable significance, 35 domestication processes (e.g. wild ancestors, timing, location, and underlying genetic mechanisms) 36 remain unsolved for many species 3. This is especially true for the only large herd animals 37 domesticated in the Americas, South American camelids 1. The domestication of South American 38 camelids played a significant role in the development of Andean societies. Their lone domesticated 39 status and ability to occupy broad domesticate roles (transport, fibre, food) exemplifies the pivotal 40 role these animals played in shaping ancient South American societies. Camelid fibres and meat 41 served both practical and ritual purposes, while pack llama (Lama glama) acted to connect distant 42 communities 4–6. Even today, South American camelids remain the most important resource for 43 many communities 7. 44 45 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint South American camelids comprise four extant species; two domesticated (llama [Lama 46 glama], alpaca [Vicugna pacos]), and two wild (guanaco [L.guanaco], vicuña [V. vicugna]). Within 47 the wild species there are two subspecies each: L.g.cacsilensis, L.g.guanicoe, V.v.mensalis, and 48 V.v.vicugna. While multiple regions show evidence of early camelid management, such as hunting, 49 herd protection, and the early stages of captive care, the most likely centers of domestication, where 50 the controlled breeding of camelids began, are currently considered to be the Central Andes, 51 Southern Bolivia, and the Circumpuna region spanning Northern Chile and NW Argentina8–14. 52 However, the timing and number of domestication centres remain debated 15. Archaeological 53 evidence suggests South American camelid domestication was a complex process and not 54 synchronous across the Andes that began with hunting, herd protection, and finally 55 captivity/selective breeding 16. However, although South American camelids can be osteometrically 56 classified as large (putatively Lama), or small (putatively Vicugna), it is difficult to designate an 57 individual as wild or domesticated 17,18. This is primarily due to the significant overlap in size 58 between wild and domesticated forms, as well as the lack of clear, consistent morphological 59 differences between them. These factors complicate efforts to definitively classify individuals as 60 wild or domesticated, limiting insights into the presence of early domesticated individuals. 61 62 Genetics is a useful tool that can complement the archaeological record to provide further 63 insights into domestication 19. However, despite attempts to address theories of South American 64 camelid domestication, the Spanish conquest of South America beginning ~1500 A.D and 65 uncontrolled post-conquest intergeneric hybridisation between domesticated species makes it 66 difficult to identify genetic affinities prior to ~1500 A.D from contemporary individuals 20. This 67 intergeneric hybridization, which is thought to have become widespread following the Spanish 68 arrival, has significantly shaped the genetic makeup of modern domesticated llama and alpaca 69 populations 20. As a result, modern populations are highly admixed, complicating efforts to trace 70 pure genetic lineages that could shed light on the early domestication process. Ancient DNA 71 (aDNA) from pre-conquest individuals is therefore a useful approach to reconstructing the 72 uncertain past of South American camelid domestication and husbandry. However, although 73 previous studies have utilised pre-conquest individuals, these have been limited to mitochondrial 74 DNA 17,21–23. Mitogenomes are maternally inherited, only represent a single locus, and are 75 confounded by hybridisation and stochastic events 24. This can be especially problematic in the case 76 of South American camelids as modern individuals are known to contain high levels of 77 hybridisation 20, limiting the power of mitochondrial reference panels for taxonomic identifications. 78 79 Here we utilise palaeogenomic data to investigate the use of camelids by ancient South 80 American communities. We analyse the sex and ancestry from pre-conquest South American 81 camelids recovered from archaeological sites from the Tulán ravine (Tulán-52, Tulán-54, Tulán-85, 82 and Tulán-94) located in the eastern border of the Salar de Atacama that descends into Northern 83 Chile (Fig 1). The Tulán ravine is a key area for the transition from hunter-gatherer societies, which 84 are characterised by subsistence through hunting, fishing, and foraging of wild plants, to initial 85 pastoral societies, where domesticated animals are herded for food, labour, and other resources, in 86 the South-Central Andes 25. The sites studied here belong to the Late Archaic period (Tulán-52; 87 5290-4840/4430-4090 cal. yBP) and to the Early Formative period (Tulán-54, Tulán-85, and Tulán-88 94; 3,360–2,370 cal. yr BP) when domesticated camelids are thought to have been incorporated into 89 the economic activities of humans of this area. Specifically, archaeological remains from ritual 90 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint (sites associated with ceremonial or religious practices) and domestic (sites used for daily living 91 and subsistence activities) contexts provide evidence of a significant cultural and economic 92 transformation from in the highlands of the Southern Andes. 26. 93 94

95 96 Mapping results 97 Our study utilised genomic data from a total of 75 ancient camelid individuals, sourced from 98 four archaeological sites: Tulán-52, Tulán-54, Tulán-85, and Tulán-94. These individuals had 99 previously undergone mitochondrial and morphological analysis (Supplementary Tables S1 and S2) 100 but have yet to be analysed from a nuclear genomics perspective. After mapping to the alpaca 101

genome, genome-wide coverages ranged between ~0 and 0.28x or 4,500 and 587 million 102 mapped base pairs (Supplementary table S2). The ancient samples showed typical ancient DNA 103 damage patterns with elevated levels of A-G and C-T transitions towards the read ends 104 (Supplementary Fig S1). 105 To build a genomic reference panel, we also incorporated previously published genomic 106 data from 80 modern individuals representing all extant Lama and Vicugna subspecies, which 107 served as a comparative basis for the ancient samples. After mapping to the alpaca reference 108 genome, all modern individuals had genome-wide coverages ranging from 8.34x to 41.11x 109 (Supplementary table S1). 110 111 Principal Component Analysis 112 To evaluate the broad relationships between the modern and ancient camelids included in 113 this study, we performed principle component analyses (PCAs) using a pseudohaploid base call 114 with ANGSD. When analysing the entire dataset (Lama sp. and Vicugna sp.), the modern 115 individuals fell into five distinct clusters, in most cases pertaining to their subspecies classification 116 (Fig 2). However, we found L.g.cacsilensis and llamas clustering together, apart from two 117 L.g.cacsilensis which clustered with L.g.guanicoe (Cacsilensis1 and Cacsilensis2). These two 118 individuals are from outside of the known range of L.g.cacsilensis and were therefore likely 119 misidentified 20. 120 When considering the ancient samples >0.01x, we see three groupings; One closest to 121 V.v.vicugna and mostly containing osteometrically small size individuals, one mostly containing 122 osteometrically large size individuals somewhat in the middle of the two modern Lama sp. clusters, 123 and a single ancient individual (C72) within the llama/L.g.cacsilensis cluster (Fig 2A). When 124 reducing the dataset to only include Lama sp. individuals, we again see C72 within the 125 llama/L.g.cacsilensis cluster, with the remaining individuals still forming their own cluster, but 126 sitting more closely to the llama/L.g.cacsilensis cluster (Fig 2B). Results are similar when 127 considering ancient individuals >0.001x; the cluster closest to V.v.vicugna remains but the 128 structuring within Lama sp. becomes less clear (Supplementary Fig S2). Individuals where the 129 osteometric size did not match with their PCA placement included a single big individual (C38), 130 and two small individuals (C69 and C70). 131 132 Phylogenetic Analysis 133 To further investigate the broad relationships between the modern and ancient camelids 134 included in this study, we computed a neighbour joining tree. Our neighbour joining tree showed 135 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint similar results as the PCA (Supplementary Fig S3). There is a split in the ancient individuals, with 136 individuals being most closely related to V.v.vicugna or the llama/L.g.cacsilensis clade. Only 137 individual C72 was found to be within the llama/L.g.cacsilensis clade. 138 139 Ancestry proportions 140 To obtain high-resolution ancestry profiles, clarify the origins of each specimen, and assess 141 possible interspecific gene flow, we calculated ancestry proportions for eight ancient individuals 142 using admixfrog 27. Investigations using simulated aDNA damaged modern individuals showed a 143 good ability to discern subspecies designation in the wild individuals using this approach as the 144 highest proportion of single subspecies ancestry in each individual was as expected (Supplementary 145 Fig S4). In the llama we saw a relatively high proportion of ancestral camelid DNA (0.11 vs. 0.01 146 in the wild species) but negligible levels of Vicugna ancestry (0.01). However, in the alpaca, we 147 saw a much higher proportion of ancestral camelid DNA (0.44) as well as relatively high levels of 148 Lama ancestry (0.13). 149 The empirical ancestry proportion results were in line with the PCA and phylogenetic tree 150

(Fig 2C and D). C72 had more L.g.cacsilensis ancestry relative to L.g.guanicoe compared to 151 the other ancient Lama sp. individuals (Fig 2C). When looking into Lama sp. specific ancestry, C72 152 contained relatively higher levels of shared ancestry with modern llama than the other three 153 investigated individuals (Fig 2D). Overall, the genetically identified Lama sp. individuals had more 154 L.g.cacsilensis ancestry than L.g.guanicoe and Vicugna sp (Fig 2C). Putative Vicugna individuals 155 have more V.v.vicugna ancestry than V.v.mensalis (Fig 2C). The small sized C69 contained mostly 156 L.guanicoe ancestry, and the big size C38 contained mostly V.vicugna ancestry. Only a very small 157 proportion of each individual contained ancestral South American camelid ancestry <0.04 (Fig 2C). 158 When using the wild species reference panel, the vast majority of the oldest individual’s, 159 C46 from Tulán-52, ancestry could not be determined beyond ancestral Lama with the highest 160 proportion of single subspecies ancestry (2%) coming from L.g.cacsilensis, the next closest being 161 only 0.1% from L.g.guanicoe. When using the Lama only reference panel we found 48.9% of its 162 ancestry as ancestral Lama, 50.1% as ancestral L.g.cacsilensis/llama, 0.58% as L.g.cacsilensis only, 163 0.25% as llama only, and 0.13% as L.g.guanicoe only. 164 165 D-statistics 166 To further test the taxonomic placement of the archaeological specimens and to identify 167 any signatures of interspecific gene flow we used D-statistics. However, very low-coverage ancient 168 DNA data has been shown to introduce biases in D-statistics analyses, especially when mapping to 169 a divergent reference genome 28. To address potential biases in our data, we conducted a simulated 170 dataset analysis. Placing an individual with simulated ancient DNA damage into the H3 position 171 produced comparable results to when placing a high quality version of the same individual in the 172 H3 position (Supplementary Figure S5). Therefore, we could proceed with confidence knowing 173 ancient DNA damage would not bias our topology/population structure results. 174 D-statistics results using the topology H1=Lama sp. H2=Vicugna sp. H3=Ancient 175 individuals returned results consistent with the PCA, phylogenetic tree, and ancestry proportion 176 analysis; individuals C75, C72, C89, C21 were more closely related to Lama sp. and individuals 177 C50, C39, C38, and C20 were more closely related to Vicugna sp. (Fig 3A). Further investigations 178 into the Lama sp. individuals showed that most had an overall closer relation to L.g.cacsilensis 179 relative to llama or L.g.guanicoe (Fig 3B). However, C72 had similar levels of L.g.cacsilensis and 180 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint L.glama ancestry, as seen by few significant results when calculating D-scores using L.g.cacsilensis 181 and L.glama in the H1 and H2 position respectively. Further investigations into the Vicugna sp. 182 individuals showed an overall closer relation to V.v.vicugna relative to alpaca or V.v.mensalis (Fig 183 3C). 184 Comparing results placing a llama individual with simulated ancient DNA damage into the 185 H2 position, a Lama sp. in the H1 position, and a Vicugna sp. individual in the H3 position to those 186 when placing a high quality version of the same llama individual in the H2 position, we saw a slight 187 discrepancy in the D-scores. The mean deviation was 0.0322 which we used to correct for ancient 188 DNA related biases in the interpretations of the empirical data (Supplementary Figure S6). 189 Investigations into the topology ((Lama sp., C21), Vicugna sp.) found that individual C21 190 had less gene flow with Vicugna sp. than all modern llama individuals, similar levels of gene flow 191 with V.v.vicugna and V.v.menalisis but more with alpaca compared to L.g.guanicoe, and less gene 192 flow with alpaca but similar levels with V.v.vicugna and V.v.menalisis compared to L.g.cacsilensis 193 (Fig 4). Investigations into the topologies ((Lama sp., C72), Vicugna sp.) found that Tulan 194 individual C72 had less gene flow with Vicugna sp. than modern llama, similar levels of gene flow 195 with V.v.vicugna and V.v.menalisis but more with alpaca compared to L.g.guanicoe, and 196 comparable levels with all Vicugna sp. compared to L.g.cacsilensis (Fig 4). 197 198 Genetic sex determination 199 To evaluate whether the sex distributions in our archaeological assemblages reflect hunting 200 or herding strategies, we genetically sexed all individuals with >5,000 mapped reads using the 201 SeXY pipeline 29. We were able to genetically determine the sex of 32 males and 29 females 202 (Supplementary table S2). We were unable to determine the sex of 14 individuals due to either low 203 number of mapping reads or ambiguous results. Further separating into each site, we found: Tulán-204 52 (1 F), Tulán-54 (24 F and 26 M), Tulán-85 (6 M and 3 F), and Tulán-94 (1 F). Further splitting 205 by taxonomic status, we found similar proportions of male and female Lama at Tulán-54, slightly 206 more females than male Vicugna at Tulán-54, and more males than females at Tulán-85 for both 207 genera (Table 1). The one individual that showed the highest proportion of ancestry shared with 208 modern llama (CL72) was a male. We note that while our findings provide insights into sex ratios 209 at these sites, the limited sample size, especially at Tulán-85, presents a challenge in making broad 210 inferences. 211 212 Table 1: Genetic sexes of the individuals with known taxonomic status. Taxonomic status was 213 based on the palaeogenomic results unless <0.001x genome-wide coverage then it was based on the 214 osteological measurements - small = Vicugna, big = Lama. 215 Site Tulán-54 Tulán-85 Genus Lama sp. Vicugna sp. Lama sp. Vicugna sp. # Males 13 8 4 2 # Females 14 10 2 1 # Total individuals 27 18 6 3 Proportion Males 0.48 0.44 0.67 0.67 Proportion Females 0.52 0.56 0.33 0.33 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 216 217

218 The Tulán ravine, occupied during the Late Archaic and Early Formative periods, contains 219 both ceremonial and residential sites with archaeological evidence for camelid exploitation, fibre 220 processing, and ritual activity. Investigating palaeogenomic data, we identify the genetic ancestry 221 and sex of camelids from these contexts and interpret these findings in light of the archaeological 222 record to better understand the composition of camelid populations present during this critical 223 period of economic and social transition in the South-Central Andes. Specifically, we investigated 224 49 individuals with mean genome-wide coverages >0.001x from the ceremonial sites of Tulán-52 225 (Late Archaic; n=1) and Tulán-54 (Early Formative; n=40) as well as the residential site of Tulán-226 85 (Early Formative; n=8). We found evidence for both genera, Vicugna and Lama, in our dataset. 227 Further filtering to 26 individuals >0.01x showed that all individuals, except one, likely represented 228 ancestry not found in modern individuals (Fig 2). Moreover, while our ancient data quality 229 undoubtedly presented some limitations, particularly due to the low coverage and poor DNA 230 preservation typical of temperate regions 30, careful evaluation of biases through simulated data 231 demonstrated that even low-coverage palaeogenomes can yield valuable insights when these biases 232 are accounted for. This highlights the potential of low quality data to provide meaningful results, 233 even from challenging archaeological contexts. 234 235 The Vicugna individuals we identified most likely belonged to the subspecies V.v.vicugna, 236 which was not unexpected given that the local subspecies of Vicugna is V.v.vicugna 31. The 237 presence of small individuals with characteristics indicative of vicuñas in the Tulán sites may 238 reflect the value of vicuña fibre, which has been appreciated since pre-Hispanic times 7, and is still 239 greatly appreciated today. This may be the reason that most of the small animals in Tulán 240 correspond to vicuñas. The use of camelid fibre is supported by microscopic camelid fibre analyses 241 32. While the overwhelming majority of fibres were light brown natural colours, a small quantity of 242 pure black fleece suggests colour variation, hinting at early colour selection of pre-domesticated 243 camelids 33. While this could suggest a lost domesticate-like lineage filling a similar role as the 244 alpaca, this could also be the result of dye so should be considered with caution 34. 245 246 Lama individuals were more difficult to place than the Vicugna. All individuals had a 247 tendency towards L.g.cacsilensis but all bar one had considerable levels of L.g.guanicoe ancestry. 248 The current distribution indicates that L.g.cascilensis are generally found further north but may 249 slightly overlap with L.g.guanicoe in Northern Chile 35,36. Therefore, the ancient Tulán individuals 250 may represent a locally extinct L.g.cacsilensis population that hybridised with the local 251 L.g.guanicoe individuals. The one individual more closely related to L.g. cacsilensis, with relatively 252 high levels of llama ancestry, is particularly interesting from a domestication standpoint. When 253 placed in the broader archaeological context of the Tulán ravine, where camelids appear to have 254 been used in a domesticated role, these genetic findings support the idea that this individual may 255 represent either an early domesticated llama or a wild population ancestral to contemporary llamas, 256 with full domestication likely occurring later. Being an early domesticated llama may have also 257 been the case for our one sample from the Late Archaic Tulán-52 individual (C46) which contained 258 mostly L.g.cacsilensis/llama ancestry with relatively little presence of L.g.guanicoe. There is a 259 cultural continuity between the Late Archaic and the Early Formative, so it is to be expected that 260 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint these populations continued to develop the domestication/specialisation process. However, the very 261 low coverage (0.008x) nature of this sample prevented us from delving deeper. Taken together, our 262

suggest individuals that were not the direct ancestors of contemporary domesticated species 263 were the main source of South American camelids used in both domestic and ceremonial activities 264 in the Tulán ravine during the Early Formative period (3,360–2,370 cal. yr BP). 265 266 Osteometric analyses of the excavated samples of Tulán-54 previously revealed the 267 presence of the two wild species: guanaco and vicuña, and the domestic species, llama 26. However, 268 the large degree of overlap between domestic and wild individuals makes clear results difficult to 269 obtain 37. Taking these difficulties together with the fact that there are other lines of evidence (e.g. 270 artefacts and fibres 38,39) supporting the notion that there were domesticated animals in the Tulan 271 ravine since the Late Archaic period 40, suggest the large animals were utilised in a domesticated 272 role despite sharing different ancestry to contemporary llama. As it has been hypothesised that 273 llama domestication was a complex process and may have occurred in several different areas of the 274 Andes 16, the Tulán individuals could represent a domestic lineage that was subsequently replaced 275 or diluted through admixture with other populations. Similar patterns of replacement through 276 introgression have been documented in other domestic mammals, such as pigs, where Near Eastern 277 domestic ancestry was largely replaced by European wild boar ancestry following introduction 41, 278 and horses, where ancient domestic lineages were replaced by later lineages of different origin 42. 279 While the exact mechanisms and timing may differ, these cases illustrate how domesticated 280 lineages can be transformed or supplanted through subsequent gene flow. 281 282 A previous study using mitochondrial genomes of the same individuals suggested that 283 individuals whose size did not match their maternal ancestry (e.g. small size but Lama ancestry or 284 big size but Vicugna ancestry) may have been hybrid individuals, which was interpreted as 285 indicating domesticated individuals 17. Our palaeogenomic analyses find that these potential hybrid 286 individuals do not contain clearly higher levels of shared ancestry between genera. This finding 287 challenges the assumption that mitochondrial and nuclear signatures, as well as osteometric 288 measurements, should align to accurately determine domestication status. As mitochondrial DNA 289 reflects only the maternal lineage, it may not fully capture genetic complexities such as 290 intergenerational hybridisation and introgression. Moreover, osteometric measurements alone, 291 though valuable in identifying size variation, do not account for the possibility of mislabelling or 292 misidentification of juveniles versus adults. As such, we propose that the apparent discrepancy 293 between mitochondrial and nuclear data regarding hybrid status is likely due the inherent 294

of mitochondrial analysis. Consequently, the assumption that size, mitochondrial 295 ancestry, and nuclear signatures should consistently align may not always hold due to the 296 complexity of domestication processes and the potential for interbreeding between different 297 camelid species and subspecies over time. While mitochondrial introgression may have occurred, a 298 single sporadic mitochondrial introgression event could lead to inferences of widespread 299 hybridisation, which, as shown by our nuclear genomic analysis, is not supported. 300 301 Despite early Spanish observers acknowledging the importance of camelid pastoralism, they 302 generally failed to distinguish between domestic species in their records. This contrasts with Lama 303 and Vicugna having long been distinct, biologically divergent genera separated approximately 2–3 304 million years ago 20 and culturally differentiated by pre-Hispanic communities 4. The subsequent 305 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint indiscriminate hybridisation between domesticates resulted in a loss of their exquisite fibre qualities 306 and colourations that made them so prized in prehistoric Andean communities 4. This uncontrolled 307 hybridisation between domesticated species and the population decline post-conquest makes 308 identification of pre-hispanic genetic ancestry of modern individuals difficult. Previous work 309 exclusively analysing genomes from modern individuals has found both llama and alpaca to contain 310 evidence of intergeneric gene flow tracing back to ~500 years ago 20. Our results show that the 311 putative llama in our dataset contained less Vicugna ancestry than modern llama, and similar 312 amounts to wild L.g.cacsilensis, the likely ancestral (sub)species of modern llama. Therefore, we 313 suggest that intermixing between Vicugna (most probably alpaca) and llama was likely not a key 314 component of the early domestication process, and only occurred within the last ~3,000 years. 315 While the absence of evidence for early hybridization in our samples does not exclude its sporadic 316 occurrence, it suggests that hybridisation was not a widespread practice during the studied period. 317 Palaeogenomic data from more domesticated individuals will undoubtedly help shed further light 318 on this issue. 319 320 In addition to the insights gained from the ancestry analysis, we were also able to infer the 321 sex of the individuals using the palaeogenomic data. South American camelids are sexually 322 monomorphic 43 showing sexual dimorphism when specific bones are available (e.g. pelvis and 323 canine teeth) 44. While it is possible to sex South American camelids, this can rarely be done due to 324 the fragmentary nature of most assemblages. Therefore, the palaeogenomic data presented a unique 325 opportunity to make inferences on sex. Phalanxes and Astragali have been used in this study, due to 326 their high frequency, completeness, and in the case of the former, as their fusion stage indicates the 327 presence of adult animals (over 24 months) 45. Thus, our results examine only the adult sex ratio in 328 archaeological assemblages. 329 330 Sex distribution in archaeological assemblages reflects both hunting and herding strategies. 331 Hunting assemblages are generally expected to be dominated by adult remains, with sex ratios 332 balanced or male-biased, whereas herding assemblages often include an adult sex ratio skewed 333 toward females 46. Herding practices can vary, for example, through selective culling of infertile or 334 undesirable males, which may produce more balanced adult sex ratios 47. In the Tulán assemblages, 335 we observe roughly similar proportions of males and females across sites and genera (Table 1), 336 although the small sample from Tulán-85 requires cautious interpretation. These patterns are 337 consistent with either hunting or herding with selective male culling and align with our genomic 338 findings showing that most individuals at these sites were wild individuals and not the direct 339 ancestors of domesticated camelids. 340 341

342 343 Archaeological contexts 344 It is proposed that the location of Archaic and Formative settlements in the Tulán area 345 reflects a combination of favourable cultural, social, and environmental conditions. These 346 conditions supported a shift from a society primarily reliant on hunting and foraging, characterized 347 by seasonal mobility, small group sizes, and exploitation of wild resources, to one increasingly 348 oriented toward pastoralism, marked by the management of domesticated camelids, more 349 permanent settlement structures, and reliance on animal products for subsistence. 350 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint The Late Archaic period (ca. 5290-4090 cal. yBP) is marked by large campsites like Tulán-351 52, the emergence of substantial architecture with clusters of circular structures, and a diversified, 352 innovative lithic industry. This period also saw the appearance of rock art, long-distance 353 interactions, and a subsistence strategy that combined hunting with early camelid domestication. 354 Altogether, these features point to a growing sociocultural complexity 48. 355 The Early Formative period (ca. 3360 - 2370 cal. yBP) is characterised by the emergence of 356 large settlements with a complex and planned ritual architecture along with the intensification and 357 expansion of productive practices, an emphasis on ritual expressions, the appearance of new 358 technologies such as pottery, and three large contemporary Formative settlements (Tu-54, Tu-122 359 and Tu- 85) have been documented for this period. These are located within a 15 km area between 360 the border of the Salar de Atacama and the Tulán ravine, at an altitude ranging from 2300–3200m. 361 These sites have different functions. At Tulán-122, located in the middle sector of the ravine, 362 residential aspects were emphasised. In contrast, Tulán-85, is located on the border of the resource-363 rich salt flat. These settlements shared a common political and religious organisation centred 364 around the Tulán 54 site, the only site that possesses monumental architecture 25. 365 The sunken temple structure identified in Tulán-54 at the centre of a refuse mound 366 contained 27 human infant graves accompanied by offerings, rock art, structured fireplaces, and 367 offering pits, among other elements, suggesting that the structure had a ceremonial function. The 368 main structure identified at the site is a large, sunken semi-oval ceremonial structure surrounded by 369 a perimeter wall made of large, vertically positioned stone blocks. The inner space is divided into 370 six precincts, with a single oval structure in the middle 25,26. The recovered camelid anatomical 371 units display a very similar temporal and spatial pattern across all occupations and among the 372 different precincts of the ceremonial structure, which was filled intentionally with the remains of 373 butchering, processing, consumption and feasting activities. The waste was likely generated outside 374 and then ritually incorporated 49. On the border of the Salar de Atacama, Tulán 85 corresponds to an 375 extensive occupation characterised by dense domestic monticulated deposits. The site also includes 376 marginal structures and a sector with burials of newborns 50. Finally, Tulan-94 site is composed of 377 groups of circular and sub-circular structures, with a total of 19 living enclosures. The evidence 378 suggests the transitional character of the site, since it has material components corresponding to the 379 Late Archaic and Early Formative phases 25. 380 381 Samples 382 All research was conducted in accordance with Chilean regulations. Excavation permits 383 were issued by the Consejo de Monumentos Nacionales de Chile (CMN) under permit Nº 4409 384 (18.11.2013). Export of osteological samples (Camelidae) to the University of Copenhagen for 385 ancient DNA analysis was authorized by CMN ordinances ORD. Nº 03642 (21.10.2016), Nº 03746 386 (27.10.2016), and Nº 03425 (20.01.2017), within the scope of FONDECYT project 1130917. All 387 exported samples were photographed, labeled, and transported via certified courier in accordance 388 with CMN requirements. Although none of the camelid remains analysed here were directly 389 radiocarbon dated, chronological placement is based on radiocarbon dates from associated 390 archaeological contexts (e.g., stratified deposits, fireplaces, and burials) 51. 391 Our study consisted of 75 ancient individuals sourced from four archaeological sites: Tulán-392 52, Tulán-54, Tulán-85, and Tulán-94, and have previously undergone mitochondrial and 393 morphological analysis 17 (Supplementary tables S1 and S2). In brief, samples were collected from 394 both ceremonial and domestic contexts across different sectors of the archaeological sites to capture 395 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint spatial and temporal variation. Bone elements (primarily first anterior phalanges and astragali) were 396 selected from stratified excavation units and structures, with care taken to sample different 397 individuals. Only adult specimens, identified by fused epiphyses or advanced ossification, were 398 included in the morphometric analysis. Measurements were taken using Vernier calipers to the 399 nearest 0.1 mm, following established osteometric protocols, and were recorded prior to aDNA 400 sampling to avoid damage-related bias. Key dimensions were measured on each bone type to assess 401 size variation, and results were compared against published morphometric data for all four modern 402 South American camelid species to evaluate patterns of morphological differentiation. The 403 threshold for distinguishing small and large individuals was determined based on a comparison of 404 key osteometric measurements from the first anterior phalanx (breadth of proximal articulation 405 [BFp] and depth of the proximal epiphysis [Dp]) and the astragalus (breadth distal [Bd] and greatest 406 length medial [GLm]). The distinction between small and large groups was tested statistically by 407 correlating measurements from both the phalanx and astragalus using PAST 52. Published 408 morphometric data for all four modern South American camelid species were used to establish the 409 thresholds for size differentiation 53,54. 410 To build a genomic reference panel, we also included published genomic data from 80 411 modern individuals (DNAzoo.org, 20,35). The modern individuals represented all extant Lama and 412 Vicugna subspecies (Figure 1). 413 414 Bioinformatic processing 415 The raw sequencing data for the ancient Tulán individuals was generated in a previous study 416 where only the mitochondrial DNA was analysed and obtained from the authors 17 while the raw 417 sequencing reads for the modern individuals were downloaded from the European Nucleotide 418 Archive (ENA) (Supplementary tables S1 and S2). For all individuals, we removed Illumina 419 adapter sequences, low-quality reads (mean q <25), short reads (<30bp), and merged overlapping 420 read pairs with Fastp v0.23.2 55. We mapped the processed reads (merged only for the ancient 421 individuals) to the alpaca reference genome (Genbank accession: GCF_000164845.3) using 422 Burrows-wheeler-aligner (BWA) v0.7.15 56 and either utilising the aln algorithm, with the seed 423 disabled (-l 999) (otherwise default parameters) for the ancient individuals or the mem algorithm 424 with default parameters for the modern individuals. We chose the alpaca genome as the reference 425 genome due to its hybrid ancestry 20, making it putatively more suitable when making genomic 426 comparisons between genera. We parsed the alignment files and removed duplicates and reads of 427 mapping quality score <30 using SAMtools v1.6 57. We checked for ancient DNA damage patterns 428 of the ancient individuals using mapdamage v2 58. As the ancient individuals had also undergone 429 targeted enrichment for mitogenomes and some selected genes, we aligned the probes to the 430

genome using BLAST v 2.15.0 59 using default parameters and we removed any regions 431 of the alpaca genome aligning to the probes from the mapped bam files using bedtools v2.29.1 60. 432 433 Genetic sex determination 434 We identified scaffolds putatively originating from the sex chromosomes in the alpaca 435 assembly by aligning the assembly to the alpaca X (included in the reference genome) and Human 436 Y (Genbank accession: NC_000024.10) chromosomes using satsuma synteny v2.1 61 with default 437 parameters. We followed the SeXY pipeline 29 to identify the sex of each of the ancient individuals. 438 We set a minimum threshold of 5,000 mapped reads to ensure the reliability of the results. If an 439 individual had a mean X:A ratio of <0.7 it was designated as a male. If an individual had a mean 440 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint X:A ratio of >0.8 it was designated a female. Individuals with mean ratios between 0.7 and 0.8 441 were deemed undetermined. 442 443 Ancient DNA Simulation 444 To assess the impact of ancient DNA damage on our analyses we simulated the damage 445 patterns and read lengths of individual C21 onto representatives of each modern subspecies (seven 446 in total) using gargammel v1.1.4 62. We chose C21 due to it having the highest coverage of the 447 ancient samples. To do this, we generated consensus fasta files from the mapped bam files of the 448 modern individuals using ANGSDv0.92163 and the following parameters -dofasta 2, minimum read 449 depth of 5 (-setmindepthind 5), and minimum mapping and base qualities of 30 (-minq 30 -450 minmapq 30). The modern individuals included L.glama (Llama2 and Llama3), V.pacos (Alpaca3), 451 L.g.cacsilensis (Cacsilensis3), L.g.guanicoe (Guanicoe3), V.v.mensalis (Mensalis2), and 452 V.v.vicugna (Vicugna1). We mapped the simulated ancient DNA reads to the alpaca reference 453 genome following the same approach as for the ancient Tulán individuals. We additionally 454 downsampled them to ~0.01x using SAMtools. 455 456 Principal Component Analysis 457 We evaluated the relationships between the modern and ancient camelids included in this 458 study by performing PCAs using a pseudohaploid base call with ANGSD. We initially did this for 459 two different sample sets, either for all modern and ancient individuals >0.01x or all modern and 460 ancient individuals >0.001x. We computed consensus pseudohaploid base calls specifying the 461 parameters: minimum mapping and base qualities of 20 and 30 respectively (-minmapQ 20 -minQ 462 30), calculate genotype likelihoods using the GATK algorithm (-GL 2), calculate major and minor 463 alleles based on genotype likelihoods (-doMajorMinor 1), remove transitions (-rmtrans 1), only 464 consider autosomal chromosomes > 10Mb (-rf), skip triallelic sites (-skiptriallelic 1), only consider 465 reads mapping to one region uniquely (-uniqueonly 1), make a consensus identity by state base calls 466 (-doIBS 2), output a covariance matrix (-doCov 1) and only consider variable positions where the 467 minor allele occurs in at least two individuals (-minminor 2). We also repeated this analysis but 468 with only individuals belonging to the Lama genus. We set the minimum individual threshold to the 469 number of modern individuals +1, i.e. 81 for the complete dataset and 57 for the Lama only. 470 471 Phylogenetic Analysis 472 To build a neighbour joining tree, we first constructed a distance matrix, which we output 473 from ANGSD when performing PCA (setting the parameter -makematrix 1). We converted the 474 distance matrix into a neighbour joining phylogenetic tree using fastME v2.1.6.1 64 and default 475 parameters. 476 477 Ancestry proportions 478 We calculated the ancestry proportions of eight ancient individuals using admixfrog v0.7.2 479 27. We included four big size camelids (C21, C38, C72, C75) and four small size camelids C20, 480 C39, C50, C69. We selected these based on their higher coverage or unexpected placement within 481 the PCA. Individuals C69 and C38 did not correspond to what was expected based on their size and 482 C75 clustered differently compared to the other big sized ancient individuals. We also ran this 483 analysis on the lone individual from Tulán-52 (C46) as it comes from an earlier time period which 484 may be closer to the onset of domestication. Initially we only considered the wild species as the 485 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint ancestry states, each represented by a single individual; L.g.cacsilensis - Cacsilensis4, L.g.guanicoe 486 - Guanicoe1, V.v.mensalis - Mensalis4, and V.v.vicugna - Vicugna2. We repeated the analysis for 487 the ancient Lama individuals using the same two Lama sp. subspecies as above plus Llama9 as 488 ancestry states. As input for the reference panel we created a multi-individual variant call file with 489 using BCFtools v1.15 65, specifying autosomal scaffolds >10Mb (-f) and minimum mapping and 490 base qualities of 20 (-q 20 -Q 20). We ran admixfrog specifying a minimum read length of 30bp, a 491 bin size of 50kb, and otherwise default parameters. We only considered a bin in interpretations if it 492 contained more than five SNPs. We pooled regions showing mixed ancestry between Lama sp. and 493 Vicugna sp. as ancestral camelid in the first run and pooled regions showing any mixed ancestry in 494 the Lama sp. only run as ancestral. Furthermore, to assess the accuracy and suitability of this 495

on our dataset, we also performed the same analysis (all wild (sub)species reference panel) 496 on the modern individuals with simulated ancient DNA damage. 497 498 D-statistics 499 D-statistics calculates shared derived alleles between non-sister branches of a given 500 topology [[[H1, H2], H3], Outgroup]. An allele is considered as derived if it is different to the 501 outgroup allele. Although it is most commonly implemented to investigate differential levels of 502 gene flow between non-sister taxon, it can also be used to investigate shared ancestral 503 polymorphisms and therefore topology/population structure 66,67. 504 We subsampled the dataset to include four modern individuals per subspecies and only the 505 eight ancient individuals also used in the ancestry proportions analysis above. We calculated the D-506 statistics using a random base call approach in ANGSD (-doabbababa 1), specifying only autosomal 507 scaffolds >10Mb and the following parameters: -minmapQ 30 -minQ 30 -blocksize 1000000 -508 rmtrans 1 -uniqueonly 1. To infer the ancestral sequence (-anc) we used a wild bactrian camel 509 (Camelus ferus; NCBI biosample SRR1947250) which we also mapped to the alpaca genome 510 following the approach mentioned above. We summarised the results using a block jackknife 511 approach with the Rscript available in the ANGSD toolsuite (ANGSD_jackknife.R). 512 We filtered the output into three different sets. First to investigate which genus each ancient 513 individual was most closely related to, we filtered for H1=Lama sp. H2=Vicugna sp. H3=Ancient 514 individuals. Based on these results we then looked into which subspecies the ancient individuals 515 were most closely related to with either H1=Lama sp. H2=Lama sp. H3=Ancient Lama sp. 516 individuals or H1=Vicugna sp. H2=Vicugna sp. H3=Ancient Vicugna sp. Individuals. To 517 investigate gene flow between the ancient Lama sp and Vicugna sp., we extracted results 518 corresponding to the topologies H1=Lama sp., H2=C21, H3=Vicugna sp. H1=Lama sp., H2=C72, 519 H3=Vicugna sp. 520 To investigate whether the ancient DNA damage could cause biases in the topology test 521 results, we compared results when placing the individuals with simulated ancient DNA damage in 522 H3 with the results generated using the high quality version of the same individual. To evaluate the 523 ability of our dataset to infer gene flow between the putative ancient llama and Vicugna species, we 524 also performed D-statistics with the two 0.01x simulated aDNA llama (Llama 2 and 3) individuals 525 in the H2 position, Vicugna sp in the H3 position and the remaining Lama sp. individuals in the H1 526 position. By comparing these results to those when the high quality versions of llama 2 and 3 were 527 in the same position, we were able to calculate a correction factor to account for potential biases. 528 We found a mean deviation between the two results of 0.0322 (Supplementary figure S6) which we 529 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint subtracted from the D-score and using the standard error output with ANGSD, we recalculated the 530 Z scores. 531 532 Figures: 533 534 535 536 Figure 1: Map indicating the location of the Tulán archaeological sites and the sampling 537 location of the previously published modern South American camelid individuals used in this 538 study. 539 540 541 542 543 544 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 545 Figure 2: Relationships between modern and ancient individuals. A) PCA performed on all 546 modern individuals and ancient individuals >0.01x using 5,163,744 transversion sites. B) PCA 547 performed on modern Lama individuals and ancient Lama individuals > 0.01x using 3,668,867 548 transversion sites. C) Ancestry states of eight ancient individuals. Ancestry states are derived from 549 a single individual; L.g.cacsilensis - Cacsilensis4, L.g.guanicoe - Guanicoe1, V.v.mensalis - 550 Mensalis4, and V.v.vicugna. Ancestral is defined as having combined ancestry between both Lama 551 and Vicugna (C) or mixed ancestry between Lama (sub)species (D). 552 553 554 555 Figure 3: Test for closer evolutionary relationships between our ancient samples and the 556 modern (sub)species using D-statistics and various input topologies. A) Comparisons between 557 Lama and Vicugna B) Comparisons within Lama C) Comparisons within Vicugna. Non-significant 558 values |Z|<3 not shown. 559 560 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 561 562 Figure 4: Test for gene flow between ancient individuals and Vicugna sp. using D-statistics. A 563 negative D-score shows relatively more gene flow between H1 and H3 compared to H2 and H3, 564 whereas positive shows relatively more gene flow between H2 and H3. Non-significant values 565 |Z|<3 are shown with open circles. 566 567

568 We would like to thank Anders Hansen for support in generating the sequencing data and Alba 569 Rey-Iglesia for help curating the data. We would also like to thank Lautaro Núñez who through 570 numerous projects (FONDECYT 1020316, 1070040 and 1130917) recovered the analysed camelid 571 bone remains. 572 573 Funding 574 MVW was supported by a Novo Nordisk Emerging Investigator grant #NNF24SA0093839 575 576 Data availability 577 Raw sequencing data is available under NCBI BioProject ID: PRJNA1215279. Coordinates of the 578 capture probes in the alpaca reference genome can be downloaded from 579 https://erda.ku.dk/archives/7d1d408b1711f5ab7e34a845bf58f9a8/published-archive.html 580 581

582 1. Larson, G. & Fuller, D. Q. The Evolution of Animal Domestication. Annu. Rev. Ecol. Evol. Syst. 45583 136 (2014). 584 2. Zeder, M. A. The Domestication of Animals. J. Anthropol. Res. 68, 161–190 (2012). 585 3. Frantz, L. A. F., Bradley, D. G., Larson, G. & Orlando, L. Animal domestication in the era of ancien586 genomics. Nat. Rev. Genet. 21, 449–460 (2020). 587 4. Wheeler, J. C., Russel, A. J. F. & Redden, H. Llamas and Alpacas: Pre-conquest breeds and post- con588 hybrids. J. Archaeol. Sci. 22, 833–840 (1995). 589 45, 115– ient conquest .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 5. Hesse, B. Archaeological evidence for camelid exploitation in the Chilean Andes. Säugetierkundliche 590 Mitteilungen 30, 201–211 (1982). 591 6. Broughton, K. L. South American camelids in Central Andean religious practices. (2010). 592 7. Vilá, B. & Arzamendia, Y. South American Camelids: their values and contributions to people. 593 Sustainability Sci. 1–18 (2020). 594 8. Goñalons, G. L. M. & Yacobaccio, H. D. The domestication of South American camelids. Documenting 595 Domestication--New Genetic and Archaeological Paradigms 228–244 (2006). 596 9. Cartajena Fasting, I. Explorando la variabilidad morfométrica del conjunto de camélidos pequeños durante 597 el Arcaico Tardío y el Formativo Temprano en Quebrada Tulán, norte de Chile. Revista del Museo de 598 Antropología 2, 199–212 (2009). 599 10. Moore, K. M. 9. Grace Under Pressure: Responses to Changing Environments by Herders and Fishers in 600 the Formative Lake Titicaca Basin, Bolivia. in Sustainable Lifeways 244–272 (University of Pennsylvania 601 Press, 2012). 602 11. Moore, K. M. Early domesticated camelids in the Andes. The archaeology of Andean pastoralism (2016). 603 12. Noe, S. J., Haas, R. & Aldenderfer, M. Hunting to herding on the Andean Altiplano: Zooarchaeological 604 insights into Archaic Period subsistence in the Lake Titicaca Basin, Peru (9.0–3.5 ka). J. Anthropol. 605 Archaeol. 77, 101658 (2025). 606 13. Le Neün, M. et al. Holocene occupation of the Andean highlands: A new radiocarbon chronology for the 607 Telarmachay rockshelter (Central Andes, Peru). Quat. Sci. Rev. 312, 108146 (2023). 608 14. Marsh, E. J. & Rademaker, K. Human–environment interactions in the Lake Junín basin: Fire, megafauna, 609 deforestation, and domestication, from the peopling of the Andes to the Inca Empire. Quat. Sci. Rev. 351, 610 109159 (2025). 611 15. Mengoni Goñalons, G. L. Camelids in ancient Andean societies: A review of the zooarchaeological 612 evidence. Quat. Int. 185, 59–68 (2008). 613 16. Yacobaccio, H. D. The domestication of South American camelids: a review. Anim Front 11, 43–51 614 (2021). 615 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 17. Diaz-Maroto, P. et al. Ancient DNA reveals the lost domestication history of South American camelids in 616 Northern Chile and across the Andes. eLife 10, e63390 (2021). 617 18. Le Neün, M. et al. Can first phalanx multivariate morphometrics help document past taxonomic diversity 618 in South American camelids? J. Archaeol. Sci. Rep. 47, 103708 (2023). 619 19. Zeder, M. A., Emshwiller, E., Smith, B. D. & Bradley, D. G. Documenting domestication: the intersection 620 of genetics and archaeology. Trends Genet. 22, 139–155 (2006). 621 20. Fan, R. et al. Genomic analysis of the domestication and post-Spanish conquest evolution of the llama and 622 alpaca. Genome Biol. 21, 159 (2020). 623 21. Abbona, C. C. et al. Were domestic camelids present on the prehispanic South American agricultural 624 frontier? An ancient DNA study. PLoS One 15, e0240474 (2020). 625 22. Westbury, M. et al. First complete mitochondrial genome data from ancient South American camelids - 626 The mystery of the chilihueques from Isla Mocha (Chile). Sci. Rep. 6, 38708 (2016). 627 23. Díaz-Lameiro, A. M. et al. Ancient DNA confirms crossbreeding of domestic South American camelids in 628 two pre-conquest archaeological sites. J. Archaeol. Sci. 141, 105593 (2022). 629 24. Westbury, M. V. et al. Hyena paleogenomes reveal a complex evolutionary history of cross-continental 630 gene flow between spotted and cave hyena. Science Advances 6, eaay0456 (2020). 631 25. Nuñez A, L., Cartajena F, I., Carrasco G, C., de Souza H, P. & Grosjean, M. Emergencia de comunidades 632 pastoralistas formativas en el sureste de la Puna de Atacama. Estud. atacameños 93–117 (2006). 633 26. Núñez, L. et al. The temple of Tulán-54: Early Formative ceremonial architecture in the Atacama Desert. 634 Antiquity 91, 901–915 (2017). 635 27. Peter, B. M. 100,000 years of gene flow between Neandertals and Denisovans in the Altai mountains. 636 bioRxiv (2020). 637 28. Yuan, J. et al. Camelus knoblochi genome reveals the complex evolutionary history of Old World camels. 638 Curr. Biol. 34, 2502–2508.e5 (2024). 639 29. Cabrera, A. A. et al. How low can you go? Introducing SeXY: sex identification from low-quantity 640 sequencing data despite lacking assembled sex chromosomes. Ecol. Evol. 12, e9185 (2022). 641 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 30. Hofreiter, M. et al. The future of ancient DNA: Technical advances and conceptual shifts. Bioessays 37, 642 284–293 (2015). 643 31. González, B. A. et al. Phylogeography and Population Genetics of Vicugna vicugna: Evolution in the Arid 644 Andean High Plateau. Front. Genet. 10, (2019). 645 32. Benavente, M. A. Análisis lanimétrico de fanéreos de los sitios Tulán-52 y 54. Ms Informe Final Proyecto 646 FONDECYT 1020316, (2005). 647 33. Dransart, P. Llamas, herders and the exploitation of raw materials in the Atacama Desert. World Archaeol. 648 22, 304–319 (1991). 649 34. Dransart, P. Camelid Fleece and the Hidden Histories of Colors in the Andes. (2020). 650 35. León, F. et al. History of diversification and adaptation from north to south revealed by genomic data: 651 Guanacos from the desert to sub-Antarctica. Genome Biol. Evol. 16, evae085 (2024). 652 36. González, B. A., Samaniego, H., Marín, J. C. & Estades, C. F. Unveiling current Guanaco distribution in 653 chile based upon niche structure of phylogeographic lineages: Andean puna to subpolar forests. PLoS One 654 8, e78894 (2013). 655 37. López M., P., Cartajena, I., Loyola, R., Núñez, L. & Carrasco, C. The use of hunting and herding spaces: 656 Stable isotope analysis of Late Archaic and Early Formative camelids in the Tulan transect (Puna de 657 Atacama, Chile): Stable isotopes analyses on camelids in tulan transect (Chile). Int. J. Osteoarchaeol. 27, 658 1059–1069 (2017). 659 38. Arias, G., Benavente, A. & Gecele, P. Identificación y variabilidad del uso del animal a través de textiles 660 arqueológicos: Contraste con patrones fanéreos actuales. in Actas del xii congreso Nacional de 661 Arqueología chilena 151–162 (1993). 662 39. Cartajena, I., Benavente, M. A., Núñez, L. & Thomas, C. La utilización de los camélidos durante el 663 Formativo Temprano: una comparación entre la cuenca del Loa medio y el Salar de Atacama. 664 Zooarqueología y Tafonomía en el Confín del Mundo 45–58 (2009). 665 40. Loyola, R. et al. Hunting and feasting in the pre-Columbian Andes: Exploring the nature and scale of early 666 ceremonial aggregations in Tulan Ravine (5300 to 2400 yr cal. BP) through the circulation of obsidian 667 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint artefacts. J. Anthropol. Archaeol. 64, 101360 (2021). 668 41. Frantz, L. A. F. et al. Ancient pigs reveal a near-complete genomic turnover following their introduction to 669 Europe. Proc. Natl. Acad. Sci. U. S. A. 116, 17231–17238 (2019). 670 42. Librado, P. et al. Widespread horse-based mobility arose around 2200 BCE in Eurasia. Nature 631, 819–671 825 (2024). 672 43. Miranda-de la Lama, G. C. & Villarroel, M. Behavioural biology of South American domestic camelids: 673 An overview from a welfare perspective. Small Rumin. Res. 220, 106918 (2023). 674 44. Kaufmann, C. A. Estructura de Edad Y Sexo En Guanaco: Estudios Actualísticos Y Arqueológicos En 675 Pampa Y Patagonia. (Sociedad Argentina de Antropología, 2009). 676 45. Kent, J. D. The Domestication and exploitation of the South American camelids: methods of analysis and 677 their application to circum-lacustrine archaeological sites in Bolivia and Peru. (1982). 678 46. Yacobaccio, H. Andean camelid herding in the South Andes: ethnoarchaeological models for 679 archaeozoological research. Yacobaccio, Hugo Daniel; Andean Camelid Herding in the South Andes: 680 Ethnoarchaeological Models for Archaeozoological Reseach; Muséum National d´Histoire Naturelle; 681 Anthropozoologica 42, 143–154 (2007). 682 47. Wheeler, J. C. Pre-conquest alpaca and llama breeding. Camelid Quarterly (2005). 683 48. De, S. P., Cartajena, I., Núñez, L. & Carrasco, C. Cazadores-recolectores del Arcaico Tardío y desarrollo 684 de complejidad social en la puna de Atacama: las evidencias del sitio Tulan-52 (norte árido de Chile). 685 Revista Werkén 13, 91–118 (2010). 686 49. Cartajena, I., Mendoza, P. L., Santander, B. & Núñez, L. Discarding practices among early herders on the 687 western watershed of the Puna de Atacama, Chile (ca. 3000-2300 cal BP). in Animals: Cultural Identifiers 688 in Ancient Societies? (eds. Peters, J., McGlynn, G. & Goebel, V.) 65–84 (Documenta Archaeobiologiae 15 689 VML Verlag Marie Leidorf, 2019). 690 50. Lopez, P., Cartajena, I. & Nunez, L. Analyses of stable isotopes in camelids collagen bones from Tulan 691 Ravine, Atacama Puna, early formative period (CA 3, 1000-2,400 BP). Revista de antropologia chilena 45, 692 237–247 (2013). 693 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 51. Núñez, L. et al. Presencia de un centro ceremonial formativo en la circumpuna de Atacama. Chungará 694 (arica) 49, 3–33 (2017). 695 52. Hammer, Ø., Harper, D. A. T. & Ryan, P. D. PAST: Paleontological statistics software package for 696 education and data analysis. Palaeontol. Electronica 4, 9 (2001). 697 53. Cartajena, I. Los conjuntos arqueofaunísticos del Arcaico Temprano en la Puna de Atacama, Norte de 698 Chile. PhD diss. , Freie Universität Berlin, Berlin (2002). 699 54. L’Heureux, G. L. Morfometría de camélidos sudamericanos modernos. La variabilidad morfológica y la 700 diversidad taxonómica. Zooarqueología a principios del siglo XXI. Aportes teóricos, metodológicos y 701 casos de estudio 39–50 (2010). 702 55. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 703 34, i884–i890 (2018). 704 56. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. 705 Bioinformatics 25, 1754–1760 (2009). 706 57. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). 707 58. Jónsson, H., Ginolhac, A., Schubert, M., Johnson, P. L. F. & Orlando, L. mapDamage2.0: fast approximate 708 Bayesian estimates of ancient DNA damage parameters. Bioinformatics 29, 1682–1684 (2013). 709 59. Mount, D. W. Using the Basic Local Alignment Search Tool (BLAST). CSH Protoc. 2007, db.top17 710 (2007). 711 60. Quinlan, A. R. BEDTools: The Swiss-Army Tool for Genome Feature Analysis. Curr. Protoc. 712 Bioinformatics 47, 11.12.1–34 (2014). 713 61. Grabherr, M. G. et al. Genome-wide synteny through highly sensitive sequence alignment: Satsuma. 714 Bioinformatics 26, 1145–1151 (2010). 715 62. Renaud, G., Hanghøj, K., Willerslev, E. & Orlando, L. gargammel: a sequence simulator for ancient DNA. 716 Bioinformatics 33, 577–579 (2017). 717 63. Korneliussen, T. S., Albrechtsen, A. & Nielsen, R. ANGSD: Analysis of Next Generation Sequencing 718 Data. BMC Bioinformatics 15, 356 (2014). 719 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint 64. Lefort, V., Desper, R. & Gascuel, O. FastME 2.0: A Comprehensive, Accurate, and Fast Distance-Based 720 Phylogeny Inference Program. Mol. Biol. Evol. 32, 2798–2800 (2015). 721 65. Narasimhan, V. et al. BCFtools/RoH: a hidden Markov model approach for detecting autozygosity from 722 next-generation sequencing data. Bioinformatics 32, 1749–1751 (2016). 723 66. Westbury, M. V. et al. Extended and Continuous Decline in Effective Population Size Results in Low 724 Genomic Diversity in the World’s Rarest Hyena Species, the Brown Hyena. Mol. Biol. Evol. 35, 1225–725 1237 (2018). 726 67. Segawa, T., Rey-Iglesia, A., Lorenzen, E. D. & Westbury, M. V. The origins and diversification of 727 Holarctic brown bear populations inferred from genomes of past and present populations. Proc. Biol. Sci. 728 291, (2024). 729 730 .CC-BY-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 6, 2025. ; https://doi.org/10.1101/2024.10.09.617380doi: bioRxiv preprint