Ancient DNA analysis

We generated powder from the skeletal remains of all individuals excavated from sites throughout the Caribbean (seeSupplementary Information section 2 for archaeological site information andFigures S1-S11 for maps showing the location of the islands and/or sites studied). Powder was produced from a cochlea^38,39, tooth, phalanx, or ossicle⁴⁰ from each individual in a clean room facility at Harvard Medical School (Boston, USA), University College Dublin (Dublin, Ireland), or the University of Vienna (Vienna, Austria); seeSupplementary Data 2 for the skeletal element used for each individual and location of powder preparation.

We extracted DNA in dedicated ancient DNA laboratories at Harvard Medical School or the University of Vienna following published protocols^41–43. From the extracts, we prepared dual-barcoded double-stranded⁴⁴ or dual-indexed single-stranded libraries⁴⁵, both treated with uracil-DNA glycosylase (UDG) to reduce the rate of characteristic ancient DNA damage⁴⁶. Double-stranded libraries were treated in a modified partial UDG preparation⁴⁴ (‘half’), leaving a reduced damage signal at both ends (5’ C-to-T, 3’ G-to-A). Single-stranded libraries were treated withE. coli UDG (USER from NEB) that inefficiently cuts the 5’ Uracil and does not cut the 3’ Uracil. For a subset of individuals, we increased coverage by preparing multiple libraries; seeSupplementary Data 2 for the number of libraries analyzed for each individual.

To generate SNP capture data, we used in-solution target hybridization to enrich for sequences that overlap the mitochondrial genome and ~1.24 million genome-wide SNPs^47–50 (“1240k”), either in two separate enrichments or simultaneously (Supplementary Data 2). We then added two 7-base-pair indexing barcodes to the adapters of each double-stranded library (single-stranded libraries are already indexed from the library preparation) and sequenced libraries using either an Illumina NextSeq500 instrument with 2×76 cycles or an Illumina HiSeqX10 instrument with 2×101 cycles and reading the indices with 2×7 cycles (double-stranded libraries) or 2×8 cycles (single-stranded libraries).

