Ancient DNA Analysis Suggests Negligible Impact of the Wari Empire Expansion in Peru’s Central Coast during the Middle Horizon

The analysis of ancient human DNA from South America allows the exploration of pre-Columbian population history through time and to directly test hypotheses about cultural and demographic evolution. The Middle Horizon (650–1100 AD) represents a major transitional period in the Central Andes, which is associated with the development and expansion of ancient Andean empires such as Wari and Tiwanaku. These empires facilitated a series of interregional interactions and socio-political changes, which likely played an important role in shaping the region’s demographic and cultural profiles. We analyzed individuals from three successive pre-Columbian cultures present at the Huaca Pucllana archaeological site in Lima, Peru: Lima (Early Intermediate Period, 500–700 AD), Wari (Middle Horizon, 800–1000 AD) and Ychsma (Late Intermediate Period, 1000–1450 AD). We sequenced 34 complete mitochondrial genomes to investigate the potential genetic impact of the Wari Empire in the Central Coast of Peru. The results indicate that genetic diversity shifted only slightly through time, ruling out a complete population discontinuity or replacement driven by the Wari imperialist hegemony, at least in the region around present-day Lima. However, we caution that the very subtle genetic contribution of Wari imperialism at the particular Huaca Pucllana archaeological site might not be representative for the entire Wari territory in the Peruvian Central Coast.


Introduction
The Central Andes in South America are defined by marked geographic contrasts between the Andean highlands and the Pacific coast, with each of the regions being characterized by unique

Mitochondrial DNA haplogroups
From a total of 52 individuals collected from the Huaca Pucllana site, we obtained 34 complete mtDNA genome sequences (9 for Lima, 10 for Wari, 15 for Ychsma) using hybridization capture with RNA baits in solution and Illumina high-throughput sequencing. The average coverage depth per position was 97.4x (range: 16.8x-177.9x; Table 1). All 34 individual sequences encompass 28 distinct haplotypes, and were assigned to one of the main four Native American "founder" mtDNA haplogroups A2, B2, C1, and D4/D1 [28,29]. A remarkable result of our study is that all ancient mtDNA sequences from Huaca Pucllana represent novel haplotypes not hitherto observed in modern-day populations, which harbor many additional variants from the common sub-haplogroups A2, B2, C1 and D4/D1 (Table 1).
Detailed sub-haplogroup frequencies for the populations analyzed in this study were calculated per cultural period and are reported in Fig 2. We then used a Fisher exact test to investigate whether the sub-haplogroup composition varied significantly between the three cultural groups. However, all three comparisons returned non-significant p-values (Table 2).

Summary statistics and haplotype sharing
Negative sequence-based genetic distances pairwise F ST between the tree cultural strata Lima, Wari and Ychsma suggest genetically homogenous groups through time: -0.05197 for Lima-Wari, -0.06111 for Lima-Ychsma and -0.00770 for Wari-Ychsma (associated p-values were also non-significant after Bonferroni correction; Table 2). Haplotype diversity (h) also did not differ between the three cultures: Lima (0.9722), Wari (0.9778), and Ychsma (0.9810). Likewise, we did not observe substantial differences for nucleotide diversity (π): Lima (0.002331), Wari (0.002492) and Ychsma (0.002159). Genetic diversity indices and neutrality tests indicated that individuals of the Ychsma culture had a higher haplotype diversity 0.9810 (p>0.0308) than the two previous occupation phases. Values for Tajima's D (neutrality test) only suggested a population expansion for individuals associated with the Lima culture (-0.13498), while for Fu's FS values, neither of the three groups Lima (1.44409), Wari (1.08278) and Ychsma (0.31237) indicated a recent population expansion (Table 2).

Huaca Pucllana mtDNA diversity in the South American context
We evaluated the relationship of cultural groups from Huaca Pucllana with modern and ancient mtDNA diversity in South America based on a genetic distance matrix and visualized as a Multidimensional Scaling Plot (MDS).
Individuals from the EIP Lima culture are genetically close to highland Laramate populations from the MH [32] and LIP [33], and to Palpa coastal samples from the LIP [32]. Interestingly, ancient Lima differed completely from the southern coastal Nasca populations (Palpa EIP) from the same cultural horizon. Moreover, individuals from the MH Wari culture from Conchopata in the highlands also fall close to the coastal Lima individuals, whereas the Wari  (Fig 4).

Discussion
We sequenced 34 complete human mitochondrial genomes from the three successive cultural periods Lima, Wari and Ychsma found at Huaca Pucllana. This dataset represents a unique chronological transect in the Peruvian Central Coast. Low and non-significant F ST values for comparisons between the three groups Lima, Wari and Ychsma indicate no differentiation of sequence haplotypes, meaning that ancient Peruvians from Huaca Pucllana did not undergo dramatic changes in genetic ancestry over approxi-mately~1000 years. These results suggest population continuity through time, with no major Additional variants from the sub-haplogroup nodes compared to RSRS [30]. We generally follow the nomenclature of Phylotree.org [31], in describing variants: Reversions (back mutations) to an ancestral state are indicated with an exclamation mark (!), while 7472d and 356.1C denotes a deletion or insertion at the given nucleotide position, respectively.
doi:10.1371/journal.pone.0155508.t001 demographic turnovers. Genetic diversity indices do not differ substantially between the three cultural periods and results from neutrality tests (Tajima's D and Fu's FS) do not attest for recent expansion in any of the three groups. Overall, samples from the latest Ychsma period are the most diverse, which could be explained by the slightly larger sample size ( Table 2).  However, mitochondrial sub-haplogroup frequencies reveal subtle shifts between the three cultural periods. We observed a loss of basal haplogroup B2 lineages during the transition from Lima to Wari, alongside a slight increase of B2b lineages, suggesting the arrival of the latter in the Huaca Pucllana region during the occupation phase of the Wari Empire. Interestingly, the fact that basal B2 lineages (non-B2b) are present in subsequent LIP Ychsma individuals, while being absent in MH Wari, suggests a potentially exogenous source of the Wari B2b lineages (Fig 2).
In addition, half of the MH Wari individuals carried the C1b sub-haplogroup, which is twice as much as the earlier EIP Lima or later LIP Ychsma groups. Since the frequency of haplogroup B2b remains constant during the Wari and Ychsma periods, we suggest that the increase in frequency of haplogroups C1b and B2b during the Middle Horizon might be linked to gene flow from the Andean region. However, we note that haplogroup A2 lineages also increase with the Wari period, but do not seem to become established in the Huaca Pucllana region during Ychsma times, following the collapse of the Wari Empire (Fig 2).
We caution that haplogroup composition represents a very broad approach to investigate signals of population structure, in particular because frequencies are subject to genetic drift or stochastic processes, population movements and gene flow [37]. Moreover, sample assignation was confirmed by radiocarbon dates for only a subset of specimens from the three cultural periods. It is therefore possible that some of the undated specimens may have been misassigned to a particular cultural layer and results may be interpreted with care. Network analysis showed that haplotype sharing occurred mostly between individuals from the same archaeological period, and less often between different cultural layers (only two instances, i.e. Wari-Ychsma, Lima-Ychsma). A shared C1b lineage between Wari and the subsequent Ychsma is very plausible, while the shared Lima-Ychsma haplotypes suggests a local continuity of particular maternal lineages over several centuries, which implies that the Wari incursion was not a replacement but rather brought new lineages to the Huaca Pucllana region (Fig 3).
Interestingly, the MDS plot shows an affinity of the Lima population during the EIP (i.e. before the arrival of Wari) in Peru's Central Coastal with ancient highland populations rather than other Coastal populations from the same period (i.e. Palpa EIP and Montegrande EIP), suggesting that the population history of the Central Coast is connected with the highlands and did not develop independently unlike Moche and Nasca populations [38]. However, this observation is independent from the archaeological record of the Central Coast and should not be over-interpreted in terms of cultural origins of the Lima Culture. The clustering of the Wari individuals from Huaca Pucllana with Huari-LIP from the highlands suggests a genetic similarity of the Wari group from the MH with the post-Wari LIP individuals in the Central Andean region in Peru. The persistence of Wari lineages in individuals from the subsequent Ychsma culture explains the intermediate position of the Wari between the Ychsma and the Lima substratum and could be considered as legacy of the Wari influence (Fig 4).
It should be noted that the overall mtDNA haplogroup composition is nonetheless homogenous across cultural periods in Peru Central Coast, even in the light of the expansion of the Wari Empire. To the limits of our resolution, the new results from ancient Huaca Pucllana therefore hint at a very subtle yet negligible genetic contribution of highlands Wari in coastal Peru during the MH.
We caution that the situation at Huaca Pucllana in the greater Lima region might not be representative of the overall genetic influx of the Wari Empire in the Central Coast of Peru, supporting the hypothesis that local settlements along the Central Coast were not entirely dominated or affected by Wari imperialism [23,24]. Moreover, it appears that the Wari did not adopt highly aggressive means of control such as direct colonization outside the political center in Ayacucho, as suggested previously [39]. In this context, we suggest that the independence of local and regional identities during the MH must have remained intact, and thus implies group interactions and a level of political and social complexity which has not yet been fully explored in the Central Coast of Peru [39]. Consequently, definitive assertions about the Wari imperialism/colonialism in Peru's Central Coast require further research.
A remarkable observation is that the ancient mtDNA sequences from Huaca Pucllana represent exclusively novel haplotypes (Table 1) with many variants additional to the common subhaplogroups A2, B2, C1 and D4/D1. This suggests that the pre-Columbian genetic diversity was previously underestimated, which can be explained by two complementary premises: i) the mtDNA diversity of pre-Columbian populations was substantially higher than in modern-day Native populations; and ii) there was a considerable loss of pre-Columbian mtDNA diversity as a result of the dramatic demise of Native Americans during and after the European colonization [40]. Alternatively, it is also possible that the modern genetic diversity is largely underestimated, and that sequencing mitochondrial genomes from relevant Central Andean populations might reveal direct matches with samples from colonial times, and today.
As exemplified by the pre-Columbian expansions of Wari and Tiwanaku and later on the Inca Empire, the Central Andes have always been a dynamic environment with considerable population movements [1]. Another large-scale internal migration movement took place in Peru in more recent times, with rural and indigenous population (e.g. Junín, Ayacucho, La Libertad, Ica, Lambayeque, Cajamarca, Piura, and to a lesser degree from other places) moving towards the urban centers like Lima since the 1940's [41]. It is therefore likely that the ancestry of people from large urban centers like today's Lima must have been constantly receiving additional genetic diversity.

Conclusion
The~1000-year transect through time of pre-Columbian history at Huaca Pucllana, allowed us to test the impact of the expansion of the highland Wari Empire on the coastal Lima culture, and its legacy in the subsequent Ychsma culture. The results from three archaeological periods indicate that genetic diversity in Huaca Pucllana shifted only slightly through time, ruling out a complete population discontinuity or replacement driven by the Wari imperialist hegemony, at least in the region around present-day Lima.
Detailed chronological sampling and high-resolution DNA sequencing efforts are required to refine the analysis of ancient genetic diversity in Peru, but importantly also from modern-day groups from Central Andean and Coastal regions to be able to contrast ancient diversity with present-day populations. Such data would allow the complete contextualization of associated meta-data, testing of demographic scenarios, and exploration of continuity/discontinuity, replacements and bottlenecks on the basis of statistical inferences or simulations. Ultimately, these analyses could elucidate the dynamic of pre-Columbian societies and the dramatic demise of Native American people that coincided with the time of the European arrival [42,43].

Materials and Methods
Ethics Statement: All necessary permits were obtained for the described study, which complied with all relevant regulations. Permissions to collect, export and analyze ancient Peruvian specimens from the Huaca
Archaeological evidence at the site defined that three cultural layers developed during slightly different times than in the classical chronology for the Central Andes established by Rowe [44]. The Lima culture (500-700 AD) occupied the site during the late EIP and early MH. The Wari culture established as a state in Ayacucho~600 AD but reached the Central Coast and developed at Huaca Pucllana only later, around 800-1000 AD [45]. Finally, the time of the Ychsma culture (1000-1450 AD) is concordant with the classical chronology [46].
The construction of the Huaca Pucllana ceremonial center was done between the 4 th and 7 th centuries AD, which represents the Lima final phase (500-700 AD) and occupation of the site. Huaca Pucllana is characterized by pyramidal platforms built with sun-dried bricks, or adobes [46,47]. The main structure at the site is the Great Pyramid, which is 400 m long, 60 m wide, 20 m high, and includes seven platforms [48][49][50]. It was used for various activities related to the worship of divinities associated with the so-called Sea Culture, with constant pilgrimage of people from surrounding regions. Human remains excavated in Huaca Pucllana were often associated with cultural elements that could be attributed to a particular layer and cultural period. In addition, remains of ceremonial activities, such as offerings of pottery vessels, ritual banquets, were also identified [51].
The Lima culture developed funerary patterns characterized either by burials, whereby single or multiple bodies were wrapped with cloth in extended position, or the practice of human sacrifices with no specific position of the bodies, with or without associated pottery and burial offerings [51]. Bioanthropological analysis of human remains at Huaca Pucllana showed a large number of violent deaths among the Lima individuals [51]. Notably, some individuals (i.e. young women), were sacrificed as part of rituals during the last stages of the construction of the Great Pyramid [52]. Huaca Pucllana appears to have been abandoned in the 7 th century AD as a result of political and population changes driven by the expansion of the Wari, who arrived from the highland region of Ayacucho [8].
The Wari culture represents the second settlement stratum at Huaca Pucllana during the MH (800-1000 AD). Wari people partly destroyed and re-structured the former ceremonial Lima center at the Great Pyramid, and used the place as an imperial elite cemetery [45]. However, most of the Wari burials were destroyed later on, which altered the distribution of the elements of each burial, and intact mummies were only found in very few cases [49]. Most of the burials from this stratum are characterized by "fardos". These funerary bundles preserved corpses before burial and are the main distinctive feature found in the upper platforms of the Great Pyramid (Fig 5).
Other diagnostic signs of elite Wari mortuary practices are the presence of clothing, household goods, ritual items and food [53,54]. Some tombs and burials only contain a single individual, while others can contain more than four bodies. Some cases include human sacrifices such as young children next to the main body, which is described as a typical Wari mortuary practice for important personalities as recently confirmed at the "Castillo de Huarmey" archaeological site [55].
The major cultural transition from the MH to the LIP (1100-1476 AD) involves the collapse of the Wari Empire [2,3,56]. From~1000 AD onward, the Wari power gradually declined as their political center increasingly lost control over cities and the vast Central Andean territories. During this period of crisis, large urban centers such as Cajamarquilla and Pachacamac in the Central Coast were abandoned [21], only to become occupied by the Ychsma culture (1000-1450 AD) soon after.
However, there is no archaeological evidence for an immediate Ychsma occupation at Huaca Pucllana. Instead, some remaining local people (farmers and fishermen) with no particular cultural affiliation settled near the ceremonial center for about 100 years, later to become the social unit called Ychsma. The Ychsma culture marks the third major period at Huaca Pucllana, during which the site was again used for offerings and burials. The Ychsma culture is associated with a characteristic pottery closely associated to funerary patterns that also used the "fardos" as mortuary elements. Similarly to other settlements/locations in the Central Coast, e.g. Pachacamac [57], several elements suggest a reconfiguration of Huaca Pucllana as religious space that was then continuously used during Ychsma times [46].

