We submitted a research proposal to the American Museum of Natural History (AMNH) requesting samples from Chaco Canyon burials classified as culturally unidentifiable following NAGPRA (Native American Graves Protection and Repatriation Act) criteria. The AMNH review committee approved our proposal in accordance with all NAGPRA legal guidelines governing research.

AMS 14C and stable isotope analyses

Methods used for AMS 14C dating crania 13 and 14 are published elsewhere11. In the present study, we directly dated nine crania recovered from above the split plank floor in room 33. Previously published radiocarbon dates for interments above the floor came from long bones15. This skeletal material was commingled, and hence we directly AMS 14C dated each of the crania to verify the date of interment.

Dentine was extracted from tooth roots for 14C and stable isotope analyses and was extracted and purified using the modified Longin method with ultrafiltration30. Tooth samples were initially cleaned of adhering sediment and the exposed surfaces were removed with an X-acto blade. Samples (200–400 mg) were demineralized for 24–36 h in 0.5 N HCl at 5 °C followed by a brief (<1 h) alkali bath in 0.1 N NaOH at room temperature to remove humates. The residue was rinsed to neutrality in multiple changes of Nanopure H 2 O, and then gelatinized for 12 h at 60 °C in 0.01 N HCl. Resulting gelatin was lyophilized and weighed to determine percent yield as a first evaluation of the degree of bone collagen preservation. Rehydrated gelatin solution was pipetted into precleaned Centriprep 30 ultrafilters (retaining >30 kDa molecular weight gelatin) and centrifuged 3 times for 20 min, diluted with Nanopure H 2 O and centrifuged 3 more times for 20 min to desalt the solution. Carbon and nitrogen concentrations and stable isotope ratios were measured at the Yale Analytical and Stable Isotope Center with a Costech elemental analyzer (ECS 4010) and Thermo DeltaPlus analyzer. Sample quality was evaluated by % crude gelatin yield, %C, %N and C/N ratios before AMS 14C dating. C/N ratios for all samples fell between 3.11 and 3.16, indicating good collagen preservation31. Shell samples were etched in dilute HCL and then hydrolysed in vacutainers with 85% phosphoric acid for 1 h to generate CO 2 .

Collagen samples (∼2.1 mg) were combusted for 3 h at 900 °C in vacuum-sealed quartz tubes with CuO and Ag wires. Sample CO 2 from collagen and marine shell was sent to KCCAMS where it was reduced to graphite at 550 °C using H 2 and a Fe catalyst, with reaction water drawn off with Mg(ClO 4 ) 2 (ref. 32). Graphite samples were pressed into targets in Al cathodes and loaded on the target wheel for AMS analysis. The 14C ages were corrected for mass-dependent fractionation with measured δ13C values,33 and compared with samples of Pleistocene whale bone or calcite (backgrounds, >48 14C kyr BP), late Holocene bison bone (∼1,850 14C BP), late AD 1800s cow bone and OX-2 oxalic acid standards for calibration. Results are presented in Supplementary Table 1.

Calibration and modelling

AMS 14C dates below the split plank floor were calibrated with OxCal version 4.2.3 (ref. 34) within a stratigraphic model using the IntCal13 northern hemisphere curve35. Using Sequence and Phase models in OxCal, a rough ordering of depositional events in room 33 below the split plank floor was used to help constrain the calibrated ages of for crania 14 and 13. The initial construction or preparation of the room for use as a crypt is modelled with a Boundary followed by a Phase (‘room 33 lower’) including crania 13 and 14, below the plank floor. Two directly AMS 14C dated abalone (Haliotis spp.) shells were also included in this model. The emplacement of the plank floor is modelled as a Boundary followed by a Phase (‘room 33 upper’). Crania about the floor were not modelled because of disturbances during subsequent interment events and they cannot be confidently placed in stratigraphic order with respect to each other. Calibrated results are presented in Supplementary Table 2 and Supplementary Fig. 1. The OxCal code used for the overall model is provided in Supplementary Note 2.

The abalone shells were modelled using the Marine13 curve35 and a reservoir (ΔR) correction of 234±23 14C year assuming collection somewhere on the southern California coast in the vicinity of San Diego where this species is concentrated36.

Radiocarbon simulations

The apparent episodic structure in the calibrated ages for the room 33 crania (n=11) could indicate several periods of interment or potentially continuous interment over the course of three centuries. When plotted against the IntCal13 curve35 (Supplementary Fig. 2), many of the dates fall on reversals and plateau in the calibration curve between 900 to 1150 CE, similar to observed patterns for directly dated macaw remains from Pueblo Bonito as reported elsewhere21. Conventional ages that occur within plateaus in the calibration curve cause the artificial stacking or overrepresentation of certain calibrated dates within certain intervals, whereas conventional dates falling on steeper parts of the curve tend to be underrepresented.