Prior to alignment, we merged paired-end sequences, retaining reads that exhibited no more than one mismatch between the forward and reverse base if base quality was ≥20, or 3 mismatches if base quality was <20. A custom toolkit was used for merging and trimming adapters and barcodes (available athttps://github.com/DReichLab/ADNA-Tools). Merged sequences were mapped to the reconstructed human mtDNA consensus sequence (RSRS)⁵¹ and the human reference genome version hg19 using the samse command in BWA v.0.7.15-r1140⁵² with the parameters -n 0.01, -o 2, and -l 16500. Duplicate molecules (those exhibiting the same mapped start and end position and same stand orientation) were removed after alignment using the Broad Institute’s Picard MarkDuplicates tool (available athttp://broadinstitute.github.io/picard/). We trimmed two terminal bases from UDG-half libraries to reduce damage-induced errors.

We evaluated the authenticity of the isolated DNA by retaining individuals with a minimum of 3% of cytosine-to-thymine substitutions at the end of the sequenced fragments⁴⁴ for double stranded libraries and 10% for single-stranded libraries, point estimates of mitochondrial DNA (mtDNA) contamination below 5% using contamMix v.1.0–12⁴⁷, and point estimates of X chromosome contamination (in males) below 3%⁵³; we also used contamLD⁵⁴ to confirm low contamination rates (<~6%) (Supplementary Data 2). Eight single-stranded libraries from Ceramic Age individuals did not reach our 10% cytosine-to-thymine substitution threshold but had at least an 8% substitution rate, and therefore assessed as authentic given the relatively recent dates for these individuals; all eight libraries also were within the expected range for the other two authenticity metrics and had <1% contamination as assessed by contamLD. Multiple libraries fromI10333 andI10334 as well as one library fromI12341 showed poor match rates to the mtDNA consensus sequence, but this is likely due to low mtDNA coverage (0.5–2.1×). Two libraries from I7977 and one fromI15596 were also slightly below this threshold (6–10% mismatch rate), but also surpassed thresholds for the other two metrics and had ~1.1% contamination as assessed by contamLD.

We determined SNPs by randomly sampling an overlapping read with minimum mapping quality of ≥10 and base quality of ≥20. Individuals with <20,000 covered SNPs were excluded from quantitative analyses. One individual from each of three pairs of first-degree relatives in the dataset was excluded from population genetics analysis; in all cases, we retained the higher coverage individual; seeSupplementary Data 1.

We also generated shotgun sequencing data for two Ceramic-associated individuals from The Bahamas,I14922 (Abaco Island) andI14879 (South Andros) using the same system of data generation and processing, although the capture step was not included (Supplementary Data 2). For shotgun data, we report thresholds of mapping quality ≥30 and base quality ≥ 20.

Radiocarbon dates

We report 45 new radiocarbon (¹⁴C) dates on bone fragments generated using accelerator mass spectrometry (AMS) (Supplementary Data 3). Most dates (n=41) were generated at the Pennsylvania State University (PSU) Radiocarbon Laboratory, and the remainder (n=4) were generated at the Center for Isotopic Research on Cultural and Environmental heritage (CIRCE). The sample preparation methodology at PSU was carried out as previously reported²², where bone collagen was extracted and purified using a modified Longin method with ultrafiltration⁵⁵ (>30 kDa gelatin); if collagen yields were low, a modified XAD process⁵⁶ (XAD amino acids) was used. Carbon and nitrogen isotope ratios were then measured (Supplementary Information section 3) as a quality control measure; all C:N ratios fell between 3.15 and 3.44, indicating good collagen or amino acid preservation⁵⁵. We also evaluated diet in these individuals (e.g., marine vs. terrestrial) and compared the results to reference data from 242 ancient Caribbean and Maya individuals (Figures S12-S14). Attenuated Total Reflectance Fourier Transform Infrared (ATR-FTIR) spectra were generated to assess postmortem changes in the apatite crystal structure of the bone samples; ATR-FTIR spectra of all samples are displayed inFigure S15 and quality control parameters are reported inTable S1. Ultimately, all calibrated¹⁴C ages were computed using OxCal v4.4⁵⁷ using the IntCal20⁵⁸ after our stable isotope analysis detected minimal consumption of marine resources. Sample preparation at CIRCE was carried out following the lab-adapted Longin method⁵⁹; isotopic information was not generated for these individuals.Supplementary Data 3 lists the preparation method used for each individual andSupplementary Information section 3 describes the generation of isotopic data in more detail and its use in calibrating the¹⁴C dates generated for the Caribbean individuals.

Dataset assembly

We merged genome-wide data for 93 previously-reported individuals⁴ with newly-generated data from 174 ancient individuals for co-analysis, retaining 89 of them for a final co-analysis dataset comprising 263 individuals (details of merging inSupplementary Information section 4). We leverage these previously published data to revisit statistics and analyses reported in that work⁴ (Tables S2,S23,S29) and carry out additional analyses using these data (Tables S3,S24,S25,S26,S27,S28,Figures S33,S34).

We merged these 263 ancient individuals that passed screening into a base dataset that included 61 previously published ancient American individuals^{16,20,60–63}, and 36 modern Indigenous American groups sourced from single nucleotide polymorphism (SNP) array genotyping datasets or whole genome sequencing datasets (Supplementary Data 5):

• ‘1240K SNPs’, whole genome sequencing data restricted to a canonical set of 1,233,013 SNPs^47-50,64,65

• ‘Human Origins dataset’, 597,573 SNPs^66–68

• ‘Illumina dataset’ (unmasked/unadmixed individuals only), 352,432 SNPs¹³

All comparative analyses involving present-day Indigenous American populations were performed on the Illumina dataset, whereas forqpAdm andqpWave’s set of outgroup populations (“Right”) we used the Human Origins dataset for increased coverage. All genome-wide analyses were performed on autosomal data.

Uniparental haplogroups

We determined mtDNA haplogroups for all individuals using bam files, restricting to reads with MAPQ ≥ 30 and base quality ≥ 20. We constructed a consensus sequence with samtools and bcftools version 1.3.1 using a majority rule and then determined the haplogroup with HaploGrep2, using Phylotree version 17. We determined Y chromosome haplogroups using sequences mapping to 1240K Y-chromosome targets, restricting to sequences with MAPQ ≥ 30 and base quality ≥ 30. We called haplogroups by determining the most derived mutation for each individual, using the nomenclature of the International Society of Genetic Genealogy (ISOGG;http://www.isogg.org) version 14.76 (April 2019). Mutational differences and corresponding mtDNA haplogroups, and Y chromosome haplogroups and their supporting derived mutations are found inSupplementary Data 9. A discussion of mtDNA and Y chromosome haplogroup distribution in the Caribbean is found inSupplementary Information section 10; seeFigures S29 for distribution of mtDNA haplogroups,Figure S30 for details of three mtDNA mutations diagnostic of previously unobserved mtDNA haplogroup which is a variant of C1d, andFigure S31 for distribution of Y chromosome haplogroups.

Kinship

We assessed kinship for every pair of individuals newly-reported here as those that we co-analyze⁴ (including individuals from different sites and islands) using a previously described method⁶⁹, and we present results for 1st-, 2nd-, and 3rd-/4th-degree (‘close’) relatives inTable S5 (Supplementary Information section 7). In our newly-reported dataset of 174 ancient individuals, we identified 49 individuals sharing 49 unique pairwise kin relationships. Three pairs of individuals were identified as 1st-degree relatives, while 21 pairs were 2nd-degree relatives, and 25 pairs were 3rd-degree or higher. For the data that we co-analyze⁴, we identified 13 individuals who were part of eight relationships (four 2nd-degree and four 3rd-degree or higher). No close relatives were identified between the datasets. Distant cousins detected using IBD analysis are presented elsewhere (Extended Data Table 2;Supplementary Data 13).

Analysis of shared genomic segments

We identified Runs of Homozygosity (ROH) within our ancient dataset using the Python packagehapROH (https://test.pypi.org/project/hapsburg/). Following a previously described method²⁴, we used 5008 global haplotypes from the 1000 Genomes Project haplotype panel³³ as the reference panel. As recommended for datasets with genotypes for 1240K SNPs, we applied our method to ancient individuals with at least 400,000 SNPs covered and ran the method on the pseudo-haploid data to identify ROH longer than 4 centimorgan (cM). We used the default parameters ofhapROH, which are optimized for ancient data genotyped at 1240K SNPs. For each individual, we group the inferred ROH into four length categories: 4–8cM, 8–12cM, 12–20cM and >20cM and report the total sum in these bins (Supplementary Data 12;Fig. S21).

To estimate effective population size N_e from ROH, we applied a maximum likelihood inference framework (for derivation of the likelihood seeSupplementary Information section 7). We fit the lengths of all genome-wide ROH lengths 4–20cM, and infer the effective population size that maximizes the likelihood for ROH lengths observed in a set of individuals. Estimation uncertainties are obtained from the likelihood profile (95% CIs correspond to values within 1.92 units down from the maximum of the log-likelihood function). Tests on simulated data confirmed the ability of our estimator to recover N_e estimates from genome-wide ROH of few individuals (Figs. S22,S23).

We also analyzed shared genomic segments on the X chromosome between pairs of male individuals (“IBD_X”). To call such IBD blocks, we paired pseudo-haploid data of two X chromosomes and ranhapROH on read counts of the resulting artificial diploid individual; seeFigure S24 for example of IBD segment shared between two individuals. We inferred population sizes from IBD with the same likelihood approach as described for ROH, applying it to all pairs of individuals between two groups of individuals. SeeSupplementary Information section 7 for details.

Conditional Heterozygosity

We used popstats⁶⁸ to compute conditional heterozygosity for all clades and sub-clades, which we compared with contemporaneous groups from continental South America, such as from the Peruvian Middle and Late Horizon periods⁷⁰. As previously described^71,72, we restricted the analysis to transversion SNPs ascertained in a Yoruba individual; seeExtended Data Fig. 5.

PCA

We performed principal component analysis (PCA) with smartpca v181623⁷³, using the 1240K + Illumina merged dataset and using the option ‘lsqproject: YES’ to project ancient individuals onto the eigenvectors computed from modern individuals in the version shown in the main manuscript. The approach of projecting each ancient individual onto patterns of variation learned from modern individuals enables us to use data from a large fraction of SNPs covered in each individual and thereby maximize the information about ancestry that would be lost in approaches that require restriction to a potentially smaller number of SNPs for which there is intersecting data across lower coverage ancient individuals. We used the option ‘newshrink: YES’ to remap the points for the individuals used to generate the PCA onto the positions where they would be expected to fall if they had been projected, thereby allowing the projected and non-projected individuals to be appropriately co-visualized. We projected 92 previously published ancient individuals^4,16,20 and 174 new ancient individuals onto the first two principal components computed using 61 individuals from 23 present-day populations (Extended Data Fig. 1b). SeeSupplementary Data 4 for all individuals included in PCA and values of PCs 1 and 2 for the main manuscript PCA. For the PCA presented asFig. S19 (Supplementary Information section 5), we used non-related, non-outlier ancient individuals fromCuba_Archaic, Venezuela_Ceramic,EasternGreaterAntilles_Ceramic, BahamasCuba_Ceramic, andSECoastDR_Ceramic with >500K SNPs to compute the eigenvectors and projected all other ancient individuals. We again used the ‘lsqproject: YES’ and ‘newshrink: YES’ options. Individuals used to compute eigenvectors are listed inSupplementary Data 4. For PCA by archaeological site, non-zoomed PCA, PCA excluding CpG sites, and PCA with axes computed using ancient individuals, seeFigs. S16-S19.

Unsupervised analysis of population structure

We used the software ADMIXTURE v1.3.0^74,75 to perform unsupervised structure analysis on a dataset comprised of autosomal SNPs that overlap between the 1240k and Illumina dataset and pruned in PLINK1.9⁷⁶ using --indep-pairwise 200 25 0.4. This left 273,245 SNPs for the analysis. We ran five random-seeded replicates for each K in the interval between 2 and 10 with cross-validation enabled (--cv flag) to identify the runs with the low cross-validation errors (Table S4). For each value of K, we plotted the replicate with the lowest cross-validation error and compared the results. We choose to present K=6 asExtended Data Fig. 1c, as we found that the model with six components had a low cross-validation error and differentiated the components in a useful way for visualization. Results for the other values of K are presented asFig. S20 inSupplementary Information section 6.

Estimation of F_ST coefficients

To measure pairwise genetic differentiation between two groups of individuals, we estimated average pairwise F_ST and its standard error via block-jackknife using smartpca v.181623 and the options ‘fstonly: YES’ and ‘inbreed: YES.’ We removed the individual with lower coverage of each pair of first degree relatives, as well as ancestry outliers (see main text); we excludedHaiti_Ceramic, which comprises only two individuals who share a second-degree relationship as well as Macao, a site in the Dominican Republic from which all four individuals analyzed are 2nd-3rd-degree relatives of at least one other individual from the site. See results inExtended Data Fig. 2.

Clade grouping framework withqpWave,TreeMix andf₄-statistics

We used a multi-step framework involvingqpWave,TreeMix, andf₄-statistics to group sites and individuals, and considered this information together with admixture profiles and proportions fromqpAdm to produceFig. 1b (detailed methodology inSupplementary Information section 8). We started by usingqpWave to identify major clades based on shared ancestry and then usedTreeMix andf₄-statistics to investigate the existence of sub-clades. Once all sub-clades were identified, we usedf₄-statistics to investigate further substructure between sites within each clade. Geographic and chronological information such as island or cultural affiliation was not considered for these analyses, ensuring all clades and subclades were based solely on genetic information. We examined the association between genetic data and archaeological cultural complexes only after considering the genetic and archaeological information separately, following a previously published example⁷⁷.

The softwareqpWave¹³ from ADMIXTOOLS v6.0⁶⁸ estimates the minimum number of ancestry sources needed to form a group of test populations (“Left”), relative to a set of differentially related reference populations (“Right”). If the “Left” group contains two populations,qpWave will evaluate if they can be modelled as descending from the same sources, and hence will determine whether they form a clade. We used 12 present-day Indigenous American populations from the Human Origins dataset⁶⁷ plus Yukpa⁶⁴ representing different language families and ancestries from the American continent as our “Right” reference population set:

Chipewyan, Zapotec, Mixe, Mixtec, Suruí, Cabécar, Piapoco, Karitiana, Yukpa, Quechua, Wayuu, Apalai, Arara

The argument ‘allsnps: NO’ was used, which restricts the analysis SNP set to intersection of all SNPs among all populations and maximizes the reliability of the analysis⁷⁸. The ‘allsnps: YES’ option was developed to increase the number of SNPs analyzed in cases where very little SNP overlap exists between all populations included in aqpWave model⁷⁹. While it is commonly used when low coverage data results in the loss of the majority of sites in the initial datasets⁷⁸, there is a risk that this option introduces unreliability in the analysis, particularly in cases where the base population is highly diverged. In this dataset, a high depth of coverage and relatively large sample sizes made it unnecessary for us to use the ‘allsnps: YES’ option. We ran two consecutive steps ofqpWave analyses, starting with the identification of major groupings (step 1;Figure S25), or clades, and then reassessed the relationships between members within those clades by running the same tests in a “model competition” approach where individuals from other sites from within the same clade were added to the “Right” set (step 2;Figure S26). A significance threshold of p>0.01 was set for accepting a clade between two sites or individuals. The range of covered SNPs was 170,927–827,039, with a median of 672,888.

After identifying the major clades and/or pairs of sites that uniquely formed a clade with one another, we ranTreeMix with these clades and 27 previously published present-day Indigenous populations¹³ (Supplementary Data 5) to identify within-clade site structure (step 3;Figures S27,S28) by generating a maximum likelihood tree. We excluded four Chibchan, Chocoan and Arawak-speaking populations possibly admixed with each other from this analysis. We ranTreeMix, grouping the SNPs in windows of 500 (flag -k 500) to account for linkage disequilibrium, setting Chipewyan as root (-root), allowing random migration events (-m), and disabling sample size correction (-noss) in order to include sites or populations represented by a single-individual. We note that single-individual populations still present artifactually long branches that do not truly represent population-specific drift. By runningTreeMix and allowing consecutive random migration/admixture events, we identified nodes and branches that maintained the same ancient Caribbean sites among the different runs. We then usedf₄-statistics to evaluate if they formed a sub-clade to the exclusion of the other sites by following the tree’s structure. For each identified intact node among allTreeMix runs we used each downstream pair of site(s) as Test1 and Test2 and investigated their relationship to upstream sites or pools of sites (step 4). If an upstream node was unchanged in all runs, the sites composing it were pooled. However, once the first inconsistency was identified in an upstream node, all sites beyond that node were pooled together. A combination of three statistics per relationship allowed us to evaluate theTreeMix structure of the sites being tested:

With Test1 and Test2 expected to be closer to each other than to Pool, the tested relationship finds support if the first test is statistically non-significant and at least one of the other two are significant. We used a Z-score threshold of 2.8 (associated with a 99.5% CI) to assess significance. These sites were then merged into a sub-clade inside the major Ceramic clade for further analysis. We did not include the sites of Cueva del Perico I, Los Indios, Punta Candelero, and Tibes in theTreeMix andf₄ due to reduced coverage, but evaluated these sites separately to see if they shared closer affinities to any sub-clades relative to the others (Supplementary Data 7;Supplementary Information section 8).

After this clading analysis, we usedf₄-statistics to further investigate potential substructure between sites within each sub-clade (step 5). For each pairwise site comparison, we randomly divided each site into two groups of individuals, and used a statistic of the formf₄(Site1_subset1, Site2_subset1; Site1_subset2, Site2_subset2) to identify positive statistics suggesting substructure within the same clade. This randomization step was repeated 10 times, and the average Z-score was calculated. If a site was composed of a single individual we instead computed statistics of the formf₄(Mbuti,Site1_subset1; Site2_singleIndividual, Site1_subset2), intended to evaluate if individuals within Site1 were closer to each other than to the single individual from Site2. No statistics were computed if both sites being tested contained only one individual.

We also usedf₄-statistics to test if any specific sub-clade within theCaribbean_Ceramic clade had more Archaic-related ancestry than another. Specifically we used the statisticf₄(Mbuti, GreaterAntilles_Archaic, Sub_Clade1, Sub_Clade2) and interpreted results as significant based on a |Z|>2.8; results are presented inTable S20.

Admixture simulations

We investigated the sensitivity ofqpWave in detecting Carib-related ancestry in theCaribbean_Ceramic sub-clades by generating artificially admixed individuals withCaribbean_Ceramic ancestry mixed with increasing amounts (1, 2, 5, 8, 10, 20, 30, 40, and 50%) of a plausibly Carib-associated ancestry. For the Carib-associated ancestry we tested Arara (present-day Indigenous Carib speakers),Venezuela_Ceramic (inhabitants of a possible region of origin for this ancient Carib migration), and alsoLesserAntilles_Ceramic (possibly representing Island Caribs), and then assessed at what admixture threshold we were able to reliably detect the latter ancestry type (Supplementary Information section 13;Fig. S32). To generate these admixed individuals, we identified common SNPs between the two sources, randomly selected genotypes from the Arara individuals from the Human Origins and Illumina SNP array datasets corresponding to each of the nine percentages to be tested, and added the remaining SNPs from a random individual fromBahamas_Ceramic,EasternGreaterAntilles_Ceramic,SECoastDR_Ceramic, andLesserAntilles_Ceramic with over 800,000 SNPs. We then ranqpWave with each of the simulated admixed individuals on the “Left” plus their correspondent sub-clade, while using the default 12 “Right” populations (excluding Arara), as described inSupplementary Information section 8, plus the Carib proxy population used to generate those individuals.

Dating admixture

We used the methodDATES (Distribution of Ancestry Tracts of Evolutionary Signals²² v3520 (Chintalapati, M., Neel, A., Patterson, N. & Moorjani, P. Reconstructing the spatio-temporal patterns of admixture in human history.In Preparation.) to estimate the dates of admixture in admixed individuals from Haiti. This method measures the decay of ancestry covariance to infer the time since mixture and estimates jackknife standard errors. Details ofDATES analysis are found inSupplementary Information section 14; results forHaiti_Ceramic are found inTable S22.

Relatedness of ancient individuals to present-day admixed Caribbean populations

We computed relative allele-sharing between present-day admixed Caribbean populations (via their Indigenous ancestry) and ancient Archaic-associated versus Ceramic-associated individuals with ADMIXTOOLS 2 (Maier R., Reich D., Patterson N. Rapid inference of demographic history using ADMIXTOOLS 2.In Preparation.) through the statisticf₄(European,Test;Cuba_Archaic,Caribbean_Ceramic). In order to evaluate statistical power, we compared results for present-day Cubans alone to results obtained by adding one ancient individual from either theGreaterAntilles_Archaic orCaribbean_Ceramic clade to the Cuban test population. Full details are found inSupplementary Information section 15.

Analysis of phenotypically-relevant SNPs

Analyzing SNPs previously known to be relevant to phenotypic traits allows us to explore their frequencies in the pre-contact Caribbean and Venezuela. We used mpileup insamtools⁸⁰ version 1.3.1 with the settings -B -q30 -Q30 to obtain information about each SNP covered by reads from the bam files of our individuals (after trimming 2 base pairs from the molecule ends) and used the fasta file from human genome GRCh37 (hg19) as a reference file for the pileup. We counted the number of reference and alternate alleles, combining counts on the forward and reverse strands. Data are provided inSupplementary Data 15, with a discussion of results inSupplementary Information section 16.

Testing for an Australasian link

We tested for a signal of relatedness to present-day Australasian populations^64,68 (“Population Y” signal), using the statisticf₄(Mbuti,Onge/Papuan;Mixe,Archaic/Ceramic) and testing all final sub-clades asArchaic/Ceramic. Here, Mixe is representative of a population that harbors no Population Y signal. When Onge was used as the Australasian proxy, several of the ancient groups showed weakly positive statistics (Z between 2 and 3), but only the Archaic individualI10126 from the site of Andrés (Dominican Republic) was significant at Z = 3.4. While this signal is significant at p=0.0030 even after performing a Bonferroni correction for the nine hypotheses tested inExtended Data Table 4, the signal is non-significant when Papuan is used as the Australasian proxy (Z=2.2). We also caution that all Population Y statistics are likely to be overinflated in their significance because the original discovery of the Population Y signal carried out extensive hypothesis testing to identify a population in the third position of the statisticf₄(Mbuti,Onge/Papuan;Mixe,Archaic/Ceramic) (Mixe) that maximized the value of the statistic when any other Native American group in was used in the fourth position; thus, there is a further multiple hypothesis testing issue for which our analysis does not correct. The lack of a clear population Y signal is consistent with prior studies that also have not found this signal in ancient individuals from this region¹⁶ and other areas of South America⁶³.

Extended Data Table 1: N_e estimates for each site.

Extended Data Table 2: Subset of cross-site relatives from different islands, identified through IBD analysis.

Extended Data Table 3: Ancestry proportion estimates withqpAdm in present-day Caribbean individuals from Cuba (and its provinces), Dominican Republic, and Puerto Rico^21,28.

Extended Data Table 4:

• 1.Rouse I. The Tainos: Rise & Decline of the People who Greeted Columbus. (Yale University Press, 1992). [Google Scholar]
• 2.Maggiolo MVLa isla de Santo Domingo antes de Colón (Banco Central de la Republica Dominicana, 1993). [Google Scholar]
• 3.Keegan WF & Hofman CLThe Caribbean before Columbus. (Oxford University Press, 2017). [Google Scholar]
• 4.Nägele K. et al. Genomic insights into the early peopling of the Caribbean. Science369, 456–460 (2020). [DOI] [PubMed] [Google Scholar]
• 5.Cook SF & Borah W. The Aboriginal Population of Hispaniola. vol. 1376–410 (University of California Press, 1971). [Google Scholar]
• 6.Henige D. On the Contact Population of Hispaniola: History as Higher Mathematics. Hispanic American Historical Review58, 217–237 (1978). [Google Scholar]
• 7.Wilson SMThe Archaeology of the Caribbean. (Cambridge University Press, 2007). [Google Scholar]
• 8.Rodríguez Ramos R. Isthmo–Antillean Engagements. in Oxford Handbook of Caribbean Archaeology (eds. Keegan WF, Hofman CL & Rodríguez Ramos R.) 155–170 (Oxford University Press, 2013). [Google Scholar]
• 9.Bérard B. About boxes and labels: A periodization of the Amerindian occupation of the West Indies. Journal of Caribbean Archaeology19, 51–67 (2019). [Google Scholar]
• 10.Callaghan RTArchaeological Views of Caribbean Seafaring in Oxford handbook of Caribbean archaeology (eds. Keegan WF, Hofman C. & Rodriguez RR) 285–295 (Oxford University Press, 2013). [Google Scholar]
• 11.Siegel PEet al. Paleoenvironmental evidence for first human colonization of the eastern Caribbean. Quaternary Science Reviews129, 275–295 (2015). [Google Scholar]
• 12.Oliver JRThe archaeological, linguistic and ethnohistorical evidence for the expansion of Arawakan into northwestern Venezuela and northeastern Colombia. (University of Illinois at Urbana-Champaign, 1989). [Google Scholar]
• 13.Reich D. et al. Reconstructing Native American population history. Nature488, 370–374 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 14.Greenberg JHLanguage in the Americas. (Stanford University Press, 1987). [Google Scholar]
• 15.Salzano FM, Hutz MH, Salamoni SP, Rohr P. & Callegari‐Jacques SMGenetic Support for Proposed Patterns of Relationship among Lowland South American Languages. Current Anthropology46, S121–S128 (2005). [Google Scholar]
• 16.Schroeder H. et al. Origins and genetic legacies of the Caribbean Taino. Proc. Natl. Acad. Sci. U. S. A115, 2341–2346 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 17.Pickrell JK & Pritchard JKInference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8, e1002967 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 18.Chinique de Armas Y., Roksandic M., Suárez RR, Smith DG & Buhay WMIsotopic Evidence of Variations in Subsistence Strategies and Food Consumption Patterns among ‘Fisher-Gatherer’ Populations of Western Cuba in Cuban Archaeology in the Circum-Caribbean Context (ed. Roksandic I.) (University Press of Florida, 2016). [Google Scholar]
• 19.Lovén SEOrigins of the Tainan Culture, West Indies. (Elanders Bokfryckeri Akfiebolag, 1935). [Google Scholar]
• 20.Nieves-Colón MAet al. Ancient DNA reconstructs the genetic legacies of pre-contact Puerto Rico communities. Molecular Biology and Evolution37, 611–626 (2020). [DOI] [PubMed] [Google Scholar]
• 21.Moreno-Estrada A. et al. Reconstructing the population genetic history of the Caribbean. PLoS Genet. 9, e1003925 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 22.Narasimhan VMet al. The formation of human populations in South and Central Asia. Science365, (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 23.Ross AH, Keegan WF, Pateman MP & Young CBFaces Divulge the Origins of Caribbean Prehistoric Inhabitants. Sci. Rep10, 147 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 24.Ringbauer H., Novembre J. & Steinrucken M. Detecting runs of homozygosity from low-coverage ancient DNA.bioRxiv.org (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 25.Ceballos FC, Joshi PK, Clark DW, Ramsay M. & Wilson JFRuns of homozygosity: windows into population history and trait architecture. Nat. Rev. Genet19, 220–234 (2018). [DOI] [PubMed] [Google Scholar]
• 26.Frankham R. Effective population size/adult population size ratios in wildlife: a review. Genet. Res89, 491–503 (2007). [DOI] [PubMed] [Google Scholar]
• 27.Browning SR & Browning BLAccurate Non-parametric Estimation of Recent Effective Population Size from Segments of Identity by Descent. Am. J. Hum. Genet97, 404–418 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 28.Fortes-Lima C. et al. Exploring Cuba’s population structure and demographic history using genome-wide data. Sci. Rep8, 11422 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 29.Toro-Labrador G., Wever OR & Martínez-Cruzado JCMitochondrial DNA Analysis in Aruba: Strong Maternal Ancestry of Closely Related Amerindians and Implications for the Peopling of Northwestern Venezuela. Caribbean Journal of Science39, (2003). [Google Scholar]
• 30.Mendizabal I. et al. Genetic origin, admixture, and asymmetry in maternal and paternal human lineages in Cuba. BMC Evol. Biol8, 213 (2008). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 31.Vilar MGet al. Genetic diversity in Puerto Rico and its implications for the peopling of the Island and the West Indies. Am. J. Phys. Anthropol155, 352–368 (2014). [DOI] [PubMed] [Google Scholar]
• 32.Benn Torres J. et al. Genetic Diversity in the Lesser Antilles and Its Implications for the Settlement of the Caribbean Basin. PLoS One10, e0139192 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 33.Consortium T. 1000 G. P. & The 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature526, 68–74 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 34.Hofman CL & Reid BAThe Saladoid. in Encyclopedia of Caribbean archaeology (eds. Reid B. & Gilmore G.) 300–303 (University of Florida Press, 2014). [Google Scholar]
• 35.Roksandic I. & Roksandic M. Peopling of the Caribbean. 199–223 (Kerns: Verlag, 2018). [Google Scholar]
• 36.Keegan W. The People Who Discovered Columbus. (University Press of Florida, 1992). [Google Scholar]
• 37.Anderson-Córdova KFHispaniola and Puerto Rico: Indian Acculturation and Heterogeneity, 1492–1550. (University Microfilms International, 1990). [Google Scholar]
• 38.Pinhasi R. et al. Optimal Ancient DNA Yields from the Inner Ear Part of the Human Petrous Bone. PLoS One10, e0129102 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 39.Pinhasi R., Fernandes DM, Sirak K. & Cheronet O. Isolating the human cochlea to generate bone powder for ancient DNA analysis. Nat. Protoc14, 1194–1205 (2019). [DOI] [PubMed] [Google Scholar]
• 40.Sirak K. et al. Human auditory ossicles as an alternative optimal source of ancient DNA. Genome Res. 30, 427–436 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 41.Dabney J. et al. Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc. Natl. Acad. Sci. U. S. A110, 15758–15763 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 42.Korlević P. et al. Reducing microbial and human contamination in DNA extractions from ancient bones and teeth. Biotechniques59, 87–93 (2015). [DOI] [PubMed] [Google Scholar]
• 43.Rohland N., Glocke I., Aximu-Petri A. & Meyer M. Extraction of highly degraded DNA from ancient bones, teeth and sediments for high-throughput sequencing. Nat. Protoc13, 2447–2461 (2018). [DOI] [PubMed] [Google Scholar]
• 44.Rohland N., Harney E., Mallick S., Nordenfelt S. & Reich D. Partial uracil-DNA-glycosylase treatment for screening of ancient DNA. Philos. Trans. R. Soc. Lond. B Biol. Sci370, 20130624 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 45.Gansauge M-T, Aximu-Petri M., Nagel K. & Meyer M. Manual and automated preparation of single-stranded DNA libraries for the sequencing of DNA from ancient biological remains and other sources of highly degraded DNA. Nature Protocols15, 2279–3000 (2020). [DOI] [PubMed] [Google Scholar]
• 46.Briggs AWet al. Removal of deaminated cytosines and detection of in vivo methylation in ancient DNA. Nucleic Acids Res. 38, e87 (2010). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 47.Fu Q. et al. A revised timescale for human evolution based on ancient mitochondrial genomes. Curr. Biol23, 553–559 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 48.Fu Q. et al. An early modern human from Romania with a recent Neanderthal ancestor. Nature524, 216–219 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 49.Haak W. et al. Massive migration from the steppe was a source for Indo-European languages in Europe. Nature522, 207–211 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 50.Mathieson I. et al. Genome-wide patterns of selection in 230 ancient Eurasians. Nature528, 499–503 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 51.Behar DMet al. A ‘Copernican’ reassessment of the human mitochondrial DNA tree from its root. Am. J. Hum. Genet90, 675–684 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 52.Li H. & Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics26, 589–595 (2010). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 53.Korneliussen TS, Albrechtsen A. & Nielsen R. ANGSD: Analysis of Next Generation Sequencing Data. BMC Bioinformatics15, 356 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 54.Nakatsuka N. et al. ContamLD: estimation of ancient nuclear DNA contamination using breakdown of linkage disequilibrium. Genome Biol. 21, 199 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 55.Kennett DJet al. Archaeogenomic evidence reveals prehistoric matrilineal dynasty. Nat. Commun8, 14115 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 56.Lohse JC, Culleton BJ, Black SL & Kennett DJA Precise Chronology of Middle to Late Holocene Bison Exploitation in the Far Southern Great Plains. Journal of Texas Archeology and History1, 94–126 (2014). [Google Scholar]
• 57.Ramsey CBBayesian Analysis of Radiocarbon Dates. Radiocarbon51, 337–360 (2009). [Google Scholar]
• 58.Reimer PJet al. The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0–55 cal kBP). Radiocarbon62, 725–757 (2020). [Google Scholar]
• 59.Passariello I. et al. Characterization of Different Chemical Procedures for 14C Dating of Buried, Cremated, and Modern Bone Samples at Circe. Radiocarbon54, 867–877 (2012). [Google Scholar]
• 60.Lindo J. et al. The genetic prehistory of the Andean highlands 7000 years BP though European contact. Sci Adv4, eaau4921 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 61.Moreno-Mayar JVet al. Early human dispersals within the Americas. Science362, (2018). [DOI] [PubMed] [Google Scholar]
• 62.Scheib CLet al. Ancient human parallel lineages within North America contributed to a coastal expansion. Science360, 1024–1027 (2018). [DOI] [PubMed] [Google Scholar]
• 63.Posth C. et al. Reconstructing the Deep Population History of Central and South America. Cell175, 1185–1197.e22 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 64.Raghavan M. et al. Genomic evidence for the Pleistocene and recent population history of Native Americans. Science349, aab3884 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 65.Mallick S. et al. The Simons Genome Diversity Project: 300 genomes from 142 diverse populations. Nature538, 201–206 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 66.Patterson N. et al. Ancient admixture in human history. Genetics192, 1065–1093 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 67.Lazaridis I. et al. Ancient human genomes suggest three ancestral populations for present-day Europeans. Nature513, 409–413 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 68.Skoglund P. et al. Genetic evidence for two founding populations of the Americas. Nature525, 104–108 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 69.Olalde I. et al. The genomic history of the Iberian Peninsula over the past 8000 years. Science363, 1230–1234 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 70.Nakatsuka N. et al. A Paleogenomic Reconstruction of the Deep Population History of the Andes. Cell (2020) doi: 10.1016/j.cell.2020.04.015. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 71.Skoglund P. et al. Genomic insights into the peopling of the Southwest Pacific. Nature538, 510–513 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 72.Harney Éet al. Ancient DNA from Chalcolithic Israel reveals the role of population mixture in cultural transformation. Nat. Commun9, 3336 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 73.Patterson N., Price AL & Reich D. Population structure and eigenanalysis. PLoS Genet. 2, e190 (2006). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 74.Alexander DH, Novembre J. & Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 75.Alexander DH & Lange K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics12, (2011). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 76.Chang CCet al. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience4, 7 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 77.Fu Q. et al. The genetic history of Ice Age Europe. Nature534, 200–205 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 78.Lipson M. Applying f4-Statistics and Admixture Graphs: Theory and Examples. Mol Ecol Resour00, 1–10 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 79.Harney É, Patterson N., Reich D. & Wakeley J. Assessing the Performance of qpAdm: A Statistical Tool for Studying Population Admixture. bioRxiv (2020) doi: 10.1101/2020.04.09.032664. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 80.Li H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics25, 2078–2079 (2009). [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials