Prehistoric mitochondrial DNA of domesticate animals supports a 13th century exodus from the northern US southwest

The 13th century Puebloan depopulation of the Four Corners region of the US Southwest is an iconic episode in world prehistory. Studies of its causes, as well as its consequences, have a bearing not only on archaeological method and theory, but also social responses to climate change, the sociology of social movements, and contemporary patterns of cultural diversity. Previous research has debated the demographic scale, destinations, and impacts of Four Corners migrants. Much of this uncertainty stems from the substantial differences in material culture between the Four Corners vs. hypothesized destination areas. Comparable biological evidence has been difficult to obtain due to the complete departure of farmers from the Four Corners in the 13th century CE and restrictions on sampling human remains. As an alternative, patterns of genetic variation among domesticated species were used to address the role of migration in this collapse. We collected mitochondrial haplotypic data from dog (Canis lupus familiaris) and turkey (Meleagris gallopavo) remains from archaeological sites in the most densely-populated portion of the Four Corners region, and the most commonly proposed destination area for that population under migration scenarios. Results are consistent with a large-scale migration of humans, accompanied by their domestic turkeys, during the 13th century CE. These results support scenarios that suggest contemporary Pueblo peoples of the Northern Rio Grande are biological and cultural descendants of Four Corners populations.


Introduction
Although archaeologists have long recognized that maize farmers left the northern US Southwest by the late 1200s CE [1] the identities of the emigrating populations, their destinations, and the causes of the depopulation have been controversial. Recent work by the Village Ecodynamics Project (VEP) in the central Mesa Verde region (CMV) (Fig 1), the most populous portion of this vast area, has shown that population size peaked in the mid-1200s, was decreasing by no later than 1260 CE, and declined to zero by about 1280 CE [2].
Additional research [3,4] has used evidence from biometry, linguistic prehistory, site architecture, etymology, and recorded social memory to argue that the ancestral Tewa population of the Northern Rio Grande (NRG) in New Mexico derives primarily (but not exclusively) from a migration of Pueblo III farmers from the CMV (Fig 1). However, discontinuities in many aspects of material culture between the pre-1280 CE CMV and post-1280 CE NRG [5] have led some to diminish both the size of the CMV influx and the importance of these immigrants for subsequent culture processes in the NRG [6,7].
Previously collected, comparative mtDNA data from two turkeys recovered from Forked Lightning (LA672, 1175-1350 CE) were used in our analyses [22] (Table 1). Note that Speller et al. [22] dated this site to 1300-1840 CE, but we associate these samples with the pre-migration period based on the tree-ring data and pottery from this site [23:42].

DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing
All pre-polymerase chain reaction (PCR) activities were conducted in the Ancient DNA Laboratory at Washington State University (WSU), a laboratory located in a separate building from wherein PCR and post-PCR activities are conducted. The Ancient DNA Laboratory is a dedicated workspace for processing degraded and low copy number (LCN) DNA samples. Precautions aimed to minimize the introduction of contamination are practiced in the laboratory [24]. DNA was extracted from the turkey and canid remains by one or more of the following three methods. These methods were employed to aid in maximizing the success rate of recovering genetic information from these remains (e.g., [25]). Under the first extraction method (noted as Extraction Method 1 in S1 and S2 Tables), approximately 14-81 mg of material was carefully removed from the whole. These portions of bones or teeth were submerged in 6% (w/v) sodium hypochlorite (bleach) for 4 min [26] and the bleach poured off. The samples were then twice submerged in DNA-free water, with the water poured off following submersion. Samples were transferred to 1.5 mL tubes, to which aliquots of 500 μL of EDTA (ethylenediaminetetraacetic acid) were added, and gently rocked at room temperature for >48 hours. Samples were extracted in batches of 6-9 with 1-2 accompanying extraction negative controls (i.e., extractions to which no bone or tooth was added). DNA was extracted following the WSU method described by Cui et al. [27].
Under the second extraction method (noted as Extraction Method 1E in in S1 and S2 Tables), approximately 26-59 mg of material was carefully removed from the whole. These portions of bones or teeth were submerged in 6% (w/v) sodium hypochlorite (bleach) for 4 min [26] and the bleach poured off. The samples were then twice submerged in DNA-free water, with the water poured off following submersion. Samples were transferred to 1.5 mL tubes, to which aliquots of 500 μL of EGTA (ethylene glycol tetraacetic acid) were added, and gently rocked at room temperature for >48 hours. EGTA is an alternative decalcifying agent that has been demonstrated to be useful in the study of aDNA [28]. Samples were extracted in batches of seven with an accompanying extraction negative control. DNA was extracted following the WSU method described by Cui et al. [27].
Under the third extraction method (noted as Extraction Method 2 in S1 and S2 Tables), approximately 63-538 mg of material was carefully removed from the whole. The portions of bones or teeth were submerged in 6% (w/v) sodium hypochlorite (bleach) for 4 min [26] and the bleach poured off. The samples were then twice submerged in DNA-free water, with the water poured off following submersion. Samples were transferred to 15 mL tubes, to which aliquots of 2 mL of EDTA were added, and gently rocked at room temperature for >48 hours. Samples were extracted in batches of 6-9 with 1-2 accompanying extraction negative control. DNA was extracted following a modified protocol of Kemp et al. [29] described by Moss et al. [25].
DNA extracts were tested for the presence of co-extracted PCR inhibitors following Kemp et al. [28]. In the case that an extract was deemed to contain sufficient inhibitors to prevent possible amplification of either turkey or dog mtDNA (if present in an extract), it was subject to repeated silica extraction along with its accompanying extraction negative (to monitor inadvertent introduction of contamination). The DNA extracts were again tested for the presence of co-extracted PCR inhibitors. This procedure was carried out until the extracts were deemed to be inhibitor "free". Refer to Kemp et al. [28] for a schematic illustration of the procedure.
Four overlapping amplicons ( 186 base pairs (bps) in length) were sequenced to cover a 460 base pair (bp) portion of turkey mtDNA from nucleotide positions (nps) 15554-16013, a part of the displacement loop, see [30], a comparable stretch of the genome investigated in previous studies [21,22,31] (Table 3). In the case that derived mutations were observed over those previously described for prehistoric Southwestern turkeys [22], re-amplification and resequencing was performed to confirm them, as they could either represent mutations that define novel lineages or be the result of postmortem nucleotide damage [32,33]. Derived mutations that were not replicable were discounted as damage or error.
To maximize comparability between turkeys maintained by humans in the CMV and the NRG prior to the hypothesized migration (i.e., predating 1280 CE), and those maintained in the NRG following the hypothesized migration (i.e., postdating 1280 CE), sequences were truncated to nps 15730-15973. Combined with previously collected (and similarly truncated) comparative data from Speller et al. [22], a total of 189 mtDNA sequences were used to compare turkeys from the three spatial/temporal units of analysis.
Six overlapping canid primer sets ( 202 bps in length), relative to the dog mtDNA reference sequence U96639 [34] and to the coyote reference sequence DQ480509 [35] were used to sequence canid mtDNA and spanned nps 15417-15775 (359 bps) for dogs and nps 15419-15773 (355 bps) for coyotes (Table 4). This stretch of the mitochondrial genome, beginning in the tRNA-Proline gene and ending in the control region [34], is comparable to that investigated in previous studies of prehistoric and extant American dogs (see studies summarized by [11]). One exceptionally common Southwest dog haplotype (derived by 15484G, 15627G, 15639A mutations) became apparent ( Table 2) and replication efforts focused on verifying derived mutations from this general motif. All mutations derived from this haplotype were verified by !2 independent sequence reads. We used failure to confirm derived mutations as an indication that the mutation observed in a single read was most likely a product of postmortem nucleotide damage [32,33]. Some of the same efforts permitted us to verify coyote haplotypes, but given that these animals were likely wild (based on their stable isotopic ratios of carbon and nitrogen), it was not our goal to authenticate the sequences from each of those samples, as they do not contribute to testing our main hypothesis about the movement of domesticate animals. To maximize comparability between sampled from the three spatial/temporal units of analysis, dog sequences were truncated to nps 15459-15691. Table 3. Primers used to amplify turkey mtDNA, their coordinates relative to the turkey (EF153719; [30]) mtDNA reference sequence, amplicon lengths produced, and annealing temperatures. While primers were previously reported by Speller et al. [22], the coordinates of T15533F are corrected here, the target regions have been renamed, and the amplicon sizes are reported.  (Tables 3 and 4), and 68˚C, and 3) a 3 min hold at 68˚C. Successful amplification was confirmed by separating 3-4 μL on 6% polyacrylamide or 2% agarose gels, which were stained with ethidium bromide and visualized under ultraviolet light. Amplicons were sequenced in both directions at either Elim Biopharm (Hayward, CA) or MC Lab (South San Francisco, CA). Sequencher (version 4.8) was used to align the sequences to comparative full mitochondrial genomes (turkey: EF153719 [30], dog: U96639 [34], coyote: DQ480509 [35]).
Fisher's exact tests, conducted using the on-line calculator in VassarStats, were used to compare observed haplotype counts. An alpha level of 0.05 was set as the cut-off for statistical significance of the tests.