The potential biases introduced by the calibration curve were explored by simulating radiocarbon ages from 750 to 1250 CE in OxCal v4.2.3 (refs 34, 37). The R_Simulate command specifies a calibrated age and a measurement error, translates the age through the curve to generate a conventional age BP and then calibrates the conventional age. A model was made to simulate 10 dates with a precision of ±20 14C years every 25 calendar years from 750 to 1250 CE for a total of 210 simulated dates. The conventional ages were extracted and frequencies plotted in 20 14C-year increments (Supplementary Fig. 3).

Comparing the frequency of simulated 14C ages with the expected (40 per century), it would be anticipated that even with more or less continuous (for example, generational) interment of individuals in room 33 that fluctuations in the calibration curve produce an apparent episodic structure on the conventional and calibrated ages through the modelled period. Given that the sample size (n=11) is fairly small compared with the overall period of interment (∼330 cal years), it is difficult to argue that certain periods are over- or underrepresented, as was possible with the procurement of macaws over a shorter interval21. For room 33, we do observe that a few gaps in the 14C ages do correspond to periods where the calibration curve is relatively steep, and therefore are likely artefacts of the curve rather than pauses in interment. Overall, the simulations indicate the observed 14C ages of room 33 are consistent with continuous rather than episodic use of the crypt from the ninth through the twelfth century CE.

Ancient mtDNA capture and sequencing

Skeletal remains from nine individuals from room 33 (Supplementary Data 1) were successfully analysed at the ancient DNA facilities at Pennsylvania State University and Harvard Medical School. Eight of the nine room 33 samples for which archaeogenomic data are analysed in this study (excepting cranium 6) were processed initially in the Penn State ancient DNA facility, with the mtDNA genome capture and sequencing performed as described in this section. The mtDNA genome sequences for three of these same individuals (crania 1, 3 and 7) were also determined at the Harvard Medical School ancient DNA facility, as part of the screening process for the nuclear genome SNP genotyping analyses that were performed there (see below). The mtDNA genome of one room 33 individual, cranium 6, was sequenced only at Harvard Medical School.

At Penn State, DNA extractions and pre-PCR library preparation steps were performed in a dedicated sterile laboratory with high-efficiency particulate arrestance-filtered air and positive air flow. Strict procedures were followed to prevent contamination and provide a sterile laboratory environment. All surfaces and workstations were sterilized using a concentrated bleach solution and 70% ethanol solution, before and after each use. Equipment and consumables were irradiated under ultraviolet light for 60 min before samples were introduced. Dedicated reagents were only opened in a sterilized ultraviolet hood located within the clean laboratory. Negative control reactions were included alongside all DNA extractions and library preparations. DNA extractions and negative control reactions were prepared and processed in tandem.

Individual ancient samples were initially pretreated in 1% bleach solution in a separate nonmolecular laboratory before entering the dedicated sterile laboratory. To further reduce surface contamination, the surface layer of each sample was removed using a Dremel tool fitted with a disposable rotating disc that had been treated in 1% bleach solution, followed by multiple washes in molecular grade water to remove traces of bleach. After surface treatment, the samples were ground into a fine powder using a mikro-dismembrator (Sartorius) ball mill with tungsten carbide balls. For each sample, we extracted DNA from 150 to 300 mg of bone powder exactly following a silica absorption/spin column method38 and eluted DNA twice in 50 μl TE buffer with 0.05% Tween-20.

Double-stranded libraries with barcoded adapters were constructed with 50 μl of each ancient DNA extract based on the approach by Meyer and Kircher39. Before the amplification step of this protocol, the libraries were split into two 10 μl aliquots, with different barcode indices added to each aliquot (Supplementary Data 1). Reactions were prepared in the ancient DNA facility, sealed and then amplified in a separate lab. All libraries were amplified in a 50 μl reaction consisting of PCR buffer, 2 mM MgSO 4 , 200 μM each dNTP, 200 nM primer IS4 (5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTT-3′), 200 nM barcoded p7 primer (5′-CAAGCAGAAGACGGCATACGAGATxxxxxxGTGACTGGAGTTCAGACGTGT-3′, where ‘xxxxxx’ represents the sample-specific barcode) and 2.5 U Platinum Taq High Fidelity DNA Polymerase (Thermo Scientific). Libraries were amplified as follows: 5 min at 94 °C; 24 cycles of 20 s at 94 °C, 15 s at 60 °C and 20 s at 68 °C; final extension 5 min at 60 °C. Indexed PCR products were purified using SPRI beads and eluted in 15 μl TET buffer.

