Cesar A Fortes-Lima

Romain Laurent

Valentin Thouzeau

Bruno Toupance

Paul Verdu
email: paul.verdu@mnhn.fr

Complex genetic admixture histories reconstructed with Approximate Bayesian Computations

Keywords: Admixture, Approximate Bayesian Computation, Inference, Population Genetics, Machine Learning

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

Acknowledgements:

We warmly thank Frédéric Austerlitz, Erkan O. Buzbas, Antoine Cools, Flora Jay, Evelyne Heyer, Margueritte Lapierre, Guillaume Laval, Nina Marchi, Etienne Patin, Noah A.

Rosenberg, and Zachary A. Szpiech for useful comments and discussions. This project was funded in part by the French Agence Nationale de la Recherche project METHIS (ANR 15-CE32-0009-01). CFL was funded in part by the Sven and Lilly Lawski's Foundation (N2019-0040).

Authors declare no conflict of interest for this work.

Figure and Table Content:

Main Text: 5 figures, 4 tables. Supplementary material: 1 note, 6 figures, 3 tables.

Novel online resources:

MetHis software package can be downloaded with manual and example dataset from https://github.com/romain-laurent/MetHis

INTRODUCTION

Hybridization between species and admixture between populations are powerful mechanisms influencing biological evolution. Genetic admixture patterns have thus been extensively studied to understand migrations and admixture-related adaptation [START_REF] Heliconius Genome Consortium | Butterfly genome reveals promiscuous exchange of mimicry adaptations among species[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF][START_REF] Skoglund | Ancient wolf genome reveals an early divergence of domestic dog ancestors and admixture into high-latitude breeds[END_REF][START_REF] Brandenburg | Independent introductions and admixtures have contributed to adaptation of European maize and its American counterparts[END_REF]. The increasing availability of genome-wide data in numerous species, and particularly humans (e.g. 1000 GENOMES [START_REF] Genomes | A global reference for human genetic variation[END_REF], further provides unprecedented opportunities to understand the genomic architecture of admixture, characterize the contribution of admixture to adaptive evolution, and infer demographic histories of admixture from genetic data.

Based on a long history of statistical developments aimed at investigating admixture patterns from genetic data [START_REF] Bernstein | Die geographische Verteilung der Bludgruppen und ihre anthropologische Bedeutung[END_REF][START_REF] Cavalli-Sforza | The genetics of human populations[END_REF][START_REF] Chakraborty | Admixture as a tool for finding linked genes and detecting that difference from allelic association between loci[END_REF][START_REF] Long | The genetic structure of admixed populations[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF][START_REF] Patterson | Ancient admixture in human history[END_REF], population geneticists recently developed methods to reconstruct the genomic architecture of admixed segments deriving from each source population, and to describe admixture linkage-disequilibrium (LD) patterns [START_REF] Sankararaman | Estimating local ancestry in admixed populations[END_REF][START_REF] Price | Sensitive detection of chromosomal segments of distinct ancestry in admixed populations[END_REF][START_REF] Lawson | Inference of population structure using dense haplotype data[END_REF][START_REF] Maples | RFMix: a discriminative modeling approach for rapid and robust local-ancestry inference[END_REF][START_REF] Guan | Detecting structure of haplotypes and local ancestry[END_REF][START_REF] Salter-Townshend | Fine-Scale Inference of Ancestry Segments Without Prior Knowledge of Admixing Groups[END_REF]. In Homo sapiens, these methods have been extensively used to infer populations' ancestral genetic origins and map local ancestry along individual genomes, often for disease-mapping purposes (e.g. [START_REF] Shriner | Mapping of disease-associated variants in admixed populations[END_REF]. Furthermore, by coupling admixture mapping approaches with natural selection scans, sometimes accounting for ancient and recent demographic history, it is possible to identify signatures of adaptive introgression or post-admixture selection having influenced genomic diversity patterns in human populations [START_REF] Jeong | Admixture facilitates genetic adaptations to high altitude in Tibet[END_REF][START_REF] Racimo | Evidence for archaic adaptive introgression in humans[END_REF][START_REF] Patin | Dispersals and genetic adaptation of Bantu-speaking populations in Africa and North America[END_REF].

In this context, several maximum-likelihood approaches have been developed to estimate the parameters of admixture models (time of admixture events and their associated intensities) that vastly improved our understanding of detailed admixture histories in particular for human populations (e.g. [START_REF] Pickrell | Inference of population splits and mixtures from genome-wide allele frequency data[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF]. The two classes of methods most extensively deployed in the past rely, respectively, on the moments of allelic frequency spectrum divergences among populations [START_REF] Reich | Reconstructing Indian population history[END_REF][START_REF] Patterson | Ancient admixture in human history[END_REF][START_REF] Pickrell | Inference of population splits and mixtures from genome-wide allele frequency data[END_REF][START_REF] Lipson | Efficient moment-based inference of admixture parameters and sources of gene flow[END_REF], and on admixture LD patterns [START_REF] Pool | Inference of historical changes in migration rate from the lengths of migrant tracts[END_REF][START_REF] Moorjani | The history of African gene flow into Southern Europeans, Levantines, and Jews[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Loh | Inferring admixture histories of human populations using linkage disequilibrium[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF][START_REF] Chimusa | Dating admixture events is unsolved problem in multi-way admixed populations[END_REF]. They allow for identifying admixture events in a given set of populations, estimating admixture fractions, and inferring time since each pulse of admixture. Notably, [START_REF] Gravel | Population genetics models of local ancestry[END_REF] developed an approach to fit the observed curves of admixture LD decay to those theoretically expected under admixture models involving one or two possible pulses of admixture from multiple source populations. This major advance significantly improved our ability to reconstruct detailed admixture histories using genetic data, for instance among several populations descending from the Transatlantic Slave Trade (TAST) across the Americas (e.g. [START_REF] Moreno-Estrada | Reconstructing the population genetic history of the Caribbean[END_REF][START_REF] Baharian | The Great Migration and African-American Genomic Diversity[END_REF][START_REF] Fortes-Lima | Genome-wide Ancestry and Demographic History of African-Descendant Maroon Communities from French Guiana and Suriname[END_REF].

Despite the unquestionable importance of these previous developments, existing admixture history inference methods somewhat suffer from inherent limitations acknowledged by the authors [START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Lipson | Efficient moment-based inference of admixture parameters and sources of gene flow[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF]. First, most likelihood approaches can only consider one or two pulses of admixture in the history of the hybrid population. Nevertheless, admixture processes in numerous species are known to be often much more complex, involving multiple admixture-pulses or periods of recurring admixture over time from each source population separately. It is not yet clear how these methods might behave when they can consider only simplified versions of the true admixture history underlying the observed data [START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Lipson | Efficient moment-based inference of admixture parameters and sources of gene flow[END_REF][START_REF] Loh | Inferring admixture histories of human populations using linkage disequilibrium[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF][START_REF] Medina | Estimating the Timing of Multiple Admixture Pulses During Local Ancestry Inference[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. Second, while it is possible to compare maximumlikelihood values obtained from fitting one or two admixture pulses to the observed data as a guideline to find the "best" scenario, formal statistical comparison of model posterior probabilities is often out of reach of these approaches [START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Foll | WFABC: a Wright-Fisher ABC-based approach for inferring effective population sizes and selection coefficients from time-sampled data[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. Finally, admixture-LD methods, in particular, rely on fine mapping of local ancestry segments in individual genomes and thus require substantial amounts of genomic data (typically several hundred thousand to several millions of SNPs), and, sometimes, accurate phasing. These still represent major challenges for most species, including humans.

To overcome these limitations, Approximate Bayesian Computation (ABC) approaches [START_REF] Tavaré | Inferring coalescence times from DNA sequence data[END_REF][START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF][START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF] represent a promising class of methods to infer complex admixture histories from observed genetic data. Indeed, ABC has been successfully used previously in different species (including humans), and using different types of genetic data, to formally test alternative demographic scenarios hypothesized to be underlying observed genetic patterns, and to estimate, a posteriori, the parameters of the winning models [START_REF] Verdu | Origins and genetic diversity of pygmy hunter-gatherers from Western Central Africa[END_REF][START_REF] Boitard | Inferring Population Size History from Large Samples of Genome-Wide Molecular Data -An Approximate Bayesian Computation Approach[END_REF][START_REF] Fraimout | Deciphering the Routes of invasion of Drosophila suzukii by Means of ABC Random Forest[END_REF].

ABC model-choice and posterior parameter inference rely on comparing observed summary statistics to the same set of statistics, calculated from a usually large number of genetic simulations explicitly parametrized by the user, and produced under competing demographic scenarios [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Wegmann | Efficient approximate Bayesian computation coupled with Markov chain Monte Carlo without likelihood[END_REF][START_REF] Blum | Non-linear regression models for Approximate Bayesian Computation[END_REF][START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF][START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF][START_REF] Sisson | Handbook of Approximate Bayesian Computation[END_REF]. Each simulation, and corresponding vector of summary statistics, is produced using model-parameters drawn randomly from prior distributions informed adequately by the user. Therefore, the flexibility of ABC relies mostly on explicit genetic data simulations set by the user. This makes ABC a priori particularly well suited to investigate highly complex historical admixture scenarios for which likelihood functions are very often intractable, but for which simulation of genetic data is feasible [START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF][START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF]. However, ABC has until now seldom been used to investigate admixture processes beyond a single admixture pulse or constant migrations [START_REF] Buzbas | AABC: approximate approximate Bayesian computation for inference in population-genetic models[END_REF][START_REF] Buzbas | Inference on admixture fractions in a mechanistic model of recurrent admixture[END_REF].

In this paper, we show how ABC can be successfully applied to reconstruct, from genetic data, highly complex admixture histories beyond exploring models with a single or two pulses of admixture. In particular, we focus on evaluating how a relatively limited number of independent SNPs can be used for accurately distinguishing major classes of historical admixture models, such as multiple admixture-pulses versus recurring increasing or decreasing admixture over time, and for conservative posterior parameter inference under the winning model. Furthermore, we show that the quantiles and higher moments of the distribution of admixture fractions in the admixed population are highly informative summary-statistics for ABC modelchoice and posterior-parameter estimation, as expected analytically [START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Buzbas | Inference on admixture fractions in a mechanistic model of recurrent admixture[END_REF].

In order to do so, and since genetic data simulation under highly complex admixture models is not trivial using existing coalescent approaches [START_REF] Wakeley | Gene genealogies within a fixed pedigree, and the robustness of Kingman's coalescent[END_REF] We exemplify our approach by reconstructing the complex admixture histories underlying observed genetic patterns separately for the African American (ASW) and Barbadian (ACB) populations from the 1000 Genomes Project Phase 3 (1000 GENOMES [START_REF] Genomes | A global reference for human genetic variation[END_REF]. Both populations are known to be admixed populations of European and African descent in the context of the TAST (e.g. GRAVEL 2012;[START_REF] Baharian | The Great Migration and African-American Genomic Diversity[END_REF][START_REF] Martin | Human Demographic History Impacts Genetic Risk Prediction across Diverse Populations[END_REF].

We find admixture histories much more complex than previously inferred for these populations and further reveal the diversity of admixture histories undergone by populations descending from the TAST in the Americas.

MATERIAL AND METHODS

We aimed at evaluating how ABC model-choice and posterior parameter estimation could allow reconstructing highly complex historical admixture processes using independent genomewide SNPs. To do so, we chose to focus on the recent admixture history of populations of African and European ancestry, descending from European colonization and the TAST in the Americas. This case-study represents an appropriate setting for empirically testing our ABC approach, since this period of history starting in the late 15th century has been extensively studied in population genetics based on the same publicly available datasets.

First, we describe the targeted case-study population and genetic datasets. Second, we present in detail the complex admixture processes here investigated and the associated demographic parameters. Third, we describe the novel simulation and summary statistics calculation software package called MetHis, here proposed to investigate these admixture processes. Fourth, we detail the Random-Forest ABC procedure used for scenario-choice inference and the performance of this approach both in general for the tested models and specifically for the real data here investigated. Finally, we detail the Neural Network ABC procedure deployed to estimate posterior parameter distributions, its parameterization, and the cross-validation procedures conducted to evaluate its power and accuracy.

Population Genetics Dataset

We considered the admixture histories of the African American (ASW) and Barbadian (ACB) population samples from the 1000 Genomes Project Phase 3 (1000 GENOMES PROJECT CONSORTIUM 2015). Previous studies identified, within the same database, the West European Great-Britain (GBR) and the West African Yoruba (YRI) population samples as reasonable proxies for the genetic sources of the admixture of both ACB and ASW populations, consistently with the macro-history of the TAST in the former British colonial empire in Africa and the Americas [START_REF] Baharian | The Great Migration and African-American Genomic Diversity[END_REF][START_REF] Martin | Human Demographic History Impacts Genetic Risk Prediction across Diverse Populations[END_REF][START_REF] Verdu | Parallel Trajectories of Genetic and Linguistic Admixture in a Genetically Admixed Creole Population[END_REF].

We excluded from our sample set, individuals previously identified to be more closely related than first-degree cousins in the four populations separately [START_REF] Verdu | Parallel Trajectories of Genetic and Linguistic Admixture in a Genetically Admixed Creole Population[END_REF]. We also excluded the three ASW individuals showing traces of Native American or East-Asian admixture beyond that from Europe and Africa, as reported in previous studies [START_REF] Martin | Human Demographic History Impacts Genetic Risk Prediction across Diverse Populations[END_REF]. This allows us to consider only two source populations for the admixture history of both admixed populations investigated here. Among the remaining individuals we randomly drew 50 individuals in the targeted admixed ACB and ASW populations, respectively, and included the remaining 90 YRI individuals and 89 GBR individuals.

We extracted biallelic polymorphic sites (SNPs as defined by the 1000 Genomes Project Phase 3) from the merged ACB+ASW+GBR+YRI data set, excluding singletons. Furthermore, we focused only on independent SNPs by LD pruning the data set using the PLINK [START_REF] Purcell | PLINK: a tool set for wholegenome association and population-based linkage analyses[END_REF]) --indep-pairwise option with a sliding window of 100 SNPs, moving in increments of 10 SNPs, and r 2 threshold of 0.1 [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]. Finally, we randomly drew 100,000 SNPs from the remaining SNP set.

Competing complex admixture scenarios

We aimed at investigating comprehensive admixture histories with, after the original foundation of the admixed population, possibly multiple pulses (>1) of admixture, or recurring monotonically increasing or decreasing admixture, from each source population separately. To do so, we chose to work under the general mechanistic model presented in [START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF], henceforth called the VR2011 model, derived from [START_REF] Ewens | The transmission/disequilibrium test: history, subdivision, and admixture[END_REF]. Briefly (Supplementary Figure S1), the VR2011 general model considers, for diploid organisms, a panmictic admixture process, discrete in generations, where M source populations Sm contribute to the hybrid population H at the following generation g + 1 with proportions 𝑠 𝑚,𝑔 each in [0,1], and where the hybrid population H contributes to itself with proportion ℎ 𝑔 in [0,1] with ℎ 0 = 0, satisfying, for each value of g ≥ 0, ∑ 𝑠 𝑚,𝑔 + ℎ 𝑔 𝑚∈ [1,𝑀] =1.

Here, we adapted the two source-populations version of the general VR2011 (M = 2), and define, next, the nine competing complex admixture scenarios considered to reconstruct the history of introgression from Africa and Europe into the gene-pool of the ACB and ASW admixed populations (see above), separately (Figure 1).

Foundation of the admixed population H

For all scenarios (Figure 1 However, MetHis does not allow simulating the source populations for the admixture process.

Instead, this can be done efficiently using coalescent-based simulations with existing software such as fastsimcoal2 [START_REF] Excoffier | fastsimcoal: a continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and SNP data[END_REF], or other forward-intime genetic data simulators such as SLIM v3 [START_REF] Haller | SLiM 3: Forward Genetic Simulations Beyond the Wright-Fisher Model[END_REF].

Simulating source populations

Here, we wanted to focus our investigation specifically on the admixture process undergone by the admixed population descending from the TAST. Therefore, we made several ad hoc simplification choices for simulating source population genetic data under the nine competing models described next.

We consider that the African and European populations at the source of the admixture processes are very large populations at the drift-mutation equilibrium, accurately represented by the Yoruban YRI and British GBR datasets here investigated. Therefore, we first build two separate datasets each comprising 20,000 haploid genomes of 100,000 independent SNPs, each SNP being randomly drawn in the site frequency spectrum (SFS) observed for the YRI and GBR datasets respectively. These two datasets are used as fixed gamete reservoirs for the African and European source population datasets separately, at each generation of the forward-in-time admixture process. From these reservoirs, at each generation separately, we build an effective individual gene-pool of diploid size Ng (see below), by randomly pairing gametes avoiding selfing. These virtual source populations provide the parental pool for simulating individuals in the admixed population H, at each generation separately. Thus, while our gamete reservoirs are fixed over the 21 generations of the admixture processes here considered, the parental genetic pools are randomly built anew at each generation of the admixture process.

Simulating the admixed population

At each generation, MetHis performs simple Wright-Fisher (FISHER 1922;[START_REF] Wright | Evolution in Mendelian Populations[END_REF] forward-in-time simulations, individual-centered, in a panmictic admixed population H of diploid effective size Ng. For a given individual in the hybrid population at the following generation (g + 1), MetHis independently draws each parent from the source populations with probability 𝑠 𝑆,𝑔 (Figure 1, Table 1), or from the hybrid population with probability ℎ 𝑔 , randomly builds a haploid gamete of 100,000 independent SNPs for each parent, and pairs the two constructed gametes to create the new individual. Here, we decided to neglect mutation over the 21 generations of admixture considered. This is reasonable when studying relatively recent admixture histories. Nevertheless, this will be improved in future versions of the software, in particular to allow studying much more ancient admixture histories. Finally, while we chose explicitly to simulate only the individuals in the admixed population H here, note that future developments of MetHis will allow to also simulate individual genetic data in the source populations in the same way.

Effective population size in the source and the admixed populations

To focus on the admixture process itself without excessively increasing the parameter space, we consider, for each nine-competing model, both source populations and the admixed population H with constant effective population size Ng = 1000 diploid individuals at each generation. Nevertheless, note that MetHis software readily allows the user to easily parameterize changes in the effective size of population H at each generation.

Sampling simulated unrelated individuals

After each simulation, we randomly draw individual samples matching sample sizes in our observed dataset: 90 and 89 individuals respectively from the African and European sources, and 50 individuals in the admixed population H. We sample individuals until our sample set contains no individuals related at the 1 st degree cousin within each population and between the admixed population and either source populations, based on explicit parental flagging during the last 2 generations of the simulations.

Simulating by randomly drawing parameter values from prior distributions

With this implementation of MetHis, we performed 10,000 independent simulations under each nine competing scenarios described above and in Figure 1, drawing the corresponding modelparameters (pulse-times and associated admixture intensities, "steepness" of the recurring admixture-increases or decreases and associated initial and final admixture intensities), in priordistributions detailed in Table 1. Although the user can perform MetHis simulations with an external parameter list, we readily provide ad hoc scripts in MetHis, which allow to easily generate parameter lists for a large number of complex admixture scenarios set by the user.

For the best models identified using Random-Forest ABC model-choice approach [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF] for the ACB and ASW admixed populations respectively (see Results), we conducted an additional 90,000 independent simulations with the same parameter priors as in the 10,000 simulations already conducted. Thus, we considered 100,000 simulations for the best scenarios for the ACB and ASW respectively, to be used for ABC posterior parameter inference (see below).

Summary Statistics

We considered 24 summary statistics for ABC model-choice and posterior parameter inference, computed on each simulated dataset with MetHis. Four statistics were strictly withinpopulations; four statistics were strictly between-populations; and 16 statistics were specifically calculated to describe the distribution of admixture among individuals within the admixed population H. Indeed, previous theoretical works have shown that this distribution and all its moments carried signatures of the underlying complex historical process [START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF]. Numerous descriptive statistical approaches have been successfully developed to estimate admixture fractions from genetic data in admixed populations (e.g. [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF][START_REF] Patterson | Ancient admixture in human history[END_REF][START_REF] Pickrell | Inference of population splits and mixtures from genome-wide allele frequency data[END_REF]. However, most methods remain computationally costly when iterated for large to very large sets of simulated genetic data. Therefore, only a few previous ABC historical inference approaches have considered the distribution of admixture fraction as a summary statistics [START_REF] Buzbas | AABC: approximate approximate Bayesian computation for inference in population-genetic models[END_REF][START_REF] Buzbas | Inference on admixture fractions in a mechanistic model of recurrent admixture[END_REF], although some admixture-related statistics have been embedded in ABC software packages [START_REF] Cornuet | DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data[END_REF].

Distribution of admixture fractions as a summary statistic

We estimated individual admixture distribution based on allele-sharing-dissimilarity (ASD) [START_REF] Bowcock | High resolution of human evolutionary trees with polymorphic microsatellites[END_REF]) and multidimensional scaling (MDS) [START_REF] Paschou | PCA-correlated SNPs for structure identification in worldwide human populations[END_REF]PRICE et al. 2009). For each simulated dataset, we first calculated a pairwise inter-individual ASD matrix using asd software (https://github.com/szpiech/asd) on all pairs of sampled individuals and all 100,000 independent SNPs. Then we projected in two dimensions this pairwise ASD matrix with classical unsupervised metric MDS using the cmdscale function R (R DEVELOPMENT CORE TEAM 2017). We expect individuals in population H to be dispersed along an axis joining the centroids of the two proxy source populations on the two-dimensional MDS plot. We projected individuals orthogonally on this axis, and calculate individual's relative distance to each centroid. We considered this measure to be an estimate of individual average admixture level from either source population. Note that by doing so, some individuals might show "admixture fractions" higher than one, or lower than zero, as they might be projected on the other side of the centroid when being genetically close to 100% from one source population or the other.

Under an ABC framework, this is not a difficulty since this may happen also on the real data a priori, and our goal is to use summary statistics that mimic the observed ones. This individual admixture estimation method has been shown to be highly concordant with cluster membership fractions as estimated with ADMIXTURE [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF] S2).

We used the mean, mode, variance, skewness, kurtosis, minimum, maximum, and all 10%quantiles of the admixture distribution obtained this way in population H, as 16 separate summary statistics for further ABC inference.

Within population summary statistics

We calculated SNP by SNP heterozygosities (NEI 1978) using vcftools [START_REF] Danecek | The variant call format and VCFtools[END_REF], and considered the mean and variance of this quantity across SNPs in the admixed population as two separate summary statistics for ABC inference. Note that, these quantities are fixed for each source population, respectively, and thus uninformative in our case study, since source populations are simulated only once and used for all subsequent simulations under the nine competing models (see above).

In addition, as we computed the individual pairwise ASD matrix for calculating the distribution of admixture fraction (see above), we also considered the mean and variance of ASD values across pairs of individuals within the admixed population H, as two within-population summary statistics.

Between populations summary statistics

In addition to previous summary statistics, we considered multilocus pairwise FST [START_REF] Weir | Estimating F-Statistics for the Analysis of Population-Structure[END_REF] between population H and each source population respectively, calculated

Prior-checking

We evaluated, a priori, if the above simulation design and novel tools can simulate genetic data for which summary statistics are coherent with those observed for the ACB and ASW as the targeted admixed population. To do so, we first plotted each prior summary statistics distributions and visually verified that the observed summary statistics for the ACB and ASW respectively fell within the simulated distributions (Supplementary Figure S3). Second, we explored the first four PCA axes computed with the princomp function in R, based on the 24 summary statistics and all 90,000 total simulations preformed for the nine competing scenarios, and visually checked that observed summary statistics were within the cloud of simulated statistics (Supplementary Figure S4). Finally, we performed a goodness-of-fit approach using the gfit function from the abc package in R, with 1,000 replicates and tolerance level set to 0.01 (Supplementary Figure S5).

Model-choice with Random-Forest Approximate Bayesian Computation

We used Random-Forest ABC (RF-ABC) for model-choice implemented in the abcrf function of the abcrf R package to obtain the cross-validation table and associated prior error rate using an out-of-bag approach (Figure 2). We considered the same prior probability for the nine competing models each represented by 10,000 simulations in the reference table. For the ACB and ASW observed data separately, we performed model-choice prediction and estimation of posterior probabilities of the winning model using the predict.abcrf function in the same R package, using the complete simulated reference table for training the Random-Forest algorithm (Figure 3, Supplementary Table S1). Both sets of analyses were performed considering 1,000 decision trees in the forest after visually checking that error-rates converged appropriately (Supplementary Figure S6), using the err.abcrf function in the R package abcrf.

Each summary statistics relative importance to the model-choice cross-validation was computed using the abcrf function (Figure 2). RF-ABC cross-validation procedures using groups of scenarios were conducted using the group definition option in the abcrf function [START_REF] Estoup | Model choice using Approximate Bayesian Computation and Random Forests: analyses based on model grouping to make inferences about the genetic history of Pygmy human populations[END_REF].

Posterior parameter estimation with Neural-Network Approximate Bayesian Computation

It Neural-Network tolerance level and number of neurons in the hidden layer

For each parameter estimation analysis, we determined empirically the NN tolerance level (i.e.

the number of simulations to be included in the NN training), and number of neurons in the hidden layer. Indeed, while the NN needs a substantial amount of simulations for training, there is also a risk of overfitting posterior parameter estimations when considering too large a number of neurons in the hidden layer. However, there are no absolute rules for choosing both numbers [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF][START_REF] Jay | An ABC Method for Whole-Genome Sequence Data: Inferring Paleolithic and Neolithic Human Expansions[END_REF].

Therefore, using the 100,000 simulations for the winning scenarios identified with RF-ABC (see above), we tested four different tolerance levels to train the NN (0.01, 0.05, 0.1, and 0.2), and a number of neurons ranging between four and seven (the number of free parameters in the winning scenarios, see Results). For each pair of tolerance level and number of neurons values, we conducted cross-validation checking of posterior parameter estimations with 1,000 randomly chosen simulated datasets in turn used as pseudo-observed data with the "cv4abc" function of the R package abc. We considered the median point-estimate of each posterior parameter (𝜃 ^𝑖) to be compared with the true parameter value used for simulation (𝜃 𝑖 ). The cross-validation parameter prediction error was then calculated across the 1,000 separate posterior estimations for pseudo-observed datasets for each pair of tolerance level and number of neurons, and for each parameter 𝜃 𝑖 , as ∑ (𝜃 ^𝑖 -𝜃 𝑖 ) 2 1000 1 (1000 × 𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝜃 𝑖 )) ⁄

, allowing

to compare errors for scenarios-parameters across NN tolerance-levels and numbers of hidden neurons, using the summary.cv4abc function in the R package abc [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF].

Results showed that, a priori, all numbers of neurons considered performed very similarly for a given tolerance level (Supplementary Table S2). Furthermore, results showed that considering 1% closest simulations to the pseudo-observed ones, to train the NN for parameter estimation, reduces the average error for each tested number of neurons. Thus, we decided to opt for four neurons in the hidden layer and a 1% tolerance level for training the NN in all subsequent NN-ABC analyses, in order to avoid overfitting in parameter posterior estimations.

Estimation of model-parameters posterior distributions for ACB and ASW

We jointly estimated model-parameters posterior distributions for the ACB and ASW admixed population separately, using 100,000 simulations for the best scenarios identified for each admixed population separately, using NN-ABC ("neuralnet" methods' option in the R package abc) based on the logit-transformed ("logit" transformation option in the R package abc)

summary statistics using a 1% tolerance level to train the NN (i.e. considering only the 1,000 closest simulations to the observed data), fitted using a single-hidden-layer neural network with four hidden neurons (Figure 4, Table 2).

Posterior parameter estimation error and credibility interval accuracy

For the ACB and ASW admixed populations separately, we wanted to evaluate the posterior error performed by our NN-ABC approach on the median point estimate of each parameter, in the vicinity of our observed data rather than randomly on the entire parameter space. To do so, we first identified the 1,000 simulations closest to the real data with a tolerance level of 1%, for the ACB and ASW respectively. Then, separately for the ACB and ASW set of closest simulations, we performed, similarly as above for the real data parameter estimation procedure, 1,000 separate NN-ABC parameter estimations using the "neural" method in the abc function with a NN trained with 1% tolerance level and four neurons in the hidden layer, using in turn the other 99,999 simulations as reference table, and recorded the median point estimate for each parameter. We then compared these estimates with the true parameter used for each 1,000 pseudo-observed target in the vicinity of our observed data and provide three types of error measurements in Finally, based on these cross-validation procedures, we evaluated a posteriori if, in the vicinity of the ACB and ASW observed datasets respectively, the lengths of the estimated 95% credibility intervals for each parameter was accurately estimated or not [START_REF] Jay | An ABC Method for Whole-Genome Sequence Data: Inferring Paleolithic and Neolithic Human Expansions[END_REF]. To do so, we calculated how many times the true parameter (𝜃 𝑖 ) was found inside the estimated 95% credibility interval [2.5%quantile(𝜃 𝑖 ^) ; 97.5%quantile(𝜃 𝑖 ^)], among the 1,000 out-of-bag NN-ABC posterior parameter estimation, separately for the ACB and ASW (Supplementary Table S3). For each parameter, if less than 95% of the true parameter values are found inside the 95% credibility interval estimated for the observed data, we consider the length of this credibility interval as underestimated indicative of a non-conservative behavior of the parameter estimation. Alternatively, if more than 95% of the true parameter-values are found inside the estimated 95% credibility interval, we consider its length as overestimated, indicative of an excessively conservative behavior of this parameter estimation.

Comparing the accuracy of posterior parameters estimations using NN, RF, or Rejection ABC With the above procedure, we aimed at estimating the posterior parameter distributions jointly for all parameters, and their errors for the scenario most likely explaining observed genetic data for the ACB and ASW respectively. Nevertheless, NN-ABC and RF-ABC parameter inference procedures also allow estimating each parameter posterior distribution in turn and separately rather than jointly. This can further provide insights into how both ABC parameter inference approaches perform in the parameter space of the winning scenarios. To do so, we performed several out-of-bag cross-validation parameter estimation analyses for the ACB and ASW separately.

We compared four methods: NN estimation of the parameters taken jointly as a vector (similarly as in the above procedure), NN estimation of the parameters taken in turn separately, RF estimation of the parameters which also considers parameters in turn and separately [START_REF] Raynal | ABC random forests for Bayesian parameter inference[END_REF], and simple Rejection estimation for each parameter separately [START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF]. For each method, we used in turn the 1,000 simulations closest to the real data as pseudoobserved data, and set a tolerance level of 1% of the 99,999 remaining simulations. We consider four neurons in the hidden-layer per neural network, and we considered 500 decision trees per random forest to limit the computational cost of these analyses at little accuracy cost a priori (Supplementary Figure S6). We then computed the mean-squared errors scaled by the variance of the true parameters, the mean-squared errors, and the mean absolute errors similarly as previously. Finally, we estimated the accuracy of the 95% credibility intervals for each method and for each parameter similarly as previously.

RESULTS

First, we present results about the ability of MetHis to simulate data close to the observed ones.

Second, we evaluate the ability of RF-ABC to distinguish, in the entire parameter space, the nine complex admixture scenarios in competition, and evaluate how each one of the 24 summary statistics contribute to distinguish among scenarios. Third, we use Random-Forest ABC to specifically predict the best fitting scenario for the history of admixture of two recently admixed populations descending from the Transatlantic Slave Trade in the Americas (African American ASW and Barbadian ACB). Fourth, we use Neural-Network ABC to estimate posterior parameter distributions under the winning scenario for the ACB and the ASW separately. Fifth, we evaluate in detail the accuracy of our posterior parameter estimation, and compare with other ABC posterior parameter inference approaches. Finally, we synthesize the complex admixture history thus reconstructed for the ASW and ACB populations.

Simulating the observed data with MetHis

With MetHis, we conducted 10,000 simulations for each one of the nine competing scenarios for the admixture history of the ASW or the ACB populations, described in detail in Figure 1 and Material and Methods, with corresponding model parameters drawn in a priori distributions described in Table 1.

We produced 90,000 vectors of 24 summary statistics each, overall highly consistent with the observed ones for the ACB and the ASW populations respectively. First, we found that each observed statistic is visually reasonably well simulated under the nine competing scenarios here considered (Supplementary Figure S3). Second, the observed data each fell into the simulated sets of summary statistics projected in the first four PCA dimensions (Supplementary Figure S4) considering all 24 summary statistics in the analysis. Finally, the observed vectors of 24 summary statistics computed for the ACB and ASW, respectively, were not significantly different (p-value = 0.468 and 0.710 respectively) from the 90,000 simulated sets of statistics using a goodness-of-fit approach (Supplementary Figure S5). Therefore, we successfully simulated datasets producing sets of summary statistics reasonably close to the observed ones, despite considering constant effective population sizes, fixed virtual source population genetic pool-sets, and neglecting mutation during the 21 generations of forward-in-time simulations performed using MetHis.

Complex admixture scenarios cross-validation with RF-ABC

We trained the RF-ABC model-choice algorithm using 1,000 trees, which guaranteed the convergence of the model-choice prior error rates (Supplementary Figure S6). Based on this training, the complete out-of-bag cross-validation matrix showed that the nine competing scenarios of complex historical admixture could be relatively reasonably distinguished using our set of 24 summary statistics and 10,000 simulations under each competing scenario, despite the high level of nestedness of the scenarios here considered (see Material and Methods).

Indeed, we calculated an out-of-bag prior error rate of 32.41%, considering each 90,000 simulation, in turn, as out-of-bag pseudo-observed target dataset and the rest of simulations (89,999) as the training dataset for RF-ABC scenario-choice. Furthermore, we found the posterior probabilities of identifying the correct scenario ranging from 55.17% (prior probability = 11.11% for each competing scenario), for the two-pulses scenarios from both the African and European sources (Afr2P-Eur2P), to 77.71% for the scenarios considering monotonically decreasing recurring admixture from both sources (AfrDE-EurDE) (Figure 2A).

Importantly, the average probability, for a given admixture scenario, of choosing any one alternative (wrong) scenario were on average 4.05% across the eight alternative scenarios, ranging from 2.79% for the AfrDE-EurDE scenario, to 5.60% for the Afr2P-Eur2P scenario (Figure 2A). This shows that our approach did not systematically favor one or the other competing scenario when wrongly choosing a scenario instead of the true one, despite high levels of nestedness among scenarios.

We find that the six summary statistics most contributing to the observed cross validation results

for RF-ABC model-choice among the 24 statistics here tested were statistics describing specifically the admixture-fraction distribution: minimum and maximum admixture fraction values, variance, skewness, as well as the 10% and 90% quantiles of the distribution (Figure 2B). Interestingly, within and between populations summary-statistics often used in population genetics (including FST, mean heterozygosity, and f3 statistics), contributed to distinguishing the competing complex admixture scenarios to a lesser extent.

Finally, note that scenarios considering monotonically recurring admixture from each source populations (AfrDE-EurDE, AfrDE-EurIN, AfrIN-EurDE, AfrIN-EurIN) can be relatively well distinguished, using our RF-ABC framework, from scenarios with at least one source population contributing to the admixed population with two possible pulses after the foundation event (Afr2P-Eur2P, Afr2P-EurDE, Afr2P-EurIN, AfrDE-Eur2P, AfrIN-Eur2P). Indeed, we found an out-of-bag prior error rate of 13.85%, and posterior cross-validation probabilities of identifying the correct group of scenarios of 86.08% and 86.23% respectively for the two groups [START_REF] Estoup | Model choice using Approximate Bayesian Computation and Random Forests: analyses based on model grouping to make inferences about the genetic history of Pygmy human populations[END_REF].

Complex admixture histories for the Barbadian and African American populations

Random-Forest ABC scenario-choice

We performed RF-ABC model-choice with 1,000 decision trees and 10,000 simulations per each nine competing scenarios (Figure 1 and Table 1, Material and Methods), separately for the admixture history of the Barbadian (ACB) and the African American (ASW) populations.

For the ACB, Figure 3 shows that the majority of votes (53.1%) went to an admixture scenario AfrDE-EurDE with a posterior probability of the winning scenario of 60.28%. This scenario encompassed monotonically decreasing recurring contributions from both the African and European source populations over the last 20 generations before present. The second most chosen scenario considered a monotonically decreasing recurring contribution from the African source population over the last 20 generations, while the European source population contributed two admixture pulses to this admixed population after the founding pulse (scenario AfrDE-Eur2P). However, this scenario is voted for 3.5 times less often than the winning scenario AfrDE-EurDE, gathering 15.1% of the 1,000 votes, only slightly above the 11.11%

prior probability for each nine-competing scenario (Figure 3; Supplementary Table S1).

Concerning the admixture history of the ASW, RF-ABC scenario-choice results were less segregating. Figure 3 shows that the AfrDE-EurDE scenario also gathered the majority of votes for the admixture history of the ASW, albeit with lower posterior probability than for the ACB (33.5% of 1,000 votes, with posterior probability = 48.0% for the ASW). The second most chosen scenario, AfrDE-Eur2P, was only slightly less chosen with 31.7% of the votes (Figure 3, Supplementary Table S1). For the ASW, considering only the two best scenarios (AfrDE-EurDE and AfrDE-Eur2P) to train the Random Forest, and re-conducting the RF-ABC scenario-choice, improved the scenario discrimination in favor of the AfrDE-EurDE scenario.

While we found only a slight majority of votes (51.8%) also in favor of the AfrDE-EurDE scenario, we found a substantially increased posterior probability for this model equal to 57.9%.

This increased posterior probability of the AfrDE-EurDE scenario compared to the previous RF-ABC scenario-choice considering the nine competing scenarios (48.0%), indicated that this scenario best explains the ASW observed genetic patterns, despite overall limited discriminatory power of our approach in the part of the summary-statistics space occupied by the ASW.

Neural-Network ABC parameter inference accuracy for the ACB and ASW populations

We performed 100,000 simulations using MetHis for the AfrDE-EurDE scenarios, in order to estimate, using Neural-Network ABC, posterior parameter distributions and the corresponding parameter prediction cross-validation errors, considering in turn the ACB and the ASW populations (Figure 4 and 5, Table 2, Table 3, andSupplementary Table S3).

For the ACB under the AfrDE-EurDE scenario (Figure 4A, S3). This shows the reliability of our method to accurately infer the three parameters in the part of the space of summary statistics occupied by the ACB observed data.

Furthermore, the two ancient admixture intensities from Africa and Europe at generation 1 immediately following the initial foundation of the admixed population H (sAfr,1 and sEur,1, respectively), also had posterior densities apparently distinguished from their prior distributions, but both had much wider 95% credibility intervals (Figure 4A, Table 2).

Consistently, we found a slightly increased posterior parameter error in this part of the parameter space for both these parameters, with average absolute error 0.121 and 0.095 respectively for sAfr,1 and sEur,1 (Table 3). Nevertheless, note that 95.8% and 94.7% of 1,000 cross-validation true values for those two parameters fell into the estimated 95% credibility intervals (Supplementary Table S3). This shows a reasonably conservative behavior of our method for these estimations, further indicating that information is lacking in our data or set of summary statistics for a more accurate estimation of these parameters, rather than an inaccuracy of our approach.

Interestingly (Figure 4A, Table 2), we found that accurate posterior estimation of the steepness of the African decrease in admixture over time (uAfr) is difficult. Indeed, the posterior density of this parameter only showed a tendency towards small values slightly departing from the prior, indicative of a limit of our method to estimate this parameter (Figure 4A, Table 2).

Finally (Figure 4A, Table 2), we found that we had virtually no information to estimate the founding admixture proportions from Africa and Europe at generation 0, as our posterior estimates barely departed from the prior and associated mean absolute error was high (0.2530, Table 3). Nevertheless, our method seemed to be performing reasonably conservatively for these two latter parameters (95.6% and 95.3% of 1,000 cross-validation true parameter values fell into estimated 95% credibility intervals, Supplementary Table S3). This indicates that information is strongly lacking in our data or summary statistics for successfully capturing these parameters, rather than inherent inaccuracy of our ABC method.

For the African American ASW under the AfrDE-EurDE model, our posterior parameter estimation accuracy results were overall quantitatively slightly less accurately estimated compared to those obtained for the ACB population, as indicated by overall larger credibility intervals and cross-validation errors (Figure 4B, Table 2, Table 3, Supplementary Table S3).

This was consistent with the more ambiguous RF-ABC model-choice results obtained for this population (Figure 3).

Comparing NN, RF, and Rejection ABC posterior parameter estimation accuracy

For posterior parameter estimations considering the ACB or the ASW population, the means of the three types of errors (scaled mean-square error, mean-square error, absolute error, see Material and Methods) were systematically lower for the two NN methods (joint or independent posterior parameter estimation) than for the RF and Rejection independent posterior parameter estimation methods (Table 4). Furthermore, we found that the means of the three types of errors were qualitatively comparable between the NN estimation of the parameters taken as a joint vector and the NN estimation of the parameters taken separately.

Altogether, these results showed that considering the NN estimation for parameters taken jointly as a vector is overall preferable for the ACB and ASW populations, since it further allowed the joint interpretation of parameter values estimated a posteriori, with little difference in accuracy between the two methods.

Finally, results showed that the lengths of 95% credibility intervals estimated with NN joint parameter estimation was, across all parameters, more accurate than all other methods with, on average, 95.1% and 95.2% of true parameter values falling within the estimated 95% credibility intervals, for the ACB and ASW respectively (Supplementary Table S3). Furthermore, we found that lengths of 95% credibility intervals estimated with NN and RF independent posterior parameter estimations were systematically under-estimated, with less than 94% of the true parameter values falling into the 95% credibility intervals estimated. Finally, we found that lengths of 95% credibility intervals estimated with the Rejection method were also rather accurately estimated although on average slightly over-estimated compared to the NN joint parameter estimation with, on average, 95.5% of the 1,000 cross-validation true parameter values within the estimated 95% credibility intervals for the ACB, and 95.8% for the ASW.

Admixture histories of the African American ASW and Barbadian ACB Figure 5 visually synthesized the estimated posterior parameters of the complex admixture scenarios reconstructed with our novel MetHismachine-learning ABC framework, and associated 95% credibility intervals (Table 2).

We found a virtual complete replacement of the ACB and ASW populations at generation 1 after foundation, thus consistent with our inability to accurately estimate the founding proportions from the African and European sources at generation 0. Furthermore, we found an increasingly precise posterior estimation of African and European contributions to the genepool of the ACB and ASW populations forward in time, with most recent estimations exhibiting narrow credibility intervals. This is also consistent with the nature of recurrent admixture processes, where older information is often lost or replaced when more recent admixture events occur.

Most interestingly, we found that the recurring contribution of the European gene pool to the admixed populations rapidly decreases after generation 1 for both the ACB and ASW albeit with substantial differences (Figure 5). Indeed, we found that the recurring contribution from the European source to the ACB gene pool falls below 10% at generation 9 until no more than 1% in the present (generation 20). Comparatively, we found that the European contribution diminished substantially less rapidly for the ASW, going below 10% only after generation 12 until roughly 2% in the present. This indicates that the European contribution to the African American gene pool was more sustained over time than for the Barbadian.

Finally, we found substantial recurring contributions from the African source population to the gene pool of both admixed populations (Figure 5). For the ACB population, we found a progressive decrease of the African recurring introgression until a virtually constant recurring admixture close to 28% from generation 10 and onward. For the ASW, our results showed a sharper decrease of the African contribution after foundation until a virtually constant recurring admixture process close to 24% from generation 5 until present (generation 20). The high overall African recurring introgression into the admixed-populations gene pools captures the importance of recurring admixture in explaining the observed patterns for both populations 767 descending from the TAST.

DISCUSSION

We evaluated how machine-learning Approximate Bayesian Computation methods can bring new insights to the reconstruction of highly complex admixture histories using genetic data. To illustrate our proof of concept and thoroughly investigate the power and accuracy of our approach using real data, we aimed at reconstructing the recent complex admixture history for the African American (ASW) and Barbadian (ACB) population samples from the 1000 Genomes project (Phase 3).

Our results demonstrated that our novel MetHis forward-in-time simulator and summary statistics calculator coupled with RF-ABC scenario-choice can often clearly infer the best class of highly complex admixture histories underlying independent SNP data diversity, in a reasonable-size sample and genetic dataset. In the two source-populations admixture models here investigated, we distinguished scenarios encompassing two pulses of admixture from each source, after the founding admixture event, monotonically increasing or decreasing admixture intensities over time, or a combination of these three scenarios. Furthermore, we found that NN-ABC provide accurate posterior parameter inference of most demographic parameters of recurring monotonically decreasing admixture processes, compared to other classes of ABC posterior parameter inference methods. Finally, we empirically demonstrated that the moments of the distribution of admixture fractions within the admixed population estimated using independent SNPs were highly informative for reconstructing the admixture history using an ABC approach, as expected theoretically [START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF]).

While we found that distinguishing among competing models is more difficult in certain parts of the parameter space due to scenario-nestedness [START_REF] Robert | Model choice versus model criticism[END_REF], our MetHis -ABC method already vastly extends the array of complex admixture models explored with most, classically used, maximum-likelihood inference approaches [START_REF] Robert | Model choice versus model criticism[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Loh | Inferring admixture histories of human populations using linkage disequilibrium[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF]. It is challenging to analytically predict genomic diversity patterns expected under realistic complex admixture histories, as likelihood calculations under such models are very often intractable [START_REF] Verdu | A general mechanistic model for admixture histories of hybrid populations[END_REF][START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Medina | Estimating the Timing of Multiple Admixture Pulses During Local Ancestry Inference[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. In turn, this makes it difficult to understand how most existing efficient maximum-likelihood admixture inference methods, which often only consider one or two pulses of admixture, behave when the observed genetic data in fact results from much more complex admixture processes (GRAVEL 2012; HELLENTHAL et al.

2014).

In this context, the proof of concept here presented more generally shows that ABC can be fruitfully attempted to explore, virtually, any other admixture model beyond the case studies here-conducted, provided that, a priori, simulation and summary statistics calculation are feasible. To these ends, other recent efficient forward-in-time genetic data simulators can also be successfully used in an ABC framework instead of MetHis [START_REF] Haller | SLiM 3: Forward Genetic Simulations Beyond the Wright-Fisher Model[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. In reality, studies investigating ABC approaches for admixture reconstruction, while allowing for exploring scenarios out of reach of other methods, will inevitably face the same difficulties as any ABC inference; such as high dimensional parameter and summarystatistics spaces, lack of information from summary statistics, and scenario nestedness [START_REF] Csilléry | Approximate Bayesian Computation (ABC) in practice[END_REF][START_REF] Robert | Model choice versus model criticism[END_REF][START_REF] Sisson | Handbook of Approximate Bayesian Computation[END_REF].

Importantly, the current MetHis -ABC approach does not make use of admixture linkagedisequilibrium patterns in the admixed population, and only relies on independent genetic markers. Nevertheless, admixture LD has consistently proved to bring massive information about the complex admixture history of numerous populations worldwide [START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Hellenthal | A genetic atlas of human admixture history[END_REF][START_REF] Medina | Estimating the Timing of Multiple Admixture Pulses During Local Ancestry Inference[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. However, existing methods to calculate admixture LD patterns remain computationally intensive and require numerous markers and accurate phasing, which is difficult under ABC where such statistics have to be calculated for each one of the numerous simulated datasets. In this context, RF-ABC [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF][START_REF] Raynal | ABC random forests for Bayesian parameter inference[END_REF] or AABC [START_REF] Buzbas | AABC: approximate approximate Bayesian computation for inference in population-genetic models[END_REF] methods allow substantially diminishing the number of simulations required for satisfactory scenario-choice and posterior parameter inference, which makes both approaches promising tools for using, in the future, admixture LD patterns to reconstruct complex admixture processes from genomic data.

Sex-biased admixture processes are known to have influenced admixed populations, and in particular populations descending from the TAST [START_REF] Moreno-Estrada | Reconstructing the population genetic history of the Caribbean[END_REF][START_REF] Fortes-Lima | Exploring Cuba's population structure and demographic history using genome-wide data[END_REF]. Future version of our MetHis -ABC framework will explicitly implement sex-specific admixture processes with, in addition to autosomal data, the possibility to investigate sex-related genetic data (X-chromosome, Y-chromosome, and mitochondrial DNA) [START_REF] Goldberg | Autosomal admixture levels are informative about sex bias in admixed populations[END_REF][START_REF] Goldberg | Beyond 2/3 and 1/3: The Complex Signatures of Sex-Biased Admixture on the X Chromosome[END_REF].

Finally, although MetHis readily allows considering changes of effective population size in the admixed population at each generation as a parameter of interest to ABC inference, we did not, for simplicity, investigate here how such changes affected our results for the African American and Barbadian admixed population. Future work using MetHis will allow specifically investigating how effective size changes may influence genetic patterns in the admixed population, a question of major interest as numerous admixed populations are expected to have experienced founding events and/or bottlenecks during their history (e.g. BROWNING et al.

2018).

For all these reasons, it is crucial, in general and in the future, to further develop novel methodological tools and evaluate how genetic patterns evolve over time as a function of each parameter of complex historical admixture models separately [START_REF] Buzbas | Inference on admixture fractions in a mechanistic model of recurrent admixture[END_REF][START_REF] Medina | Estimating the Timing of Multiple Admixture Pulses During Local Ancestry Inference[END_REF][START_REF] Ni | MultiWaver 2.0: modeling discrete and continuous gene flow to reconstruct complex population admixtures[END_REF]. MetHis can help to this task since it allows the users to investigate how parameters of the complex admixture process can influence, over time, a large number of population genetics summary-statistics calculated in the simulated admixed population at each generation.

Concerning the specific admixture history of the two admixed populations descending from the TAST here reconstructed, note that several competing scenarios can clearly be discarded for explaining the observed genetic patterns. In particular, the Afr2P -Eur2P scenario considering two possible pulses of introgression after the founding event, separately from the African and European source, does not significantly exceeds the prior probability of choosing any ninecompeting scenario (4.6% and 11.2% of the 1,000 votes for the Afr2P -Eur2P, respectively for the ACB and the ASW, Figure 3). Note that this scenario embeds models analogous to the most complex admixture scenarios that have been previously tested for these populations with maximum-likelihood approaches based on extensive genome-wide data and admixture-LD based statistics [START_REF] Gravel | Population genetics models of local ancestry[END_REF][START_REF] Baharian | The Great Migration and African-American Genomic Diversity[END_REF]. Interestingly, very recent migrations from either Africa or Europe to the Americas are known to have been intense demographically in the 19 th and 20 th century (BERLIN 2010). However, the recent increased demographic migrations do not seem to have left the equivalent signature in the genetic admixture process of both the ACB and ASW populations, as monotonically recurring increasing admixture scenarios can here be rejected confidently.

Nevertheless, we found that genetic admixture of African origin in both admixed populations, although decreasing since foundation, retained high levels in the present day (between 20% and 30%). These results could stem from the known importance of African recurring forced migrations during the TAST into the Americas; further prompts the influence of African slave descendants forced migrations within the Americas after the initial crossing of the ocean (often called the Middle Passage); and highlights the major importance of post-slavery migrations of TAST descendant populations within the Americas [START_REF] Berlin | The making of African America : the four great migrations[END_REF][START_REF] Eltis | The Lewis Walpole series in eighteenth-century culture and history[END_REF][START_REF] Baharian | The Great Migration and African-American Genomic Diversity[END_REF]. For instance, intense migrations from Haitian slave-descendants in the 19 th century have already been shown to possibly have contributed to the admixture patterns of other populations in the Caribbean and continental America [START_REF] Moreno-Estrada | Reconstructing the population genetic history of the Caribbean[END_REF][START_REF] Fortes-Lima | Exploring Cuba's population structure and demographic history using genome-wide data[END_REF].

Finally, we found that the genetic contribution from Europe rapidly decreases, after the foundation of both admixed populations, to marginal amounts during the 20 th century.

Therefore, it seems that neither sustained European migrations, nor the relaxation of social and legal constraints on admixture between descendant communities subsequent to the abolition of slavery and the end of segregation, have translated into increased European genetic contribution to the gene-pool of admixed populations descending from European and African forced or voluntary migrations into the Americas after the TAST.

Altogether, our results for the two recently-admixed human populations illustrated how our MetHis -ABC framework can bring fundamental new insights into the complex demographic history of admixed populations; a framework that can easily be adapted for investigating admixture history in numerous populations and species, particularly when maximum-likelihood methods are intractable. Summary statistics' respective importance in the RF-ABC model-choice out-of-bag crossvalidation [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF]. 1). 1,000 decision trees were considered in the model-choice prediction, respectively for each population.

Supplementary note S1

We used the rectangular hyperbola class of functions to obtain increasing/decreasing patterns using only one shape parameter. We give here the derivation of the equations used, giving the example of a decreasing pattern.

A decreasing hyperbola is given by the function: The intersection between the hyperbola and y = x is given by x = y = -a + a + a 2 thus, we can sample an uniform deviate u ∈ [0; 1 2 ] and set parameter a:

f (x) = a(1 -x) a + x (S1.1) with x ∈ [0; 1], f (
a = u 2 1 -2u
to obtain all hyperbola shapes.

We then transformed equation S1.1 to rescale the ranges of x and f (x): With the notation used in the main text for contributions, and considering 20 generations of admixture, we obtain:

f (x) = a(y max -y min )(1 - x -x min x max -x min ) a + x -x min x max -x min + y min (S1 .
s S,g = a(s S,1 -s S,20 )(1 - g -1 20 -1 ) a + g -1 20 -1 + s S,20 (

  the alpha version of the software -Conducted benchmarking and data analyses -Helped writing the article RL: Built the beta version of the software -Conducted benchmarking and data analyses -Helped writing the article VT: Conducted benchmarking and data analyses -Helped writing the article BT: Helped building the beta version of the software -Conducted benchmarking and data analyses -Helped writing the article PV: Designed and supervised the project -Conducted benchmarking and data analyses -Wrote the article

  , we propose a novel ad hoc forward-in-time genetic data simulator and a set of parameter-generator and summarystatistics calculation tools embedded in an open source C software package called MetHis. It is adapted to conduct primarily ABC inferences with existing ABC tools implemented in the R (R DEVELOPMENT CORE TEAM 2017) packages abc (CSILLÉRY et al. 2012) and abcrf (PUDLO et al. 2016; RAYNAL et al. 2019).

Figure 3 :

 3 Figure 3: Random-Forest Approximate Bayesian Computation model-choice predictions for the ACB (left panel) and ASW (right panel) populations. Nine competing models were compared, each with 10,000 simulations (Figure 1, Table1). 1,000 decision trees were

Figure 4 :

 4 Figure 4: Neural-Network Approximate Bayesian Computation posterior parameters estimated densities under the winning scenario AfrDE-EurDE, for (A) the ACB and (B) the ASW populations. Median posterior point estimates are indicated by the red vertical line, 95% credibility intervals are indicated by the colored area under the posterior curve (Table2). All

Figure

  Figure S1.1: Influence of a on equation S1.1

Figure S1. 2 :

 2 Figure S1.2: Influence of a on equation S1.2 with x min = 1, x max = 20, y min = 0.1 and y max = 0.9

  Figure 2

, Table 1) we

	years per generation and the sampled generation born in the 1980s. Note that simulations
	considering a parameter sAfr,0 close to 0, or alternatively 1, correspond to foundations of the
	population H from either one source population, therefore delaying the first "real" genetic
	admixture event to the next, more recent, demographic event. Following foundation, we
	consider three alternative scenarios for the admixture contribution of each source population S,
	African or European in our case, separately.
	Admixture-pulse(s) scenarios
	For a given source population S, African (Afr) or European (Eur), scenarios S-2P consider two
	possible pulses of admixture into population H occurring respectively at time tS,p1 and tS,p2
	distributed in [1,20] with tS,p1 ≠ tS,p2, with associated admixture proportion sS,tS,p1 and sS,tS,p2 in
	[0,1] satisfying, at all times t, ∑	𝑆∈(𝐴𝑓𝑟,𝐸𝑢𝑟)	𝑠 𝑆,𝑡	≤1 (Figure 1, Table 1). Note that for one of
	either sS,t parameter values close to 0, the two-pulse scenarios are equivalent to single pulse
	scenarios after the foundation of H. Furthermore, for both sS,t values close to 0, scenarios S-2P
	are nested with scenarios where only the founding admixture pulse 21 generations ago is the
	source of genetic admixture in population H. Alternatively, sS,t parameter values close to 1
	consider a virtual complete genetic replacement of population H by source population S at that
	time. Finally, certain S-2P scenarios with two consecutive pulses from a given source S (tS,p1 =
	tS,p2 -1), may be strongly resembling single-pulse scenarios (after foundation).
	Recurring decreasing admixture scenarios
	For a given source population S, scenarios S-DE consider a recurring monotonically decreasing
	admixture from source population S at each generation between generation 1 (after foundation
	at generation 0) and generation 20 (sampled population) (
		chose a fixed time for the foundation (generation 0,
	forward-in-time) of population H occurring 21 generations before present, with admixture
	proportions sAfr,0 and sEur,0 from the African and the European sources respectively, with sAfr,0
	+ sEur,0 = 1, and sAfr,0 in [0,1]. This corresponds to the first arrival of European permanent settlers
	in the Americas and Caribbean in the late 15 th and early 16 th centuries, considering 20 or 25

Figure 1, Table 1). In these scenario,

	immediately after foundation of intensity sS,1, followed by constant recurring admixture at each
	generation until present of intensity sS,20. Alternatively, parameter uS values close to 1/2 create
	scenarios with a linearly decreasing admixture between sS,1 and sS,20 from source population S
	at each generation after the foundation of population H.
	Recurring increasing admixture scenarios
	Finally, for a given source population S, scenarios S-IN mirrors the S-DE scenarios by
	considering instead a recurring monotonically increasing admixture from source population S
	(Figure 1, Table 1). Here, sS,g, with g in [1..20], are the discrete numerical solutions of the
	same function as in the S-DE decreasing scenarios (see above), flipped over time between
	generation 1 and 20. In these scenarios, sS,20 is defined in [0,1] and sS,1 in [0,sS,20/3], and u in
	[0,1/2] parametrizes the "steepness" of the curvature of the increase. Note that in this case,
	parameter u values close to 0 create pulse-like scenarios occurring in the present of intensity
	sS,20, preceded by constant recurring admixture of intensity sS,1 at each generation since
	foundation. Alternatively, parameter uS values close to 1/2 create scenarios with a linearly
	increasing admixture between sS,1 and sS,20 from source population S at each generation after
	the foundation of population H.
	Combining admixture scenarios from either source populations
	We combine these three scenarios to obtain nine alternative scenarios with two source
	populations, African (Afr) and European (Eur) respectively, for the admixture history of
	population H (
	sS,g, with g in [1..20], are the discrete numerical solutions of a rectangular hyperbola function
	over the 20 generations of the admixture process until present as described in Supplementary
	Note S1. In brief, this function is determined by parameter uS, the "steepness" of the curvature
	of the decrease, in [0,1/2], sS,1, the admixture proportion from source population S at generation
	1 (after foundation), in [0,1], and sS,20, the last admixture proportion in the present, in [0,sS,1/3].
	Note that we chose the boundaries for sS,20 in order to reduce the parameter space and nestedness
	among competing scenarios, and explicitly force scenarios S-DE into a substantially decreasing
	admixture process. Indeed, defining sS,20 in [0,sS,1] instead would have also allowed for both
	decreasing admixture processes and relatively constant recurring admixture processes.
	Furthermore, note that parameter uS values close to 0 create pulse-like scenarios occurring

Figure 1, Table 1), the ASW or ACB alternatively, with the only condition that,

	at each generation g in [1..20], parameters satisfy sAfr,g + sEur,g + hg = 1, with hg, in [0,1], being
	the remaining contribution of the admixed population H to itself at the generation g. Four
	scenarios (Afr2P-EurDE, Afr2P-EurIN, AfrDE-Eur2P, and AfrIN-Eur2P) consider a mixture
	of pulse-like and recurring admixture from each source. Three scenarios (Afr2P-Eur2P, AfrDE-
	EurDE, and AfrIN-EurIN), consider symmetrical classes of admixture scenarios from either
	source. Two scenarios (AfrIN-EurDE and AfrDE-EurIN) consider mirroring recurring
	admixture processes. Importantly, this scenario design considers nested historical scenarios in
	specific parts of the parameter space, as exemplified above.

Forward-in-time simulations with MetHis

	2012). In this context, we developed MetHis, a C open-source software package available at
	https://github.com/romain-laurent/MetHis, to simulate large amounts of genetic data under the
	two-source populations VR2011 model and calculate summary statistics of interest to
	population geneticists interested in complex admixture processes. MetHis, in its current form,
	can be used to simulate any number of independent SNPs in the admixed population H.
	Simulation of genome-wide independent SNPs under highly complex admixture histories is
	often not trivial under the coalescent and using classical existing software (WAKELEY et al.

  in real data analyses (e.g.[START_REF] Verdu | Parallel Trajectories of Genetic and Linguistic Admixture in a Genetically Admixed Creole Population[END_REF]. Considering the real data here investigated, we confirm these previous findings since we obtain a Spearman correlation (calculated using the cor.test function in R), of

	rho = 0.950 (p-value < 2.10 -16 ) and rho= 0.977 (p-value < 2.10 -16 ) between admixture estimates
	based on ASD-MDS and on ADMIXTURE, for the ACB and ASW respectively
	(Supplementary Figure

  Furthermore, we calculated the mean ASD between individuals in population H and, separately, individuals in either source population. Finally, we computed anew from[START_REF] Patterson | Ancient admixture in human history[END_REF]) the f3 statistics based on allelic frequencies obtained with vcftools[START_REF] Danecek | The variant call format and VCFtools[END_REF]). In the two-source population case, this statistic is extensively employed to test the original source of the admixture of a target admixed population, infer the time since admixture, and estimate admixture intensities using maximum-likelihood approaches.

	Approximate Bayesian Computations
	MetHis has been designed to operate under an ABC framework for model choice and parameter
	inference. Thus, it allows simulating genetic data under numerous possible models by drawing

using vcftools

[START_REF] Danecek | The variant call format and VCFtools[END_REF]

. Note that the FST between the source populations is fixed, since simulated source populations are themselves fixed (see above), and thus uninformative in our case study.

parameter values in a priori distributions set by the user in a flexible way. In addition, MetHis allows for calculating numerous summary statistics a priori of interest to admixture processes, and provides, as outputs, scenarios-parameter vectors and corresponding summary-statistics vectors in reference tables ready to be used with the machine-learning ABC abc

[START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF] abcrf (PUDLO et al. 2016;[START_REF] Raynal | ABC random forests for Bayesian parameter inference[END_REF]

) R packages (R DEVELOPMENT CORE TEAM 2017).

Table 3 .

 3 The mean-squared error scaled by the variance of the true parameter

	∑	(𝜃 ^𝑖 -𝜃 𝑖 ) 1000 1	2	(1000 × 𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝜃 𝑖 )) ⁄	as previously (Csilléry et al. 2012); the mean-
	squared error ∑	1000 1	(𝜃 ^𝑖 -𝜃 𝑖 )	2	1000 ⁄	, allowing to compare estimation errors for a given
	scenario-parameter between the ACB and ASW analyses; and the mean absolute error
	∑	|𝜃 ^𝑖 -𝜃 𝑖 | 1000 1	1000 ⁄	, which provides a more intuitive parameter estimation error.

Table 2 )

 2 , we found that the two recent admixture intensities fromAfrica and Europe (sAfr,20 and sEur,20, respectively) and the steepness of the European decrease in contribution over time (uEur) had sharp posterior densities clearly distinct from their respective priors. Note that the cross-validation error on these parameters in the vicinity of our real data were low (average absolute error 0.02744, 0.0044, and 0.1084, respectively for sAfr,20, sEur,20, and uEur) (Table3), and lengths of 95% credibility intervals reasonably accurate (96.4%, 94.4%, 94.1% of 1,000 cross-validation true parameter values fell into estimated 95% credibility intervals, Supplementary Table

Table 2

 2

	). All

Table 1 .

 1 Parameter prior distributions for simulation with MetHis and Approximate Bayesian Computations historical inference. Parameter list correspond to the nine competing historical admixture models described in Figure1and Material and Methods.

Table 2 .

 2 Neural-Network Approximate Bayesian Computation posterior parameter weighted distributions under the winning scenario AfrDE-EurDE, for the ACB and ASW populations.

	All posterior parameter estimations were conducted using 100,000 simulations under scenario
	AfrDE-EurDE (Figure 1, Table 1), a 1% tolerance rate (1,000 simulations), 24 summary
	statistics, logit transformation of all parameters, and 4 neurons in the hidden layer (see Material
	and Methods).

Table 3 .

 3 Neural-Network Approximate Bayesian Computation posterior parameter errors under

	the winning scenario AfrDE-EurDE, for the ACB and ASW populations. For each target
	population separately, we conducted cross-validation by considering in turn 1,000 separate NN-
	ABC parameter inferences each using in turn one of the 1,000 closest simulations to the
	observed ACB (or ASW) data as the target pseudo-observed simulation. All posterior parameter
	estimations were conducted using 100,000 simulations under scenario AfrDE-EurDE (Figure
	1, Table 1), a 1% tolerance rate (1,000 simulations), 24 summary statistics, logit transformation
	of all parameters, and four neurons in the hidden layer (see Material and Methods). Median
	was considered as the point posterior parameter estimation for all parameters. First column
	provides the average absolute error; second column shows the mean-squared error; third
	column shows the mean-squared error scaled by the parameter's observed variance (see
	Material and Methods for error formulas).

Table 4 .

 4 Approximate Bayesian Computation mean posterior parameter errors under the winning Scenario AfrDE-EurDE, for the ACB and ASW populations separately, using four different methods: NN estimation of the parameters taken jointly as a vector, NN estimation of the parameters taken separately, Random Forest (parameters taken separately), and Rejection the 1,000 simulations closest to our real data, separately for the ACB and ASW, as pseudoobserved datasets for four separate methods : NN estimation of the parameters taken jointly as a vector, NN estimation of the parameters taken independently, Random Forest (parameters are taken independently), and Rejection (parameters are taken independently).

	Table 1.

Table 2 .

 2

		AfrDE-EurDE
		parameters	Median	Mean	Mode	95% Credibility Interval
	ACB	sAfr,0	0.3097	0.3747	0.1121	[0.0116 ; 0.9347]
		sAfr,1	0.6797	0.6769	0.6813	[0.4577 ; 0.8880]
		sAfr,20	0.2707	0.2655	0.2788	[0.1985 ; 0.2967]
		uAfr	0.1409	0.1684	0.0508	[0.0041 ; 0.4507]
		sEur,1	0.1807	0.2160	0.1158	[0.0542 ; 0.5525]
		sEur,20	0.0100	0.0102	0.0093	[0.0018 ; 0.0200]
		uEur	0.4858	0.4627	0.4929	[0.1886 ; 0.4992]
	ASW	sAfr,0	0.5258	0.5124	0.7015	[0.0262 ; 0.9758]
		sAfr,1	0.6006	0.6026	0.6081	[0.3506 ; 0.8581]
		sAfr,20	0.2352	0.2286	0.2385	[0.1222 ; 0.2714]
		uAfr	0.0662	0.1105	0.0253	[0.0025 ; 0.4393]
		sEur,1	0.2917	0.3080	0.2203	[0.1048 ; 0.5951]
		sEur,20	0.0180	0.0189	0.0157	[0.0022 ; 0.0389]
		uEur	0.4250	0.3966	0.4567	[0.1077 ; 0.4950]

EurIN Scenario AfrIN -Eur2P Scenario AfrIN -EurIN Scenario AfrIN -EurDE

  S1.3) and an example of the patterns obtained for different u values is given in figure S1.3.

	1.0	Scenario Afr2P -Eur2P SOURCE EUROPE SOURCE AFRICA ADMIXED POP. H Scenario Afr2P -g = 0 Scenario Afr2P -EurDE s Afr,0 s Eur,0	u=0.00 u=0.05 u=0.10
	0.8				u=0.15 u=0.20
					u=0.25
					u=0.30
	0.6				u=0.35 u=0.40
	S,g				u=0.45
	s	Scenario AfrDE -Eur2P	Scenario AfrDE -EurDE	Scenario AfrDE -EurIN	u=0.50
	0.4
	0.2
	0.0
	0	5	10	15	20
			g
	Figure S1.3: Influence of u on equation S1.3 with s S,20 = 0.1 and s S,1 = 0.9
	Literature Cited in Supplementary Materials
	Rosenberg, N. A., 2004 DISTRUCT: a program for the graphical display of population structure.
	Molecular Ecology Notes 4: 137-138.
			2 3

h 1 h 20 s Eur,t Eur,p1 s Eur,t Eur,p2 t Eur,p2 t Eur,p1 g = 20 s Afr,t Afr,p2 t Afr,p2 t Afr,p1 s Afr,t Afr,p1

	SOURCE	ADMIXED
	AFRICA	POP. H
	g = 0
	s Afr,1

s Afr,20 g = 1 g = 20 s Afr,0 h 1 h 20 u Afr SOURCE EUROPE s Eur,0 s Eur,t Eur,p1 s Eur,t Eur,p2 t Eur,p2 t Eur,p1 g = 0

	SOURCE	ADMIXED
	AFRICA	POP. H

s Afr,1 s Afr,20 g = 1 g = 20 s Afr,0 h 1 h 20 u Afr SOURCE EUROPE s Eur,0 s Eur,t Eur,p1 s Eur,t Eur,p2 t Eur,p2 t Eur,p1

	SOURCE	ADMIXED	SOURCE
	AFRICA	POP. H	EUROPE
	g = 0	s Afr,0	s Eur,0

h 1 h 20 g = 20 s Afr,t Afr,p2 t Afr,p2 t Afr,p1 s Afr,t Afr,p1 s Eur,1 s Eur,20 u Eur

	SOURCE	ADMIXED	SOURCE
	AFRICA	POP. H	EUROPE
	g = 0	s Afr,0	s Eur,0

h 1 h 20 g = 20 s Afr,t Afr,p2 t Afr,p2 t Afr,p1 s Afr,t Afr,p1 s Eur,1 s Eur,20 u Eur

	ADMIXED	SOURCE
	POP. H	EUROPE
		s Eur,0

h 1 h 20 s Eur,1 s Eur,20 u Eur

	SOURCE		SOURCE	ADMIXED	SOURCE
		AFRICA		AFRICA	POP. H	EUROPE
	g = 0	s Afr,1 s Afr,0	g = 0	s Afr,1	h 1	s Eur,0
	g = 1
		u Afr
					h 20
	g = 20	s Afr,20

s Afr,20 g = 1 g = 20 s Afr,0 u Afr s Eur,1 s Eur,20 u Eur

	ADMIXED	SOURCE
	POP. H	EUROPE

s Eur,0 h 1 h 20 s Eur,1 s Eur,20 u Eur g = 0 SOURCE AFRICA g = 1 g = 20 s Afr,0 s Afr,1 s Afr,20 u Afr

	ADMIXED	SOURCE
	POP. H	EUROPE

s Eur,0 h 1 h 20 g = 0 SOURCE AFRICA g = 1 g = 20 s Afr,0 s Afr,1 s Afr,20 u Afr s Eur,1 s Eur,20 u Eur Figure 1 True model RF-ABC Predicted model

	AfrIN -EurIN	1.2%	2.7%	2.9%	2.9%	0.1%	10.8%	2.8%	9.7%	61.4%
	AfrDE -EurIN	1.6%	6.6%	0.8%	1.7%	1.7%	1.4%	9.5%	73.2%	17.7%
	Afr2P -EurIN	5.9%	2.2%	0.0%	5.1%	0.3%	0.0%	76.9%	8.8%	0.3%
	AfrIN -EurDE	2.0%	1.5%	9.9%	6.0%	2.0%	72.2%	0.9%	1.5%	18.4%	FST (Eur -Pop H)
	AfrDE -EurDE	5.9%	15.8%	2.0%	15.6%	77.7%	4.8%	1.4%	4.8%	1.7%	FST (Afr -Pop H)
	Afr2P -EurDE	11.2%	1.8%	0.4%	58.2%	7.6%	0.6%	6.9%	1.1%	0.1%
											f-3 (Pop H; Afr, Eur)
	AfrIN -Eur2P	5.7%	4.6%	76.3%	1.8%	0.4%	9.1%	0.0%	0.0%	0.2%
											Mean ASD (Pop H)
											Variance Heterozygosity (Pop H)
	AfrDE -Eur2P	11.2%	57.3%	6.7%	1.8%	7.6%	0.9%	0.6%	0.6%	0.2%	Mean Heterozygosity (Pop H)
											Mean ASD (Eur -Pop H)
	Afr2P -Eur2P	55.2%	7.6%	0.9%	7.0%	2.7%	0.2%	0.9%	0.2%	0.1%	Mean ASD (Afr -Pop H)
		A fr 2 P -E u r2 P	A fr D E -E u r2 P	A fr IN -E u r2 P	A fr 2 P -E u rD E	A fr D E -E u rD E	A fr IN -E u rD E	A fr 2 P -E u rI N	A fr D E -E u rI N	A fr IN -E u rI N

Mean Admixture from Afr in Pop H Variance Admixture from Afr in Pop H Skewness Admixture from Afr in Pop H Kurtosis Admixture from Afr in Pop H Mode Admixture from Afr in Pop H Minimum Admixture from Afr in Pop H 10%-Quantile Admixture from Afr in Pop H 20%-Quantile Admixture from Afr in Pop H 30%-Quantile Admixture from Afr in Pop H 40%-Quantile Admixture from Afr in Pop H 50%-Quantile Admixture from Afr in Pop H 60%-Quantile Admixture from Afr in Pop H 70%-Quantile Admixture from Afr in Pop H 80%-Quantile Admixture from Afr in Pop H 90%-Quantile Admixture from Afr in Pop H Maximum Admixture from Afr in Pop H

Mean Admixture from Afr in Pop H Variance Admixture from Afr in Pop H Skewness Admixture from Afr in Pop H Kurtosis Admixture from Afr in Pop H Mode Admixture from Afr in Pop H Minimum Admixture from Afr in Pop H 10%-Quantile Admixture from Afr in Pop H 20%-Quantile Admixture from Afr in Pop H 30%-Quantile Admixture from Afr in Pop H 40%-Quantile Admixture from Afr in Pop H 50%-Quantile Admixture from Afr in Pop H 60%-Quantile Admixture from Afr in Pop H 70%-Quantile Admixture from Afr in Pop H 80%-Quantile Admixture from Afr in Pop H 90%-Quantile Admixture from Afr in Pop H Maximum Admixture from Afr in Pop H Nine competing models Observed ACB Observed ASW AfrIN -EurIN AfrIN -EurDE AfrDE -EurIN AfrIN -Eur2P Afr2P -EurIN AfrDE -EurDE Afr2P -EurDE AfrDE -Eur2P Afr2P -Eur2P

	Supplementary Figure S3
	FST (Afr -Pop H)	FST (Eur -Pop H)	Mean ASD (Afr -Pop H)	Mean ASD (Eur -Pop H)	f-3 (Afr -Pop H -Eur)
	Mean ASD (Pop H)	Mean Heterozygosity (Pop H)	Variance Heterozygosity (Pop H)

kernel, and are constrained to the prior limits. Posterior parameter densities are indicated by a solid line; prior parameter densities are indicated by black dotted lines. and ADMIXTURE for the Barbadian (ACB) and the African American (ASW). 100,000 independent SNPs were considered from the 1000 Genome Project Phase 3 for 279 unrelated individuals (90 Yoruba (YRI), 89 British (GBR), 50 Barbadian (ACB), 50 African American (ASW)). (A) Allele Sharing Dissimilarity was computed between all pairs of individuals and the resulting matrix projected on the first two dimensions of a metric MDS. The twodimensional centroid of the Yoruba (YRI) and, respectively, the British (GBR) are indicated in red and connected by a black dotted line. ACB and ASW individuals are projected orthogonally onto this line and their relative distance to the Yoruba centroid is calculated to obtain ASD-MDS based individual admixture estimates. (B) A single run of unsupervised ADMIXTURE [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]) has been computed using the 279 individuals and 100,000 SNPs and results were plotted using DISTRUCT (Rosenberg 2004). Individual membership proportions to the "orange" cluster mostly represented by Yoruba (YRI) genotypes was considered as an estimate of African admixture for the ACB and ASW respectively. (C) Spearman correlation between ASD-MDS and ADMIXTURE-based estimates of African admixture for the ACB and ASW individuals separately. S3: Summary statistics prior-distribution densities for each nine competing models considered (Figure 1). 10,000 simulations were performed for each ninecompeting scenario and prior densities plotted with a different color indicated for each scenario.

Supplementary Figure

Corresponding statistics observed from the ACB and ASW population separately are represented, on each plot, by vertical doted-lines (red and blue respectively for ACB and ASW).

The 24 separate summary statistics considered are described in Material and Methods.

Supplementary Figure S4: Four first axes of the principal component analysis for the 90,000 sets of 24 summary statistics computed on simulated data under each nine-competing scenario (Figure 1). The 24 same statistics calculated for the observed ACB and ASW population samples, respectively, are then projected on the PCA and represented by, respectively, a red and blue star. All two-dimensional projections are orthonormal.

Supplementary Figure S5:

Histogram of the goodness-of-fit for the observed set of 24 summary statistics computed for (A) the ACB population, and (B) the ASW population, in turn serving as the observed admixed population H considering the YRI population sample as the African source and the GBR population sample as the European source (see Material and Methods). Goodness-of-fit statistics were calculated as the mean distance between observed and accepted summary statistics. Observed statistics are fitted to the full 90,000 sets of the same statistics calculated from 10,000 simulations performed under each nine-competing models (Figure 1). Goodness-of-fit was obtained considering 1,000 repetitions and a tolerance value of 0.01. S6: RF-ABC out-of-bag prior error rate as a function of the number of trees considered to build the forest for the model-choice procedure considering nine-competing scenarios (Figure 1).

Supplementary Figure

(parameters taken separately). For each target population separately and for each method, we conducted an out-of-bag cross validation by considering in turn 1,000 separate parameter inferences each using one of the 1,000 closest simulation to the observed ACB (or ASW) data as the target pseudo-observed dataset. All posterior parameter estimations were conducted using the other 99,999 simulations under the AfrDE-EurDE scenario (Figure 1, Table 1), a 1% tolerance rate (i.e. 1,000 simulations), 24 summary statistics, logit transformation of all parameters, four neurons in the hidden layer per neural network and 500 trees per random forest.

Median was considered as the point posterior parameter estimation for all parameters. First column provides the average absolute error; second column shows the mean-squared error;

third column shows the mean-squared error scaled by the parameter's observed variance (see

Material and Methods for error formulas).

Supplementary Table S1. Random-Forest Approximate Bayesian Computation model-choice predictions for the ACB and ASW populations. 1,000 decision trees were considered for RF prediction for the ACB and ASW respectively. Corresponding results are plotted in Figure 3. S2. Parameter prediction cross-validation error as a function of the number of neurons in the hidden layer and the rejection tolerance rate under the AfrDE-EurDE scenario. We considered, 1,000 random simulations in turn as pseudo-observed data to estimate posterior parameter distributions, considering 4, 5, 6, or 7 neurons in the hidden layer ("NN-HL" row), and 100,000 total simulations. Tolerance levels of 0.01, 0.05, 0.1 and 0.2 were considered ("Tolerance" row). The median values of posterior parameter distributions were used as point estimates for the error calculation. S3. Accuracy of the 95% credibility interval estimated for posterior parameters in the vicinity of the observed ACB and ASW datasets. We provide the empirical coverage of the estimated 95% credibility interval, i.e. how many times (in percentage) the true parameter (𝜃 𝑖 ) is found inside the estimated 95% credibility interval [2.5%quantile(𝜃 𝑖 ^) ; 97.5%quantile(𝜃 𝑖 ^)], among the 1,000 posterior parameter estimations conducted using in turn

Supplementary Table

Supplementary Table