Stable isotopic analysis
University of California-Davis. Subsamples of material from 23 canid remains were sent to the University of California-Davis (UCD) for stable carbon ( 13 C/ 12 C) and nitrogen ( 15 N/ 14 N) isotopic analysis. This was a blind test; whereas genetic species identification of these samples was known to the co-author who shipped the remains (B.M.K.), they were unknown to the co-author who processed them for stable isotopic ratios (J.W.E.).
Collagen was extracted by following a modified Longin procedure [36]. Bone samples were cleaned of any adhering soil using a small brush and dental pick, and the outer layer (approximately 0.5mm) removed using a Fordham microdrill. This step removed portions of the bone directly exposed to soil, and hence, most susceptible to diagenetic effects. The sample was then rinsed and sonicated in deionised water (dH 2 O) to remove any loose dust or other adhering material.
To isolate collagen, approximately 1g of cortical bone was demineralized with a solution of 0.5M hydrochloric acid (HCl) set in a refrigerator set to 5˚C. HCl was replaced approximately at 48-hour intervals until the bone was soft and no longer visibly reacted with HCl solution, typically 1-2 weeks. The bone was rinsed three times with dH 2 O and soaked in 0.125M NaOH (sodium hydroxide) for 24 hours to remove humic acids and rinsed again in dH 2 O. To solubilize collagen, pH3 water was added to the vial and the sample placed in an oven set to 80˚C. Every 24 hours, this water was removed, saved, and new pH3 water added until the bone completely solubilized or only residual bone powder remained, typically within 3 days. The vial containing the solubilized collagen was then centrifuged and the fluid pipetted into a clean vial and freeze-dried, isolating the collagen fraction. δ 13 C col and δ 15 N were measured by continuous-flow mass spectrometry (PDZ Europa ANCA-GSL elemental analyzer interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer) at the Stable Isotope Facility, University of California Davis. Carbon isotopes ratios, δ 13 C col , are reported in permil notation (parts per thousand) relative to the PeeDee Belemnite standard (arbitrarily set at 0‰), while N isotope ratios, δ 15 N, are expressed against N 2 in modern atmospheric air (also arbitrarily set to 0‰). The long-term standard deviation for samples in the lab is 0.2‰ for δ 13 C and 0.3‰ for δ 15 N.
Sample quality was evaluated using collagen yield (% weight of mechanically cleaned bone) and atomic C/N ratios. Previous studies suggest that collagen yields greater than 2% produce radiocarbon dates within expected ranges [37]. Assuming similar diagenetic processes for non-radiocarbon ( 13 C and 12 C), we use this figure (! 2% collagen yield) as a cutoff for acceptable samples. Likewise, collagen samples with atomic C/N ratios lower than 2.9 or higher than 3.6 are typically considered degraded and unsuitable for archaeological interpretation [38,39]. We use these same figures to indicate acceptable samples.
Washington State University. Subsamples of material from an additional 12 canid remains were processed for stable carbon ( 13 C/ 12 C) and nitrogen ( 15 N/ 14 N) isotopic analysis at WSU. This was also a blind test; genetic species identification of these samples was unknown to the co-author (L.H.) who processed them for stable isotopic ratios. Sample preparation used recently applied methodologies in current literature [40][41][42]. Approximately 200 mg of cortical bone were removed from the whole. The samples were washed with distilled (DI) water three times and submersed in 1 M HCl for 48 hours at 4˚C for demineralization. After demineralization, the bone fragments were washed three times with DI water, and then submerged in 0.125 M NaOH for 24 hours at 4˚C. The samples were again rinsed 3 times with DI water, and then were immersed in 2 mL of DI water and gelatinized for 24 hours at 90˚C. The soluble "collagen" solution (supernatant) was removed, frozen, and lyophilized using Virtis FM5SL Freezemobile. Approximately 1.25-1.5 mg of the collagen extracts were weighed (Mettler-Toledo #XP56DR), and run in an elemental analyzer (ECS 4010, Costech Analytical, Valencia CA) with a GC column and mass spectrometer to determine δ 13 C and δ 15 N ratios. Normalization of the samples was conducted with acetanilide and corn as quality controls. Standards for comparison of the generated values were C VPDB (Vienna PeeDee Belemnite) for carbon and N AIR (nitrogen content of the atmosphere) for nitrogen.
Once the isotopic readings were generated by isotope-ratio mass spectrometry (IRMS), raw delta values were regressed relative to the standard, to adapt the values to make them comparable to a standard (Eq 1). This was done using the raw delta equation: where R sample is the ratio of an element's normal to heavy isotopes within the sample and R standard is the ratio of an element's normal to heavy isotopes within the standard. The regressed values were then used to create a line from which a linear equation was derived and used to calculate the standard deviations of the standards so the accuracy of the runs could be examined. The normalized delta values of the samples were then compared to the adapted standards using Eq 2: Two tailed t-tests were employed in StatPlus (AnalystSoft Inc.) to compare average δ 13 C and δ 15 N values, setting the alpha level at 0.05 for determining statistical significance of the tests.

Results
One hundred and eighty DNA extractions were performed from 173 turkey samples (S1 Table). An average of 1.03 (SD 1.24) repeat silica extractions were required to sufficiently remove PCR inhibitors from these extractions. Complete or partial mtDNA sequences were recovered from one hundred and forty-four of the one hundred and seventy-three samples (83.2%) (S1 Table). It is notable that turkey mtDNA was neither observed in any extraction negative controls nor in any PCR negative controls.
Complete sequences ranging from nps 15554-16013 were obtained from twenty-three of the one hundred and forty-four of the samples (16%) that yielded analyzable mtDNA (S1 Table). Following established nomenclature [22], 19 of these turkeys belong to the aHap1 lineage, three to the aHap2 lineage, and one to the aHap2b lineage. Thus, 19 of 23 (82.6%) belong to the haplogroup H1 (of which aHap1 is a member), with the remaining four (17.4%) belonging to haplogroup H2 (of which aHap2 and aHap2b are members). Partial sequences from one hundred and four of the one hundred and forty-four samples (72.2%) could be used to confidently categorize individual turkeys as members of either the aHap1 or aHap2 lineage (noting that eventual resolution of these missing data may reveal them to belong to one or more derived lineages). These samples were recorded with question marks as "aHap1?" (n = 102) and "aHap2?" (n = 2) in S1 Table. Combining these complete and partial sequence data reveals that one hundred and twenty-one of the one hundred and twenty-seven samples (95.3%) belong to haplogroup H1 (i.e., aHap1 and aHap1? turkeys) and six of the one hundred and twenty-seven samples (4.7%) are members of haplogroup H2 (i.e., aHap2 and aHap2? turkeys) Partial sequences that do not permit confident assignment to either the aHap1? or aHap2? lineage or to haplogroup H1 or H2 were obtained from the remaining seventeen of the one hundred and forty-four samples (11.8%) that yielded analyzable mtDNA. Turkey mtDNA sequences have been deposited in Genbank (accession numbers MF197043-MF197186).
Our results are combined with those of Speller et al. [22] in Table 1, and Fisher's exact tests of these combined results are given in Table 5. The tests reveal that mitochondrial haplogroup counts in the turkey populations from the CMV and the NRG prior to 1280 CE are significantly different (p = 0.0194) and are highly unlikely to represent a single panmictic population. In addition, there is a statistical difference in haplogroup frequencies between pre-1280 and post-1280 samples in the NRG (p = 0.0212); thus reflecting substantial change in haplogroup frequencies in NRG turkey populations over time. Finally, haplogroup frequencies among turkeys maintained in the pre-1280 CMV and the post-1280 NRG are practically identical (p = 1.0000). These results demonstrate that NRG turkey populations were initially distinct, but became indistinguishable from CMV turkey populations across the hypothesized period of migration from the CMV. One hundred and ninety-one extractions were performed across the 128 canid samples (S2 Table). An average of 0.43 (SD 0.84) repeat silica extractions were required to sufficiently remove PCR inhibitors from these extractions. Mitochondrial DNA was recovered from eighty-one of the one hundred and twenty-eight samples (63.3%). Of these 81 samples, 19 (23.5%) were identified as coyotes ( Table 2, S2 Table). Coyotes were observed at three of the seven pre-1280 CE CMV sites (Albert Porter Pueblo, Shields Pueblo, and McPhee Pueblo). No coyotes were observed in the Northern Rio Grande samples pre-dating 1280 CE and all canid remains sampled from the post-1280 CE NRG (at Abiquiu Pueblo, Arroyo Hondo, and Tsama Pueblo) were identified as coyotes, despite the fact that several of these samples were originally recorded as "dog burials".
Stable isotopic data were generated for all 23 canid remains processed at UC Davis with C/ N ratios of between 1.9 and 3.4 ( Table 6). Two of the 12 samples, 5MT4475-10 and 5MT4475-16, processed at WSU yielded little collagen, produced unacceptable C/N ratios of 4.1 and 5.2, and were eliminated from the study. An average δ 13 C of -15.9 (SD 2.0) and δ 15 N of 10.0 (SD 1.9) was observed for the 15 coyotes samples and an average δ 13 C of -9.1 (SD 3.1) and δ 15 N of 9.1 (SD 1.4) was observed for the 18 dog samples. There is little overlap in the δ 13 C values for the two species (Fig 2), except for dog samples 5MT123-03 (δ 13 C: -18.6) and 5MT5-01 (δ 13 C: -15.9). Both of these dogs lived in the CMV prior to 1280 CE (S2 Table). With regards to 5MT123-03, only a rather short mtDNA sequence was generated (spanning nps 15441-15594), one that matches the reference sequence. However, lacking the 15484G mutations, this individual cannot belong to the most widespread dog haplotype observed in this study (i.e., 15484G, 15627G, and 15639A, see Table 2). While 5MT5-01 exhibits a unique haplotype (i.e., 15625C, 15627G, 15639A), it is one mutational difference (15625C) from a dog lineage observed at Albert Porter Pueblo (sample 5MT123-07). A two-tailed t-test indicates that δ 13 C average values of coyotes and dogs are indeed different (d.f. = 31, t-value = 7.31, p = 3.14x10 -8 ). Based on dietary differences, these genetically identified coyotes were unlikely to have been domestic animals. Thus, only the hypothesized movement of dogs from the CMV to the NRG after AD 1280 was addressed.
Of the 62 dogs that yielded sequences, 56 were complete from nucleotide positions 15459-15691 (Table 2, S2 Table). Dog mtDNA sequences have been deposited in Genbank (accession numbers MF196981-MF197042). A Fisher's exact test of haplotypes observed in the pre-1280 CE CMV and pre-1280 CE NRG provides no evidence of differentiation (p = 0.284) ( Table 7). Since no dogs were identified in the post-1280 CE NRG samples, no comparison was made.

Discussion
Our study confirms the observation of low mtDNA diversity among prehistoric, domestic Southwest US turkeys [22]. Importantly, our increased sampling of turkeys from some archaeological sites over those investigated by Speller et al. [22] (e.g., Shield's Pueblo and Sand Canyon Pueblo (Table 1)), lends support to the notion that there are not frequent haplotypes beyond those initially described. No novel mtDNA lineages were identified in this study.
Turkey mtDNA haplogroup frequencies are consistent with the hypothesis that enough of these domestic birds were moved from the CMV (southwestern Colorado) to the NRG (northcentral New Mexico) during the 13 th century CE to change the mtDNA haplogroup profile of the NRG population. The mtDNA composition of the turkeys maintained in the NRG prior to 1280 CE differs statistically from that of the turkey population post-dating 1280 CE (p = 0.0212), whereas post-1280 CE NRG turkeys are indistinguishable from those maintained in the CMV prior to 1280 CE (p = 1.0000). Notably, however, one matriline exhibited by two turkeys from the post-1280 CE NRG population (i.e., those exhibiting 15782T, 15793T, 15894T, 15953C (Table 1)) was exclusively observed from that unit of analysis and, thus, appears to have not arisen in the CMV prior to 1280 CE. As one specimen, AHT-21, dates to 1310-1345 CE, and the other dates to 1300-1845 CE [22], it is possible that it evolved in the post-1280 CE NRG turkey population. Future aDNA studies of turkeys from the region should be able to address this hypothesis.
It is also important to remember that previous studies have found that turkey remains are, in general, quite rare in NRG sites that pre-date 1200 CE, and that the incidence of turkey within the archaeological record increased substantially in the following century; concurrent Table 6. Results of the stable carbon ( 13 C/ 12 C) and nitrogen ( 15 N/ 14 N) isotopes ratios preserved in the subset of canid remains subjected to analysis. These data are also depicted in Fig 2.  with rapid human population growth [4:292-294, 43]. This archaeological pattern reinforces our finding that the turkey population of the post-migration NRG does not reflect an expansion of the pre-migration turkey population, but rather an influx of birds that accompanied human immigrants. Our results suggest the CMV was a plausible source area for both the turkeys and, by proxy, their human keepers. Prehistoric dogs in the US Southwest exhibit a larger number of matrilines than turkeys. While seven matrilines were detected among the 56 dog samples, the most common lineage  Table 7.

Species
https://doi.org/10.1371/journal.pone.0178882.g002 (i.e., 15484G, 15627G, 15639A) is exhibited by~85.7% of them (forty-eight of fifty-six individuals). Interestingly, approximately 92.6% (one hundred and seventy-five of one hundred and eighty-nine; Table 1) of the turkeys belong to a single matriline. Given the strong genetic evidence for intensive management of turkey populations in the prehistoric US Southwest [18,21,22], it is possible to begin to evaluate a similar scenario for dog populations. Similar patterns of reduced genetic diversity have been observed from dogs recovered from archaeological sites in the US Southeast, and is interpreted as having arisen from deliberate prehistoric breeding practices [11]. Intensive dog breeding practices in the prehistoric US Southwest could also help to explain: 1) the differential diets ascertained for dogs studied here (Fig 2), and 2) the existence of two distinct-looking dog breeds (as ascertained from their mummified remains) recovered from White Dog Cave, Arizona dating to 400 BCE [44]. Future genetic analyses of these specific specimens and other mummified canids, especially in cases where the phenotype of "breed" can be ascertained, would be useful to ascertain whether matriline has relevance to breed, if at all. As only coyotes were identified in the post-1280 CE Northern Rio Grande samples, possible changes in the NRG dog mitochondrial gene pool across this apparent migration boundary could not be addressed. However, unlike observations from turkeys, we found no evidence that pre-1280 CE dogs from the two regions are distinguishable (p = 0.284). This is primarily due to the genetic diversity of the sampled individuals, with a series of relatively rare haplotypes being represented. So it may be that turkeys, with their lower genetic diversity, are a better proxy for the identification of human population movements in any case.
The coyotes detected in our study appear to have been wild based on their stable carbon ( 13 C/ 12 C) isotopic values relative to the majority of dogs. Barring two dogs that have coyotelike δ 13 C values (samples 5MT123-03 and 5MT5-01 (Table 6 and Fig 2)), isotopic variation among coyotes is larger than among dogs. This suggests some degree of dietary plasticity with regards to trophic level among coyotes, a finding that is demonstrated in modern ecological studies of coyote as well [45,46]. Coyotes are able to significantly alter their diet based on food availability and the presence of competitors. By contrast, domesticated dogs appear to have had a more homogenous diet, particularly with respect to accessing protein with C4-derived carbon. The smaller δ 15 N values of most dogs also indicate consumption of lower trophic-level foods, relative to coyotes. Why two dogs stand out isotopically and fall in the coyote cluster in Fig 2 is unclear, but will be the subject of additional analysis. Nevertheless, this demonstrates that genetic species identification of canid remains can improve respective species tabulations in archaeological reports. Genetically identifying species based on DNA would be critical if canid diets are to be used as "proxies" for human diets (which was not our aim here, but see e.g., Guiry [47]).
Overall, however, the findings are consistent with human management or breeding and dietary provisioning of dogs. This specific diet may have put selective breeding pressures on dogs as well. Specifically, dogs that responded well to the δ 13 C-enriched but δ 15 N-depleted diet would have been bred, while those that did not may have been culled. Given recent findings that domesticated dogs have many more copies of genes involved in the digestion of starches and fatty acids relative to wolves [48], we propose that the domestic canine diet included significant amounts of maize. In other words, dogs that responded well to a maize-rich diet were preferentially survived and were bred.

Conclusions
Scholars have debated the role of migration in the collapse of ancestral Pueblo societies of the Four Corners region (including Southwest Colorado) and the coincident emergence of the Rio Grande Pueblos. Previous studies of biological variation using morphometric traits [3, 4:87-124] have found evidence of population movement from the CMV to the NRG during the 13 th century CE. Due to repatriation and native population discontinuity in the four corners region it is not currently possible to confirm these morphometric findings through analyses of contemporary or ancient human DNA. We have therefore considered mtDNA from remains of domestic species as a proxy for human population relationships. Domestic dogs and turkeys were sampled from three temporal/geographic groups for analysis: the pre-migration CMV (pre-1280 CE), the pre-migration NRG (pre-1280 CE), and the post-migration NRG (post-1280 CE).
Our analyses show that many of the canid remains were of coyotes, and isotopic analyses confirm that these coyotes consumed sources of carbon comparatively less enriched in 13 C than those consumed by dogs. Although these results are tangential to our primary research question, they are still informative in that they are consistent with a scenario in which the dogs consumed a diet rich in maize, perhaps supplemented by some turkeys that also ate maize, whereas the coyotes likely subsisted on a wide range of wild resources rich in C3-derived carbon. As a result of no dogs being identified from post-migration NRG sites, we lack the necessary samples to test for discontinuity between pre-migration and post-migration NRG dogs. However, we were not able to detect statistically-significant genetic differences between premigration CMV and NRG dogs.
In contrast, our analyses of turkey mtDNA, combined with results of previous studies, provide evidence that the turkey population of the post-migration NRG exhibits genetic continuity with the pre-migration CMV turkey population, and both exhibit genetic discontinuity with the pre-migration NRG turkey population. These results are consistent with an influx of turkeys to the NRG during the 13 th century, and the CMV was not excluded as a possible source for these birds. They also suggest the post-1280 CE NRG turkey population could not have derived wholly from the pre-1280 population from this region. Our results are consistent with the hypothesis that migrating human populations from the CMV brought domestic turkeys with them to the NRG, leading to a change in the mtDNA composition of the post-1280 CE NRG turkey population. This study thus provides the first direct genetic evidence, by way of a domestic animal proxy, in support of the hypothesis that Tewa communities occupying the NRG region today derive in large measure from the migration of 13 th century CMV populations. If true, this would support a scenario in which migration was a primary solution of CMV populations in the face of significant social and environmental problems [2,4,[49][50][51], as it often is for human groups under extreme duress today.
Supporting information S1 Table. Turkey samples studied with their provenience, dates, and other pertinent archaeological information, and the full results of their genetic analyses. Under the amplicon columns, an "X" in a cell indicates successful amplification and sequencing of the amplicon, whereas "O" indicates results from the amplicon were not obtained. Partial sequences that could be used to confidently categorize individual turkeys as members of either the aHap1 or aHap2 lineage, but for which eventual resolution of these missing data may reveal them to belong to one or more derived lineages, are indicated with question marks as "aHap1?" or "aHap2?". Sequence reads and mutational positions are relative to the turkey mtDNA reference sequence (EF153719; [30]). (XLSX) S2 Table. Canid samples studied with their provenience, dates, and other pertinent archaeological information, and the full results of their genetic analyses. Under the amplicon columns, an "X" in a cell indicates successful amplification and sequencing of the amplicon, whereas "O" indicates results from the amplicon were not obtained. Cells colored black indicate that no attempt to amplify the amplicon was made. Dog sequence reads and mutational positions are relative to the dog mtDNA reference sequence (U96639; [34]). Coyote sequence reads and mutational positions are relative to the coyote mtDNA reference sequence (DQ480509; [35]). (XLSX)