Archaeological sites background

The four archaeological sites analyzed in this study date to the period A.D. 700–1450 (Table 1). The two Greater Nicoya sites in northwest Costa Rica, Jícaro, and La Cascabel (Fig. 2a), correspond to the Sapoá Period (A.D. 800–1250)104. Excavations at Jícaro were conducted over three seasons between 2005 and 2008105 (Fig. 2b) while excavations in La Cascabel were conducted in 2007106.

Samples analyzed from the northern Mexican Casas Grandes site of Paquimé date to the Medio Period (AD 1200–1450), and one sample from the Convento site dates to the Viejo Period (A.D. 700–1200). Paquimé and Convento are located in the cultural borderlands between the American Southwest (Greater Southwest and Aridoamerica) and Mesoamerica and were systematically excavated from 1958 to 196132 (Fig. 2c,d).

Genetic investigation of skeletal material from the sites of Jícaro, La Cascabel, Paquimé, and Convento, additionally present an opportunity to observe the impact of climate on ancient DNA preservation in northern Mexico and Central America. Jícaro and La Cascabel are found in the tropical dry forest of northwest Costa Rica, 10.6° latitude from the Equator, with an elevation of 70 meters above sea level (masl) and annual temperatures 26–33 °C. The average annual precipitation for this area is >1500 mm107. In contrast, Paquimé and Convento are found in the dry climate of Chihuahua, at latitude of 30.3° and elevation of 1457 masl, with annual temperatures between −10 °C and 42 °C. The average annual precipitation for this area is approximately 326 mm108. Given the differences between these two sites, DNA preservation is expected to be poorer at the Costa Rican sites, where higher humidity, lower altitude, and equatorial latitude are likely to accelerate DNA decay.

Archaeological samples

A total of fourteen individuals were selected for this study (Table 1), with one tooth analyzed per individual, except for Jícaro sample ID_7 for which two teeth were tested (Table 2). Thus, fifteen teeth were processed in total. Nine individuals from Greater Nicoya were sampled at the National Museum of Costa Rica. The sites of Jícaro and La Cascabel were selected because of their notable skeletal preservation105,106, which included the preservation of aging indicators on the pelvic auricular surface and dental occlusal surfaces109. The five Casas Grandes individuals were sampled from Museo de las Culturas del Norte, INAH Chihuahua, as part of a larger project on Casas Grandes bioarchaeological research at the University of Calgary. A map with the approximate locations from where the samples were excavated (Fig. 2) was produced using Esri’s ArcMapTM 10. The site maps were redrawn and overlaid on the satellite imaginary basemap provided by Esri.

Permission to sample the Casas Grandes material was provided by the INAH-Chihuahua and Museo de las Culturas del Norte (Paquimé), and ethical review for genetic analysis was performed and approved by the University of Calgary Conjoint Faculties Research EthicsBoard (REB15-2680). For the Costa Rican material, permission for sampling and analysis was granted by the National Museum and Minister of Culture of Costa Rica.

DNA extraction

Archaeological dental samples were processed in dedicated aDNA facilities at the University of Calgary’s Ancient DNA Laboratory and at the University of Oklahoma’s Laboratories of Molecular Anthropology and Microbiome Research (LMAMR). At both facilities, all sample decontamination steps and DNA extraction methods followed established procedures to prevent contamination and specific details are provided below. A non-template extraction control (negative control) was processed alongside the experimental samples during all analytical steps to monitor for the presence of contamination.

Ancient DNA Laboratory, University of Calgary

Prior to beginning this study, 26 samples from Jícaro and La Cascabel and one sample from Paquimé were screened for DNA preservation. Only six of these samples (Sample IDs_2, 6, 7, 8, 9, and 10 in Table 2) were determined to be suitable for further study. DNA extraction was performed as follows. Samples were submerged in 6% sodium hypochlorite solution for 10 minutes, rinsed twice with ultrapure water, and irradiated with 254-nm ultraviolet light for 30 min per side to reduce surface contamination. Afterward, whole samples were pulverized, without removing the crown, and a total of 1–2 g of tooth powder was transferred to a 15 ml tube with 5 ml of extraction solution (0.5 M EDTA pH8, 0.12% SDS, 0.5 mg/ml proteinase K) and incubated at 50 °C overnight with agitation. Following this, DNA was purified and concentrated using a MinElute silica spin column protocol110 and eluted twice into 50 μl aliquots of TET buffer.

LMAMR, University of Oklahoma

DNA was extracted from twelve samples in a dedicated aDNA laboratory at the University of Oklahoma, including four samples that were previously extracted at the University of Calgary (Table 2). Prior to extraction, dental calculus was removed from all dental samples, and the tooth surface was gently cleaned with a 2% sodium hypochlorite solution followed by ultrapure water. The tooth root surface was removed by mechanical abrasion using a Dremel rotary tool, separated from the crown, UV-irradiated on each side for 1 min, and then pulverized to a coarse powder. Approximately 100 mg of tooth powder was pre-digested in a 0.5 M EDTA solution for 15 min and decanted to further remove contaminants. The cleaned tooth powder was then incubated in a solution of 1 ml of 0.5 M EDTA for 24 h at room temperature under agitation, followed by the addition of 100 µl of Proteinase K (Qiagen) and further incubation until digestion and decalcification were complete. DNA was isolated following a silica adsorption protocol111 using a Qiagen MinElute PCR Purification kit and twice eluted into 30 µl of EB buffer (10 mM Tris-Cl, pH 8.5) for a combined elution volume of 60 µl.

PCR amplification of the mitochondrial hypervariable region

At the University of Calgary, amplification of approximately 600 bp of the mitochondrial hypervariable regions I and II (HVR-I and HVR-II) was attempted for 6 samples (Table 2) extracted there. The region was amplified by targeted PCR in four overlapping segments of approximately 140 bp each, with the following primer pairs: F15989 - R16158, F16112 - R16251, F16190 - R16322, F16268 - R16410. Primer design followed Gabriel and colleagues (2001) except for R16251 [5′-GGA GTT GCA GTT GAT GT-3′]. The PCR reaction final volume was 30 µl and contained 2.5 mM MgCl 2 , 0.2 mM dNTP, 1.0 mg/ml Bovine Serum Albumin (BSA), 0.3 µM forward and reverse primer, and 2.5 U AmpliTaq GoldTM enzyme. A total of 3–5 µl of DNA extract was added to each PCR reaction, which was amplified according to the following conditions: initial denaturing at 95 °C for 12 min, followed by 40 cycles at 95 °C for 30 sec, 50 °C for 30 sec, 72 °C for 30 sec, followed by a final extension of 72 °C. Electrophoresis on 2% agarose gels was used to visualize positive amplifications of targeted fragments. PCR products were sequenced using forward and/or reverse primers at Eurofins Genomics (Louisville, KY). The resulting chromatograms were edited using ChromasPro 2.6.2 software to remove primer sequences and low quality bases and aligned to the revised Cambridge reference sequence (rCRS NC_012920; Andrews et al.112 using BioEdit 7.2.5 software to determine variant positions.

HTS library construction

Illumina sequencing libraries were prepared for a total of fifteen samples (Table 2). All initial library preparation steps were performed in the University of Oklahoma’s LMAMR dedicated ancient DNA facility following previously published protocols113, with minor modifications. In brief, up to 30 μl of sample extract (to a maximum of 100 ng input DNA) (Supplementary Table 2) was constructed into Illumina libraries using the NEBNext DNA Library Prep Master set for 454 (E6070; New England Biolabs). Library preparation was modified from the manufacturers’ instructions, decreasing the total volume to 50 µl, and also replacing SPRI bead purification with silica column purification (Qiagen MinElute PCR Purification kit). Blunt end adapters (IS1/IS3 and IS3/IS2) were prepared and used for ligation at a final concentration of 0.6 µM in a final volume of 50 µl. The Bst polymerase fill-in reaction was inactivated by heating at 80 °C for 20 min, followed by freezing overnight.

For each archaeological sample, the library was amplified in triplicate using 0.75 µl of each indexed primer at 10 µM, 12.5 µl of KAPA HiFi Uracil+ (Kapa Biosystems), 1 µl of BSA at 2.5 mg/ml, 4 µl of library template, and 6 µl of water for a final volume of 25 µl. PCR conditions were as follows: initial denaturation at 95 °C for 5 min, followed by 10–16 cycles of denaturation at 98 °C for 20 s, annealing at 60°C for 15 s, and elongation at 72°C for 30 s, followed by a final elongation at 72°C for 1 min. The number of cycles necessary for the library amplification was determined by Real-Time PCR, selecting the midpoint of the exponential phase of amplification for each sample. All negative extraction controls were amplified for 16 cycles, equal to the greatest cycle number for an archaeological sample. The success of library preparation was determined by visualization on a 2% agarose gel and the triplicate reactions were pooled. The reactions were purified using the AMPure substitute protocol114 and eluted into 10 ul RNase free water.

Mitogenome in-solution capture and sequencing

Mitochondrial in-solution capture was performed on all shotgun libraries following the protocol described by MYbaits V3.0, with minor modifications. Briefly, each pool of libraries was incubated with RNA probes and buffers for 48 h in a touchdown enrichment that started at 65 °C for 16 h, then 60 °C for the next 16 h, and finally 55 °C for the last 16 h. Libraries and beads were incubated for 30 min at 55 °C and Wash Buffer 2.2 was heated to 55 °C. The subsequent steps of beads washing and releasing from the RNA bait were conducted as described in the MYbaits V3.0 protocol. Real-Time PCR was used to determine the cycles necessary for post-capture amplification, and post-capture PCR was performed in triplicate reactions of 25 µl each containing the following: 1X KAPA HiFi Uracil+, 0.1 mg/ml BSA, 0.3 µM of each Primer-Illumina ext, 4 µl of template, and water. PCR conditions were as follows: initial denaturation at 98 °C for 2 min, followed by 14–25 cycles of denaturation at 98 °C for 20 s, annealing at 60°C for 30 s, and elongation at 72 °C for 30 s, followed by a final elongation at 72°C for 5 min.

The amplified post-capture libraries were quantified using a Bioanalyzer 2100 DNA 1000 assay (Agilent) and pooled into equimolar ratios. Size selection to remove primer dimers was performed using a PippinPrep (Sage Science), and the combined pool was sequenced using Illumina HiSeq v2 2 × 100 bp chemistry at the Yale Center for Genome Analysis (YCGA).

HTS data analysis

Sequence reads were processed using a bioinformatics pipeline optimized for ancient DNA analysis. Post-capture paired reads were merged and adapters were trimmed with the program leeHom using the ancient DNA flag115. The reads were mapped to the rCRS112 using a modified version Burrows Wheeler Aligner (-n 0.01, -o 2, -l 16500)116 (https://github.com/mpieva/network-aware-bwa). Reads were filtered to those that were mapped and merged or properly paired (https://github.com/grenaud/libbam) and duplicate molecules collapsed on the basis of identical 5′ and 3′ mapping positions (https://bitbucket.org/ustenzel/biohazard). Mapped reads were further filtered to those with a minimum length of 35 bp and minimum quality of 30117. Schmutzi75 was utilized to estimate endogenous mitochondrial consensus sequences and contamination (−qual 10) using the comparative Eurasian database. Automated indel calls were manually confirmed directly from the read pileups and corrected to avoid alignment errors. The assignment of mitochondrial haplogroups was performed using Haplogrep 2.0118,119,120. Ancient DNA damage patterns were assessed using MapDamage v.274.

Data availability

Genetic data are available in the NCBI Short Read Archive (SRA) under the BioProject accession PRJNA360145 (SRA accessions SAMN06204047 to SAMN06204061).