Ancient genomes reveal long-range influence of the pre-Columbian culture and site of Tiwanaku.

Tiwanaku civilization flourished in the Lake Titicaca basin between 500 and 1000 CE and at its apogee influenced wide areas across the southern Andes. Despite a considerable amount of archaeological data, little is known about the Tiwanaku population. We analyzed 17 low-coverage genomes from individuals dated between 300 and 1500 CE and demonstrated genetic continuity in the Lake Titicaca basin throughout this period, which indicates that the substantial cultural and political changes in the region were not accompanied by large-scale population movements. Conversely, the ritual center of Tiwanaku revealed high diversity, including individuals with primarily local genetic ancestry and those with foreign admixture or provenance from as far as the Amazon. Nonetheless, most human offerings associated with the Akapana platform exhibited pure Titicaca basin ancestry and dated to ca. 950 CE—the onset of Tiwanaku’s decline as a sociopolitical center. Our results strengthen the view of Tiwanaku as a complex and far-reaching polity.

INTRODUCTION

The motivations and intentions shaping the transition from small-scale societies to settled life—i.e., villages and cities—are one of the primary questions in archaeology. Located 3850 m above sea level near the shores of Lake Titicaca, Tiwanaku represents the unlikeliest case of the spontaneous and primary rise of social complexity on par with a handful of other locations on Earth (Fig. 1A). For almost half a millennium (500 to 1000 CE), Tiwanaku was one of the most influential centers in the southern Andes (Fig. 1B). More than 100 years of archaeological research has revealed how cultural and demographic changes in the Lake Titicaca basin preceded Tiwanaku’s emergence as the primary ritual center around 500 CE. The southeastern basin enjoys a milder climate than the other parts of the Altiplano and was densely inhabited from at least the Early Middle Formative Period (1500 to 100 BCE), during which small, politically independent villages existed. Trade intensification and the appearance of a common religious tradition led to the formation of multicommunity polities and culminated with the initial formation of Tiwanaku around 300 CE.

Fig. 1.

Provenience of the individuals and groups analyzed in this study.

(A) Geographic locations of the groups in South America and (B) map focusing on the Central Andes; present-day political borders provided for reference. (C) Extent of Tiwanaku influence, showing main sites and locations from where TIW, LUK, and ORU individuals were sampled for this study (red stars) [redrawn from (); map to scale around the lake, not to scale in outlying regions]. (D) Location of monuments at the site of Tiwanaku [contours redrawn from ()].

The emergence of Tiwanaku was apparently accompanied by socioeconomic, ideological, and demographic transformations (–). For example, Khonkho Wankane, a major center in the Desaguadero Valley, was largely abandoned by the end of the Late Formative Period (100 to 500 CE), whereas Lukurmata remained an important settlement in the adjacent Katari Valley (). Tiwanaku is remarkable in that it was the first archaic state in the Titicaca basin to expand beyond its core area and, in this process, subordinated or otherwise influenced a large territory in the southern Titicaca region [Fig. 1C; ()].

Tiwanaku was the first large urban center in the Titicaca area, with an estimated population of as many as 20,000 inhabitants during peak ceremonial events (, ). The monuments of Tiwanaku are spread over several kilometers of a flat plain that lies to the southeast of the colonial and modern town of Tiwanaku (Fig. 1D) and include a half dozen primary ceremonial buildings accompanied by numerous minor constructions. As in other Andean cultures, the ritual practices of Tiwanaku forged a degree of political unity among diverse groups (). The ritual core of the city of Tiwanaku was designed according to beliefs concerning the unity between cosmological, mythical, and physical spaces (). The earliest monument at the site, the semi-subterranean temple, was complemented in the third to fifth century CE by the adjacent Kalasasaya platform and Kheri Khala complex (). Around the seventh century CE, a substantial transformation at the site took place with the construction of the Akapana platform (). A few centuries later, earlier buildings were razed, provisioning the construction of the Putuni complex (). The construction of another stepped platform, the Pumapunku temple, started to the southeast in the mid-seventh century CE, and the structure continued to be modified for centuries (Fig. 1D). Between and around these monuments’ plaza areas were walled compounds that included homes, kitchens, and large patios ().

The site is dominated by Akapana, a monumental 17-m-tall platform constructed in the shape of a half Andean cross. This construction consisted of seven stepped terraces composed of stone revetments and earthen fill (, ). Archaeological evidence suggests that it was a place of ritual ceremonies, including commensal feasting and offerings. At the platform’s base, deposits of human remains were found alongside mixed camelid bones and ceramic sherds; these human remains were interpreted by the excavators as ritual offerings (–). The analysis of these remains found evidence of dismemberment at the time of, or soon after, death (). According to the authors of this analysis, these sacrifices were used by elites to attest their power by manifesting their connection with the deities. Another interpretation is an emerging elite class dismembering sacred ancestors of recently conquered and now subjugated populations (, ).

The Tiwanaku culture was preliterate, which means that archaeological research is the sole source of primary information. The past century of research has provided layout, size, and distribution of the buildings present at the site. However, fundamental questions regarding the site’s population demography remain. Scholarly opinions range from a densely occupied city to a near-empty ceremonial center that cyclically pulsated with life from seasonal pilgrims (). Regardless of the habitation pattern of the city of Tiwanaku, the archaeological and bioarcheological investigations suggest the presence of diverse groups of individuals at the site (, , ). These findings were recently supported by the first genomic data from three individuals excavated at the site, one of whom exhibited nonlocal ancestry ().

Explanations of the mechanism behind the apparent population diversity at the site range from bellicose to benign—from forced migration to voluntary pilgrimage and trade. Tiwanaku had the capacity to incorporate other polities and mobilize labor in its heartland territory and selected enclaves (Fig. 1C). Its political economy was based not only on local production, including raised field agriculture, camelid herding, and lake exploitation, but also on extensive trade and colonial relationships throughout the South-Central Andes to secure access to lowland resources, especially maize (, ).

In this study, we use a combination of archaeological and paleogenomic approaches to gain insight into the population structure and dynamics in the Central Andes, with a focus on the Tiwanaku ceremonial center. We examine genome-wide information for 13 ancient individuals from the Bolivian Lake Titicaca region associated with the Tiwanaku culture (500 to 1000 CE) and four from the region of Coropuna volcano in southern Peru associated with the Wari (500 to 1000 CE) and Inca (ca. 1400 to 1540 CE, in this region) cultures. We analyze the genetic makeup of these groups, compare their affinities to other ancient and modern populations, and review the absolute ages of the examined remains. The results of our study reveal the complex demographic history of the population at Tiwanaku.

Ethics statement

All samples for the individuals included in this study were obtained under permission from the local authorities. Individuals from Bolivia were sampled for analyses under permission from La Unidad de Arqueologia y Museos (UDAM) Ministerio de Culturas y Turismo de Bolivia nos. 052/2016 and 086/2016. The samples from Peruvian individuals were collected under permission granted by the Peruvian Ministerio de Cultura (formerly Instituto Nacional de Cultura). The results obtained in this study will be provided to authorities in the Centro de Investigaciones Arqueológicas, Antropológicas y Administración de Tiwanaku (CIAAAT) and the Museo Arqueológico de la Universidad Católica de Santa María in Arequipa (UCSM). The results of this study will be made available and disseminated among the local communities through these institutions.

RESULTS

Sequencing statistics and authenticity

We screened 93 specimens sampled from pre-Columbian sites near Lake Titicaca in Bolivia and southern Peru (dataset S1, A and B) for DNA preservation. We used various approaches to increase the very low content of endogenous DNA, including whole-genome capture (), predigestion (), and drilling only the outer layer of the root (cementum) in the case of teeth (). We generated low-coverage genome sequences for 18 individuals with a depth of coverage ranging from 0.15× to 2.56×. The sequencing data showed deamination patterns at 5′ and 3′ ends and mean fragment lengths that were characteristic of ancient DNA (dataset S1B). The authenticity of the data was corroborated by a lack of detectable contamination in the nuclear DNA, low estimates of contamination for the mitochondrial DNA, and a lack of heterozygosity on the X chromosome in male specimens (all estimated to be below 5%) (dataset S1A). One individual (TW098) did not meet the quality control threshold for ancient DNA authenticity, with nuclear contamination estimated at approximately 9%, and was removed from further analyses. We estimated slightly higher contamination for the X chromosome (5.5%) for individual TW059, but considering that it was based on a limited number of single nucleotide polymorphisms (SNPs) and all other contamination estimates met the threshold, we decided to retain this individual for downstream analyses.

Archeological context and radiocarbon dating

Dating of the ancient individuals based on archaeological context and stratigraphy was confirmed by direct radiocarbon dating of each individual’s skeletal remains (Table 1, dataset S1C, and figs. S1 and S2). Our dataset contains 13 individuals from the Lake Titicaca basin, spanning the period between ca. 300 and 1500 CE. The residential group from outside the Tiwanaku ritual core site consisted of five individuals. Of these, four originated from the site of Lukurmata in the nearby Katari Valley (TW013, TW020, TW027, and TW028; hereafter LUK). Individual TW033 (hereafter ORU) was excavated in Totocachi, Oruro, separated from Katari Valley by 200 km.

Table 1.

Radiocarbon dates and main genetic indices of the individuals analyzed in this study. n.a., not applicable.

ID	Group	Site	cal CE*	Genome coverage	Genetic sex	mtDNA haplogroup	Y haplogroup
TW004	TIW	Akapana	870–990	0.32	M	C1b	Q1b1a1a1
TW008		Akapana	880–990	0.15	M	B2	Q1b1a1a1w
TW056		Monolito Descabezado	890–990	1.20	F	C1d1	n.a.
TW059		Pumapunku	680–740 (21.3%),760–880 (74.1%)	0.19	M	B2	Low cov
TW060		Akapana	900–930 (10.9%),960–1030(84.5%)	0.19	F	D1	n.a.
TW061		Akapana	1020–1150	0.14	F	B2b	n.a.
TW063		Putuni	670–780 (93.5%),820–840 (1.9%)	0.58	F	C1c	n.a.
TW097		Akapana	880–970	0.21	M	C1c	Q1b1a
TW013	LUK	Lukurmata	200–370	2.56	F	C1b	n.a.
TW020		Lukurmata	1020–1150	0.17	F	B2	n.a.
TW027		Lukurmata	980–1050(95.0%),1100–1110(0.4%)	0.34	M	B2	Q1b1
TW028		Lukurmata	1430–1510(82.7%),1580–1620(12.7%)	1.23	M	B2	Q1b1a1a1i
TW033	ORU	Totocachi	1390–1440	0.14	F	B2	n.a.
CO001	COR	Culcunche	770–820 (5.1%),840–980 (90.3%)	0.43	F	B2	n.a.
CO066		Maucallacta	1450–1510(61.2%),1570–1630(34.2%)	0.4	M	B2b	Q1b1a1a1i
CO154		Antaura	1490–1630	0.9	F	A2d1	n.a.
CO193		Cotahuasi	1280–1320(48.7%),1350–1390(46.7%)	0.45	F	B2	n.a.

*Dates (95.4% probability) calibrated in OxCal v4.3.2 () using SHcal13 calibration curve ().

Radiocarbon dating revealed that the LUK group contained individuals representing three periods. Individual TW013 was dated to ca. 300 CE and predated the incorporation of the Katari Valley into the later Tiwanaku polity. Two other individuals (TW020 and TW027) were dated to 1100 and 1010 CE, respectively, i.e., to the period in which Tiwanaku lost its status as a sociopolitical center. Last, individual TW028 was dated to 1470 CE—to the period of the occupation of the area by the Inca Empire. The ORU individual is dated to after the abandonment of the Tiwanaku site and perhaps already represents a later culture.

The other eight individuals originated from five different locations within the ritual core of the Tiwanaku site (TIW): TW060, a part of a complex collection of human and animal bones interpreted as an offering, and TW061, an infant, were placed in alluvia accumulating along the base of the Akapana platform. TW059 is a full individual placed face down and covered in construction fill of the Pumapunku and thus a strong candidate for a sacrifice marking a substantial modification to this platform. TW063 is a solitary skull found in the disturbed and heavily quarried south exterior wall of the Putuni platform. TW097 individual was placed in a fill for a pebbled surface constructed between the semi-subterranean temple and the Akapana, with the position of the limbs consistent with a person wrapped in a burial bundle. TW056 is a partial individual mixed with a midden of ceramics, bones (human and animal), and ash adjacent to a monolith (Monolito Descabezado) along the northeastern corner of the Kalasasaya platform (Fig. 1D). A position of a more complete individual from this location, consistent with someone that was bound, suggests a sacrificial character of this deposit. However, DNA preservation was insufficient for pursuing downstream analyses for this individual (dataset S1A). Two individuals do not have a precise provenience (TW004 and TW008), but they are suspected of having come from excavations at the Akapana platform. In the “TIW” group used in downstream analyses, we also include three previously published individuals from Akapana: I0976, I0977, and I0978 ().

Two individuals (TW059 and TW063) lived in a period of strong Tiwanaku influence throughout the basin and the southern Andes (700 to 800 CE). Six are dated to approximately 950 CE, during what is considered to be the onset of decline, and one dates to around 1100 CE, the period after active maintenance of the Tiwanaku site. In addition to the individuals directly associated with the Tiwanaku culture, genomic data were generated for four individuals from the Late Intermediate Period (1000 to 1450 CE) to Late Horizon (1400 to 1540 CE in this region) sites surrounding the Coropuna volcano in southern Peru (hereafter COR) (Fig. 1A, Table 1, and dataset S1B). This region was an important ritual and pilgrimage center not only during but also before the Inca period, with traces of Wari and Tiwanaku influence (, ). Individual CO001, dated to 920 CE, comes from a pastoralist burial at Culcunche and most probably represents a local population that later hosted the ceremonial center of Maucallacta. Individual CO066 (1500 CE) originated from Maucallacta, whereas individual CO154 (1560 CE) was from Antaura, settlement located near Maucallacta. The last individual (CO193; 1320 CE) was a mummy from the nearby Cotahuasi Valley.