Sample description
A total of 115 samples-mostly teeth-from 52 individuals were collected at the Museo del Sitio Huaca Pucllana (S1 Table). Individuals could be associated to either one of the three successive cultures found in Huaca Pucllana: Lima (n = 35), Wari (n = 47) and Ychsma (n = 33) [46,58]. We followed standard aDNA sampling guidelines, and collected at least two independent samples per individual for replication and authentication purposes [59,60].

Radiocarbon dates from Huaca Pucllana samples
Nine samples excavated at the site were sent for radiocarbon analysis at the University of Oxford Radiocarbon Accelerator Unit, England. The dates were then used to confirm the assignment to cultural periods based on archaeological context for all other Huaca Pucllana samples. Radiocarbon dates were calibrated using the Oxcal computer program (v4.2), applying the`IntCal13' dataset-Northern Hemisphere [61]. Dates estimates were also calibrated with the Southern Hemisphere correction (SHCal13) [62]. For each individual sample analyzed, we observed slight differences between dates calibrated with IntCal13 and SHCal13, although there was a strong overlap. However, these differences between calibration methods did not lead to ambiguities regarding cultural assignment (Table 3).

Sample preparation and DNA extraction
Samples were processed in a dedicated aDNA facility at the University of Adelaide's Australian Centre for Ancient DNA (ACAD), Australia. The laboratory employs standardized protocols and infrastructure for aDNA analysis (positive air pressure, UV irradiation and regular cleaning with oxidating agents, e.g. commercial bleach and Decon 1 to minimize contamination) [63][64][65]. Samples were decontaminated upon entry in the aDNA laboratory by exposure to UV light. The surface of the samples was gently wiped with 3% bleach and then physically removed by abrasion using a Dremel 1 drill. A Mikro-dismembrator ball mill (Sartorius) was used to pulverize the sample and 0.2 g of bone powder were subsequently used in DNA extractions. For the extractions, samples were decalcified by incubation in 4 mL of 0.5 M EDTA (pH 8.0) overnight at 37°C and mixed constantly on a rotor mixer. Next, 70 μL Proteinase K (Invitrogen) was added and the lysis mix was incubated for 2 hours at 55°C under constant rotation. DNA was isolated using silicon dioxide solubilized in a Guanidinium buffer (Qiagen), as described previously [66,67]. DNA was resuspended in 200μL of TE buffer including 0.05% Tween-20 and stored at -20°C until further use.
Library preparation and hybridization capture of mtDNA All DNA extracts were initially screened by sequencing the mitochondrial Hypervariable Region (HVR-I), and by genotyping using two multiplex mitochondrial SNP assays, the Geno-CoRe22 [68] and the AmericaPlex26 [69]. We selected one DNA extract for genomic libraries preparation from all individuals for which replications of HVR-I sequences and SNP genotypes were successful. DNA libraries and enzymatic clean-up steps were performed in a dedicated contamination-free room at ACAD. Library preparation with barcoded truncated adapters and hybridization capture of mtDNA are described in [27].

DNA sequencing and sequence assembly
Pooled DNA libraries were sequenced on one Illumina flow cell lane on an Illumina HiSeq2000 at the Australian Cancer Research Foundation (ACRF) Cancer Genomics Facility, Adelaide, South Australia. Post-sequencing processing was performed using a custom pipeline. Raw Illumina reads were processed and filtered by means of demultiplexing of DNA sequences according to their index sequence, using the Illumina program Casava. Further demultiplexing and trimming of internal barcodes was performed using the program Sabre 1.0 (https://github. com/najoshi/sabre). Residual adapter sequences were trimmed and read pairs were merged using Adapter Removal v1.5 [70] with default parameters. Merged reads were mapped to the reference sequence RSRS (Revised Sapiens Reference Sequences) [30] using the mapping program bwa 0.7.5a-r405 [71], and duplicate reads were removed using FilterUniqueSAMCons.py [72]. Ancient DNA damage patterns were assessed using MapDamage v0.3.6 [73].

Sequence analysis
Read pileups were visually inspected using the Geneious Pro 1 Software V.6 (Biomatters Ltd) [74] and variants were called against the Reconstructed Sapiens Reference Sequence RSRS (16,569 bp) [30]. Variant calling was based on a minimum coverage of 3x and an initial automated majority call of 75%, followed by independent verification by eye by two researchers (GV and BL) from the read pile-ups to avoid errors due to DNA damage or to locally re-align regions around homopolymeric stretches.
Haplogroup determination was performed using the online database phylotree (mtDNA tree Build 16 [19 Feb 2014]) [31]. We excluded substitutions at nucleotide positions 16182 and 16183, because they are dependent on the presence of C at the position 16189 [75]. For phylogenetic reconstruction and network analysis, we did not consider positions 309.1C(C), 315.1C AC indels at 523 and 524, deletion 3107 and position 16519 according to [31]. The resulting consensus sequences from Huaca Pucllana individuals were also included in a parallel study surveying the wider pre-Columbian mitochondrial genome diversity across South America [43].

Population genetic and statistical analyses
Haplotype diversity (h), nuclear diversity (π), pairwise genetic distances F ST , neutrality test Tajima's D [76] and Fu's FS [77] were calculated using the software Arlequin 3.5 [78]. P-values of genetic distances were corrected for multiple comparisons using the p.adjust function in R 3.2.3 (The R Project for Statistical Computing 2011, https://www.r-project.org/). Median Joining (MJ) networks [79] were built using the program Network version 4.6 (http://www.fluxusengineering.com) to evaluate genetic relationships among all individuals and to illustrate the extent of haplotype sharing within and among temporal periods and cultures.
We generated multidimensional scaling (MDS) plots to visualize genetic similarities and dissimilarities as measured by fixation indices (Slatkin's F ST ) [80] in bi-dimensional space using a customized R script (https://www.r-project.org/). Since the mtDNA genome data for modern-day populations is heavily biased towards phylogenetic and phylogeographic questions that result in an excess of specific haplogroups, we had to restrict our analyses to mtDNA (HVR-I) sequences compiled from published sources to describe the genetic relationships between ancient samples from Huaca Pucllana and other ancient individuals (n = 352) from 16 populations, as well as present-day South Americans (n = 1090) from 28 populations (S2 Table).
We used a Fisher exact test to investigate the variation in sub-haplogroup composition between the three cultural groups. Individuals were assigned to sub-haplogroups A2, B2, B2b, C1b, C1c, C1d, D1, and D4/PreD1/D1. The Fisher extact test was computed using the function fisher.test in R 3.2.3 (http://www.r-project.org) and 10,000 permutations per comparison.
Supporting Information S1