DNA libraries were enriched for mtDNA fragments using an in-solution biotinylated RNA bait hybridization method40. We designed a probe set of 100-mer baits with 10 bp tiling (that is, a new bait starting every 10 bp) complementary to the sequences of five human mitochondrial reference sequences representing the different mtDNA haplogroups observed among Native Americans (GenBank accessions EU095194.1 (A2), EU095219.1 (B2), EU095222.1 (C1), EU095232.1 (D1) and EU095242.1 (haplotype X2a)). The baits were synthesized by MyCroarray, Inc. (probe design: 140429).

The two separately indexed libraries from each sample were combined before hybridization capture. DNA captures were performed using the MYcroarray MYbaits system based on the manufacturer’s protocol (Version 2.3) with two modifications: (1) up to 1 μg of the amplified DNA library was used for each sample, and (2) twice as many wash steps were completed following hybridization to remove more nontargeted DNA fragments. Post-capture libraries were amplified using IS5 and IS6 primers39 and Kapa Hifi DNA Polymerase, otherwise following the MYbaits protocol. After amplification, libraries were purified with SPRI beads and eluted in 15 μl TET buffer.

Post-capture barcoded libraries were pooled and sequenced at the Penn State Huck Institutes Genomics Core Facility on an Illumina HiSeq 2500 platform using 76-bp paired-end reads in Rapid Run mode. Forward and reverse sequence reads were trimmed of adapter sequences and overlapping reads were merged41, enforcing a minimum 11-nt overlap and base quality score of 20. Merged sequences were mapped to a human reference mitogenome (GenBank accession EU256375.1) using the BWA-backtrack algorithm in the Burrows–Wheeler Aligner, version 0.7.5 (ref. 42). Some samples were sequenced on multiple lanes to increase sequence coverage. All BAM files from these samples were merged to a single per-sample alignment using the SAMtools merge command. Potential PCR-duplicated reads were then removed using the rmdup command from SAMtools-0.1.19 (ref. 43. Reads <20 bp were removed to limit error from nonspecific mapping of exogenous DNA fragments.

Ancient DNA damage analysis

Characteristic damage patterns in DNA degradation over time can be used to assess the authenticity of ancient DNA samples24,44. Ancient DNA damage is detected by observing increases in specific nucleotide misincorporation patterns diagnostic of mismatches in the single-stranded overhanging ends of sequenced DNA fragments. We analysed each alignment using mapDamage 2.0 (ref. 45) (Supplementary Fig. 4), and extracted mean lambda values from the Bayesian estimation of damage characteristics, describing the per-base probability of terminating an overhang creating a geometric distribution of overhang lengths. We used the cumulative geometric distribution within the R statistical environment46 to calculate the overhang length encompassing 95% of inferred overhangs per value of lambda (6–14 nt; Supplementary Data 1). We hard-masked all 5′ T and 3′ A residues within this interval and called a consensus sequence using the SAMtools mpileup command and a perl script, enforcing 2 × nonredundant coverage and 80% site identity. Mitochondrial genome haplogroups (Supplementary Data 1) were identified using the program Mitomaster47 that uses the Mitomap information system and Haplogrep48. Individuals with the mtDNA B haplogroup have a characteristic 9-bp deletion compared with the reference genome sequence at positions 8,271–8,279 (9-bp sequence CCCCCTCTA)49,50,51,52. In this region, the consensus sequences for haplogroup B individuals were erroneous because of misalignment and were manually corrected based on visualization of the sequence reads.

To verify support for the consensus haplotype among only reads with deaminated cytosine residues, we used PMDtools53 to generate BAM files from fragments having a minimum post-mortem degradation (PMD) score of 3, retaining ∼15–43% of reads per alignment. We re-ran mapDamage to generate new lambda values, masked the PMD BAM files and called consensus sequences as above. We called consensus sequences as above with 80% site identity, but with an increased 4 × minimum nonredundant coverage to account for the preponderance of damaged sites in the PMD alignments. PMD consensus sequences matched the originals at 100% of represented sites.

Comparison with modern mtDNA genomes from the Americas

We curated a reference database of 171 modern complete mitochondrial DNA sequences from the Americas54,55,56,57,58,59,60,61,62,63, restricted to genomes with ethnicity or region of origin information (Supplementary Data 2). The ancient mtDNA genome consensus sequences from the Chaco region generated in this study were aligned to the modern reference sequences using MAFFT64 and manually verified using Geneious v8.1.8.

Genetic sex estimation

The genetic sex of each ancient individual was estimated using the R y approach25 that is based on the ratio of chrY to chrX+Y mapped reads. For this analysis we used two approaches. First, for the results shown in Fig. 2 we used the sequence reads obtained following the mitochondrial DNA capture enrichment that did not map to the human mtDNA genome (that is, ‘bleed-through’ reads, or those reads that were not removed from the library during the washing steps of the DNA capture protocol (DNA capture is not 100% efficient at reducing library representation). The sequence reads from each sample were aligned to the human reference sequence hg18 (NCBI36/hg18 human genome assembly) using the BWA-backtrack algorithm in the Burrows-Wheeler Aligner, version 0.7.5 (ref. 42) and processed as above to remove potential PCR duplicate reads. Sequence alignments were restricted to those with a mapping quality score of at least 30. The genetic sexes of six samples could be estimated confidently, whereas for two other individuals our results were consistent with an assignment of male but the hypothesis that these individuals were female could not be rejected based on the 95% confidence interval (Supplementary Data 1).

Second, we applied the same the same R y genetic sex estimation approach25 to data from the nuclear genome SNP genotyping experiments described below for the six individuals with the highest levels of endogenous DNA preservation. This analysis was limited to sequence reads that mapped to the targeted chromosome X and Y SNPs with a mapping quality score of at least 37, where clusters of duplicate reads (as identified by orientation, and start and end position) were represented by the read with the highest sequence quality. The sex estimates from this analysis were 100% concordant with those based on the mtDNA capture bleed-through data (Supplementary Data 1).

We also compared the genetic sex estimates for the eight room 33 individuals with two previous osteological sex determinations18,65, showing 100% concordance with the more recent assessment18 and resolving the discrepancy for one individual between the two osteological studies (Supplementary Data 1).

Nuclear genome SNP genotyping

In a dedicated ancient DNA clean room at Harvard Medical School, we generated nine ancient DNA libraries, either from powder drilled directly from teeth end extracted using a protocol optimized for the recovery of ultra-short DNA fragments66 or using extracts made in the ancient DNA clean room at Penn State. We converted the extracts into sequencing libraries, most of which were treated with uracil-DNA-glycosylase67.

We first used in-solution capture method to enrich the DNA library for human mitochondrial DNA fragments5,23 to screen for samples with sufficient quality for nuclear genome SNP genotyping. We sequenced the enriched DNA libraries on an Illumina NextSeq500 instrument using 2 × 75 base pair reads. We used Seqprep68 to identify paired sequences overlapping by at least 15 base pairs, and then mapped these reads to the mitochondrial DNA genome RSRS56 using BWA-0.6 (ref. 42).

We took forward to genome-wide analysis samples with appreciable amounts of DNA and no evidence of heterogeneity of mitochondrial sequences based on contamMix-1.0.9 (ref. 69). One accepted strategy for ancient DNA authentication in genome-wide ancient DNA analysis is to require at least 10% damage in the first nucleotide for a non-UDG-treated library and at least 3% in the first nucleotide for a partially UDG-treated library, and all libraries we took forward to genome-wide analysis met this requirement67.

We enriched the libraries taken forward to genome-wide analysis for 1,237,207 SNP targets (SNP panels 1 and 2 from Fu et al.23 that include all SNPs on the Affymetrix Human Origins, Illumina 610-Quad and Affymetrix 50k arrays). We sequenced the enriched products on a NextSeq500 instrument using 2 × 75 base pairs reads, and processed the sequences as we did for the mitochondrial DNA analysis except that we mapped instead to the human genome reference sequence hg19. We restricted our analysis to six samples for six distinct libraries passing our quality control. These samples had 14,948–144,227 SNPs covered by at least one sequence read.

Relatedness coefficient estimates

We extracted 856,473 SNPs on chromosomes 1–22 that are not in CpG dinucleotide context in order to avoid the potential influence of differences in postmortem damage rates among individuals. For each pair of individuals, we computed the average mismatch rate across all autosomal SNPs covered by at least one sequence read for both of the two individuals being compared (when >1 sequence read was present for one individual at a given site, a random read was sampled for the analysis) and computed s.e. using a weighted block jackknife. There were too few X-chromosome SNPs with sequence data for high-confidence relatedness results.

We observed that the greatest mismatch rates x between pairs were ∼21%, but several were significantly lower (Z>3), providing putative evidence of relatedness. To estimate the relatedness coefficient r, we need to account for the fact that we are downsampling the data to a single sequence at each position. Thus, if we were comparing two identical individuals, the mismatch rate would still be only half of that observed between unrelated individuals. We thus correct our estimate of r using the expected mismatch rate because of random sequence sampling b. We choose b=0.21/2=0.105 based on the approximate maximum mismatch rates observed. Our estimator is thus

We also computed a 95% confidence interval using block jackknife standard errors.

Data availability

Sequence reads from the mtDNA capture experiments and all nuclear genome SNP genotyping BAM files have been deposited in the NCBI Sequence Read Archive, Accession no. SRP094965 (SRA BioProject no. PRJNA353635). All BAM and FASTA files of both the original and damaged read version consensus mtDNA sequence for each individual have been deposited in the Dryad Digital repository at http://dx.doi.org/10.5061/dryad.3344d. The multisequence alignment has been deposited in the Dryad Digital repository at http://dx.doi.org/10.5061/dryad.3344d.