Genetic affinities—PCA and ADMIXTURE

We first performed a qualitative assessment of the genetic affinities of the individuals using principal components analysis (PCA) and unsupervised ADMIXTURE analyses. Principal components (PCs) were computed using data from present-day South American individuals without European ancestry. Genomic datasets for the present-day individuals were generated using different techniques, and their final intersection resulted in 199,175 SNPs in common (dataset S1D). Ancient individuals from this study and other available ancient South American individuals were projected onto the computed PCs (Fig. 2, A and B). PCA plot showed that most of the individuals clustered within their populations. We found a split between the Amazonian and the Andean populations, which corroborates the results of previous studies (, ). COR individuals from this study clustered together with other ancient Peruvian individuals, closest to those from the SouthernPeruvianHighlands (Laramate) and SouthernPeruvianCoast (Ullujaya and Palpa) ().

Fig. 2.

Genetic affinities of the ancient and present-day South American individuals.

(A) PCA plot of SNP-based genomic variation of 217 present-day South American individuals (empty markers) with ancient individuals projected (filled markers). Label for each marker is provided in (C). The covered part of this panel does not contain any data points. (B) Close-up of the most relevant section of the plot. (C) ADMIXTURE plot for K = 5. (D) MDS plot for the outgroup f3-statistic [1-f3(IND, Test, Mbuti)] matrix.

Present-day Bolivians made up a tight cluster with ancient individuals from the Lake Titicaca basin. However, two of the TIW individuals fell outside this ancient Titicaca cluster. TW056 was projected within the dispersed cluster of modern Amazonians, Columbians, and ancient individuals from the NorthernPeruvianCoast, whereas TW059 was projected between the ancient northern and southern Peruvian groups (Fig. 2, A and B). Unsupervised ADMIXTURE analysis showed that, for K = 5 [the lowest cross-validation (-CV) error; dataset S1E], two of the ancestral components were dominant and characteristic for two Amazonian populations from Brazil: Karitiana and Surui (red and blue, respectively; Fig. 2C). One component (yellow) was characteristic for Wichi, a population from the Gran Chaco (northern Argentina). The last two components were dominant in the Peruvian Amazonian (green) and Andean (violet) populations.

LUK individuals, together with other ancient Lake Titicaca individuals (ORU, Rio Uncallane, Iroco, and Miraflores), trace their ancestry mostly (>95%) to a single (violet) component. On the other hand, most of the TIW individuals showed at least some proportion of other components. The exceptions are TW004, TW0097, and I0977 () individuals, which comprise the violet component almost exclusively (Fig. 2C).

Ancestry modeling—f-statistics

To investigate the genetic makeup of the studied individuals in greater depth, we computed outgroup f3- and f4-statistics. In all calculations, Mbuti, a hunter-gatherer population from the Central Africa, was used as an “outgroup.” Shared genetic drift was measured between each of the 17 individuals from this study (“Ind”) and all other individuals/populations (“Test”) from the dataset in the form f3(Ind, Test, Mbuti). We generated a multidimensional scaling (MDS) plot and a neighbor-joining tree based on outgroup f3-statistics (Figs. 2D and 3). We performed exhaustive calculations of f4-statistics in the form f4(Mbuti, Population; Ind1, Ind2) and f4(Mbuti, Ind; Group1, Group2). Here, Ind represents an individual from this study, whereas Population or Group represents a population or a group of ancient individuals from a specified region and period as defined in () or a present-day population. Significantly negative Z score values (Z < −3) suggest that Ind1 shared more alleles with test Population than with Ind2 (dataset S2A) or with Group1 than with Group2 (dataset S2B).

Fig. 3.

Neighbor-joining tree based on inverted outgroup f3-statistics.

Distance: 1/f3(IND, Test, Mbuti). Individuals from the current study are shown in bold, TIW in red, LUK and ORU in blue, and COR in black. Previously published ancient individuals are marked with the colors used for their responding markers in Fig. 2C.

According to both f3- and f4-statistics, all LUK individuals showed the highest affinities with each other or with other ancient individuals from the Lake Titicaca area (Figs. 2D and 3 and dataset S2B). The homogeneity of the Lukurmata individuals was further tested using f4-statistics in the form f4(Mbuti, Test; LUKind1, LUKind2), where LUKind1 and LUKind2 represented each possible combination of Lukurmata individuals (dataset S2C). Neither of these tests, nor any in the form f4(Mbuti, Test; LUK, ancientTiticaca), were significant. This result shows that none of the Lukurmata individuals shares significantly more alleles with either Test group than with either individual from Lukurmata or the group of ancient individuals from the Titicaca area. We used qpWave/qpAdm, which summarizes multiple f4-statistics to determine the minimal number of distinct ancestry sources required for two tested populations, and found that, to the limits of our statistical resolution, the Lukurmata population was genetically homogeneous and, together with ancient Titicaca population, could be explained with a single wave of ancestry (Fig. 4 and dataset S3A).

Fig. 4.

Ancestry modeling for Tiwanaku site and Lukurmata individuals.

Ancestral source of the Tiwanaku (TIW) and Lukurmata (LUK) individuals, as estimated by qpWave and qpAdm using all possible South American populations as potential sources. Multiple columns for a single individual represent all of the potential single sources of ancestry that fit the data. For admixed individuals with more than one source of ancestry, we display significant results (P > 0.05) of the fitted two-source model from qpAdm (including SE). Provenience of each of the Tiwanaku ritual core remains is labeled with a silhouette of the corresponding monument. nd, not dated.

Using the same strategy, we found that the group of individuals from the Tiwanaku ritual core was heterogeneous, and we could recognize various patterns of ancestry (datasets S2 and S3). Individuals TW004, TW008, TW060, and TW097, as well as the published Akapana individuals I0977 and I0978 (), together with ancient Titicaca, could be modeled using one source of ancestry. Individuals TW059 and TW063 showed strong evidence of admixture, as they could not be modeled using a single source of ancestry and instead required at least two sources. The additional source of ancestry is likely Amazonian, because only models with ancient Titicaca and either Amazonian or Gran Chaco sources fit the data (dataset S3C). Individual TW061’s very low genomic coverage limited the resolution of our analysis, and we could not distinguish between ancient Titicaca and north Chile as a single source of ancestry (dataset S3B). The highest values of outgroup f3-statistic between TW061 and two north Chilean groups (Pukara and Pica Ocho) suggested that this individual might have had north Chile–related ancestry (Fig. 3 and fig. S4). Individual I0976 () could be modeled with a single source of ancestry, although we were not able to discriminate between southern Peru highlands [consistent with ()] and north Chile ancestry (dataset S3B). However, outgroup f3-statistics showed close affinities between this individual and the Peruvian groups (Figs. 2D and 3). The greatest outlier was TW056, which showed exclusively Amazonian-related ancestry, as indicated by both f4-statistics (dataset S2A) and qpWave (Fig. 4 and dataset S3B). Outgroup f3-statistics and f4-statistics indicated that the Peruvian Amazonian populations showed the greatest degree of shared genetic drift with this individual (Fig. 3 and dataset S2B). The distinctiveness of this individual was also expressed in the analysis of stable nitrogen (δ15N) and carbon (δ13C) isotopes, given that her profile was outlying from other Titicaca basin individuals tested here and previously (dataset S1C and Supplementary Text Information).

The single individual from Totocachi (TW033) was modeled with a single source of ancestry, although we were not able to discriminate between the SouthernPeruHighlands and NorthChile groups. The f4-statistics also showed that TW033 had greater affinity to SouthernPeruHighlands than to the ancientTiticaca group (fig. S5 and dataset S2F). We did not find any significant differentiation among the Coropuna individuals, as none of the f4-statistics computed in the form f4(Mbuti, AncientGroup, COR1, COR2) was significant (dataset S2G). Outgroup f3-statistics, f4-statistics, and qpWave modeling indicated closeness to South Peruvian populations (Fig. 3; datasets S2, B and H, and S3D; and fig. S5). Only individual CO066 could not be modeled with a single source, and it was best modeled as a combination of SouthernPeruCoast with a minor component from one of the groups to the north (dataset S3E).

DISCUSSION

Our results show that the population of the Titicaca basin was rather genetically homogeneous from at least 300 CE until the arrival of Europeans. This is well illustrated for the Lukurmata site, for which we have the widest time transect consisting of four individuals from between ca. 300 and 1500 CE. These individuals exhibit genetic similarity indicating that, over at least 12 centuries, no major genetic turnovers took place, despite the occurrence of substantial cultural and political changes. Earlier genetic studies of present-day South American populations, using either uniparental markers (, ) or genomic data (), suggested that both precontact and present-day Central Andean populations were genetically homogeneous and characterized by high levels of gene flow between them. In contrast, we find little overlap between precontact individuals from central/southern Peru and those from the Lake Titicaca basin (Figs. 2 and 3). Our results corroborate the findings in () that also suggested homogeneity within distinct regions of the Central Andes during the past ca. 2000 years. However, contact between various parts of the Andes was not entirely absent, as shown by the individual from the Coropuna region who had affinity to northern Peru and an individual from the Totocachi site who showed genetic affinity to the populations from the south Peruvian Andes or northern Chile.

Given the genetic homogeneity of the Bolivian Altiplano populations, the group of individuals from the Tiwanaku ritual core stands out because of its diversity. We found that several individuals show a close affinity to the Titicaca basin population, whereas others show distinct or admixed ancestry. At the peak of its influence, the Tiwanaku state reached outside the Titicaca basin to the fertile valleys of Moquegua (southern Peru), Cochabamba (Bolivia), and Azapa (northern Chile) (Fig. 1C) (, ). A pronounced presence of foreign ceramic clusters found at the Tiwanaku site supports a scenario in which ethnic groups settled in particular compounds around the complex (, ). Previous studies based on strontium isotopes () and cranial modifications () suggested the presence of individuals from neighboring areas, probably from Tiwanaku-affiliated sites at the peripheries of the polity. The evidence for contacts with populations from farther away, such as Amazonia, has been much more limited and based solely on isolated artifacts like a jaguar canine necklace (), iconographic representations of tropical lowland fauna (), tropical botanical remains (), and the consumption of lowland hallucinogenic substances (). The on-site presence of an individual with entirely Amazonian ancestry (TW056), although it is only one individual in the rather special circumstance of sacrificial offering, demonstrates that contact between the Tiwanaku site and Amazonia was not limited to the exchange of goods but included the physical movement of people. Contacts with remote regions were further confirmed by the finding of two individuals (TW059 and TW063) who showed a combination of Lake Titicaca basin and Amazonia and/or Gran Chaco ancestry. On the basis of the genomic data, we could not determine whether these individuals were already admixed incomers to Lake Titicaca basin or descendants of incomers and people with local ancestry. The latter scenario, in which the admixture between incomers and Titicaca individuals occurred locally, seems more probable, considering that no signs of foreign genetic ancestry have hitherto been found in Titicaca basin individuals outside the Tiwanaku ritual core.

A crucial element for understanding of the organization of the Tiwanaku society and the eventual abandonment of the site is the interpretation of numerous coeval human remains surrounding the Akapana platform. The accumulation of ritual offerings can be interpreted as the public destruction of captives and ancestors of vanquished peoples by an emerging elite, as exemplified by the historically described sacrifices in the Aztec capital after successful wars of conquest (). This model carries the expectation that the sample will contain victims from various origins. TIW dataset consisted of eight individuals associated with Akapana, five of whom had confirmed context at the base of the platform. All but one (I0976) of them trace their full ancestry to the Titicaca basin, although the genetic homogeneity of Titicaca basin populations makes it impossible to track the provenance of these individuals precisely within the region. However, it is worth noting that they differ from individuals excavated at other parts of the ritual core, who show traces of more distant ancestries from outside the Titicaca basin. Direct radiocarbon dating puts most of the Akapana-associated individuals toward the end of the Tiwanaku phenomenon in the mid-10th century. This dating is too late to support the use of organized violence as a means of establishing and consolidating a hierarchy. Placed on alluvia that had accumulated over a prepared surface associated with the base of the platform, they present a pattern more reminiscent of individual visitations and offerings to local sacred sites known as huacas (). Perhaps, as the sociopolitical organization associated with the monuments lost power and/or the elites at Tiwanaku no longer had the capacity to hold big formal feasts and clean up afterward, groups of people could leave offerings at their own discretion.

There are documented cases where the intensification of human sacrifice, such as that observed at Tiwanaku, indicates a society that is in crisis and grasping for a solution (). The manner and process of the abandonment of the site of Tiwanaku are a subject of debate, but some scholars consider the last portion of the 10th century as the beginning of a gradual or sharp slide toward the abandonment of the site (, ). Some ascribe the end of the Tiwanaku site and, later, Tiwanaku culture to a long-term drought that started ca. 1100 CE, collapsing the raised field agricultural system (, ). Others find this explanation too environmentally deterministic and point to the fact that there is little evidence of Tiwanaku presence after 1000 CE, although the raised fields persisted for several centuries ().

The new genomic and radiocarbon data reported here demonstrate that, in the mostly static genetic landscape of the Central Andes, in which each region seems to have maintained an inner homogeneity with little exchange with outside populations, people from both neighboring valleys and places much further afield—even outside the direct influence of the polity—found their way to the Tiwanaku site. Although we may never know the reason why certain individuals took or were taken for such a journey, our data suggest that they, and/or their descendants, left a trace at Tiwanaku. The contrasting local ancestry of the human remains associated with Akapana is an intriguing finding. It remains uncertain whether this offering selection was based on some provenance-based social stratification within the Tiwanaku polity, a result of a change in its sociopolitical organization, or whether the offered individuals were members of some adversarial political power within the Titicaca basin. Resolving this question would require more genetic, isotopic, and archaeological data from within the Tiwanaku site, as well as from Tiwanaku-influenced regions.

MATERIALS AND METHODS

DNA extraction and library preparation

Samples from a total of 93 individuals from modern-day Bolivia and Peru were selected for genetic analyses. These human remains originated from archaeological sites with contexts and chronologies relating to Tiwanaku or Inca cultures (dataset S1A).

All work with human remains, DNA extraction, and library preparation were performed in facilities dedicated to working with ancient DNA. The laboratory was ultraviolet (UV)–irradiated when not in use. We followed all protocols recommended to prevent sample contamination with modern DNA. DNA was extracted from teeth or long bones. The samples were washed with ultrapure water and UV-irradiated (245 nm) for 10 min on each side in a laminar flow cabinet and powdered in a cryogenic mill (SPEX CentriPrep). In addition, for some teeth, we tried drilling only cementum as described in (). We performed extraction following the protocol commonly used for retrieving short fragments of ancient DNA (). Extractions were performed for up to 15 samples accompanied by negative controls.

Double-indexed sequencing libraries were constructed using 20 μl of DNA extract according to the protocol proposed in () with the following minor modification: After blunt-end repair and fill-in steps, the enzymes were heat-inactivated for 20 min at 75° and 80°C, respectively. Every library had a unique pair of P7 and P5 indexes. The number of cycles in the indexing PCR was determined using quantitative polymerase chain reaction (qPCR) with primers I7 and I8 (). Indexing PCR was performed using PfuTurbo Cx DNA polymerase (Agilent) and 10 μl of DNA. Three independent PCRs were performed for each library to increase complexity. PCR products were pooled and purified using SPRI beads. For a few samples, we performed whole-genome enrichment using the myBaits WGE Human Kit (Arbor Biosciences) according to the manufacturer’s protocol (dataset S1B). Sequencing libraries were pooled in equimolar ratios and were sequenced on Illumina platforms: NextSeq550 (HighOutput, 75 cycles, single-end; MidOutput, 150 cycles, paired-end), HiSeq4000 (50 cycles, single-end), or NovaSeq6000 (S1, 100 cycles, single-end).

Radiocarbon dating

Radiocarbon dating was performed at the University of Waikato Radiocarbon Dating Laboratory (Wk) and the Poznan Radiocarbon Laboratory (dataset S1C). The quality and purity of the extracted collagen were assessed using its chemical composition (%N, %C, and C:N ratio). In the case of three individuals (TW020, TW027, and TW098), the amount of extracted collagen was sufficient only to perform accelerator measurement but not assessment of the chemical composition of collagen nor measurement of stable isotope ratios. For these three individuals, the sole criterion of collagen quality was the gelatin yield, and this was above 0.5% in all three cases. The %C and %N values were in the accepted ranges for well-preserved ancient collagen () in all cases. The C:N ratio was in the acceptable range (2.9 to 3.6) in all but one individual (3.66 in CO193). In the case of this individual, we used a phalanx with mummified tissue for dating. The laboratory (Wk) reported that the slightly elevated C:N ratio may have resulted from the incomplete removal of tissue fragments during pretreatment. All dates were corrected for isotopic fractionation using accelerator mass spectrometry (AMS)–measured δ13C values.

We measured the carbon and nitrogen isotopic ratios (δ13C and δ15N) to obtain information on the diet of the individuals studied (Supplementary Information Text and fig. S1). Measurements were made using isotope-ratio mass spectrometry at three laboratories (University of Waikato Radiocarbon Dating Laboratory; Stable Isotope Laboratory at the Institute of Geological Sciences, Polish Academy of Sciences; and Institute of Geosciences at Goethe University). We did not find evidence of substantial consumption of freshwater fish, which could affect the calibration of radiocarbon dates. For calibration, we used the SHCal13 curve () in OxCal 4.3.2 (fig. S2) (). In addition, because of air currents that probably affect the carbon concentrations by carrying air from the Northern Hemisphere into the Titicaca region (), we calibrated the dates using a mixed-curve model (dataset S1C). The median ages of three samples differed by more than 10 years and up to 60 years between the SHCal13 and the mixed-curve calibration. Throughout this study, we report SHCal13-calibrated dates, rounded to the nearest 10 years.

Data processing

Adapter sequences were trimmed, and paired-end reads were collapsed using AdapterRemoval2 (). The sequencing reads were mapped to human reference genome h37d5, applying the default parameters of the Burrows-Wheeler Aligner (BWA) mem algorithm (). We used SAMtools () to remove duplicates. Only reads longer than 30 base pairs (bp) and with mapping quality over 30 were used in subsequent analyses. To minimize the impact of deamination on genotyping, we trimmed 7 bp from both ends of all reads using the trimBam tool from bamUtils ().

Authentication

We used mapDamage 2.0 () to assess the damage and fragmentation patterns of the obtained sequencing reads. Contamination with present-day human DNA was estimated using schmutzi () and contamMix (). In addition, for male individuals, we investigated nuclear contamination based on polymorphic sites on the X chromosome using ANGSD (). We set the minimum base quality to 30 and the minimum mapping quality to 30; other parameters were left at the default. Last, for both females and males, we estimated autosomal contamination using ContamLD version 1.0 () (estimate based on breakdown of the linkage disequilibrium). For ContamLD analyses, default settings were used with PEL as a reference population for all samples. All contamination results are reported in dataset S1B.

Sex determination

The genetic sex of each individual was determined by calculating the ratio of sequence reads aligning to the X and Y chromosomes, as described in (). We also compared the coverage of the X and Y chromosomes with the coverage of autosomal chromosomes, as described in ().

Genotyping and datasets

Genotypes were called by choosing one random read for every SNP from the 1240K SNP list () using the script pileupCaller, a part of sequenceTools (https://github.com/stschiff/sequenceTools). The analyzed individuals yielded between 142,180 and 926,465 SNPs overlapping with the 1240K set (dataset S1B). These data were merged with the available modern and ancient genome data from South America (dataset S1D) (, , –). Because present-day genotypes were obtained using various techniques, including shotgun sequencing and different Illumina microarrays, the final intersection resulted in 199,175 SNPs in common. The number of SNPs intersected in the samples from this study ranged from 24,601 to 161,748 (dataset S1B). In downstream analyses, we used published ancient genomes that yielded more than 20,000 SNPs and were dated to the past 2000 years (dataset S1D). The exception was an individual I0977 from Tiwanaku () with 11,710 overlapping SNPs.

Genetic affinities

We computed PCA using the smartpca script from the EIGENSOFT package () applying lsqproject=YES and shrinkmode=YES options. Ancient individuals were projected on a PC plot computed using 37 modern unadmixed South American populations (dataset S1D). To investigate the genetic structure of the Tiwanaku populations, we performed an unsupervised admixture analysis using ADMIXTURE (). Before the ADMIXTURE analysis, we pruned genotypes for minor allele frequency below 0.01 (--maf 0.01) and linkage disequilibrium using a window size of 200, a step size of 5, and an R2 threshold of 0.5 using plink (--indep-pairwise 200 5 0.5) (). Following these pruning steps, we retained 100,290 SNPs. Every individual from this study met the requirement of minimum 10,000 SNPs and was used in this analysis. We ran five replicates for each value of K (K = 2 to K = 15), and the optimal K was chosen on the basis of the lowest cross-validation error.

Ancestry modeling

Outgroup f3-statistics were estimated using qp3pop from ADMIXTOOLS () in the form f3(Ind, Test, Mbuti). In all calculations, we ascribed each of the individuals analyzed in this study as Ind and every other individual/population from the dataset as Test (dataset S1D). Higher f3-value means a higher genetic affinity between two tested individuals or between an individual and a population. All obtained values were plotted using the ggplot2 package in R (fig. S4). We created a dissimilarity matrix from the f3-statistics, calculating the 1-f3 values. We performed MDS using the cmdscale and plotted the principal coordinates. We generated a neighbor-joining tree using PHYLIP () and used the ancient Beringian individual USA_USR1_AncientBeringian_11400BP.SG as an outgroup ().

The R package admixR (), which uses the ADMIXTOOLS software suite, was used to calculate f4-statistics and for qpWave and qpAdm modeling. SEs were computed using a jackknife block size of 0.050. We used qpWave (version: 600) to determine the minimum number of ancestry sources for each individual and groups of individuals using ancient and modern populations (dataset S1D). The same populations were used as a left (“Source”) in qpAdm modeling (version: 1000, allsnps:YES). A “rotating” strategy was used to set the right (“Reference”) populations, as suggested in (). In this strategy, the basic set of reference populations included Argentina_ArroyoSeco2_7700BP, Peru_Cuncaicha_4200BP, Peru_Lauricocha_3500BP, Brazil_LapaDoSanto_9600BP, and Wichi as well as all the populations from the Source list, except for the one population used as a Source in the particular estimation. Under this rotating approach, populations are consistently moved from the Source set to the set of Reference populations.

Mitochondrial DNA and Y-chromosome haplotyping

For every studied individual, we obtained a complete sequence of mitochondrial genome. The sequencing reads were mapped to a reference mitochondrial DNA sequence (rCRS, NC_012920) using BWA mem. Several programs from SAMtools () were used to remove sequences shorter than 30 bp that had mapping quality under 30 and that were duplicates. Variants and consensus sequences were called using bcftools (). Positions with coverage below 3 were masked. BAM files were manually checked using Tablet (). Nucleotide positions 309 to 311, 523 to 524, 3107, 8271 to 8279, 16,182, and 16183 were masked and omitted in haplogroup calls. The obtained mitochondrial genomes were assigned to haplogroups using Haplogrep2 (Phylotree 17) (). Y-chromosome haplogroups were determined using Yleaf v2.1 ().