Spatially explicit paleogenomic simulations support cohabitation with limited admixture between Bronze Age Central European populations

The Bronze Age is a complex period of social, cultural and economic changes. Recent paleogenomic studies have documented a large and rapid genetic change in early Bronze Age populations from Central Europe. However, the detailed demographic and genetic processes involved in this change are still debated. Here we have used spatially explicit simulations of genomic components to better characterize the demographic and migratory conditions that may have led to this change. We investigated various scenarios representing the expansion of pastoralists from the Pontic steppe, potentially linked to the Yamnaya cultural complex, and their interactions with local populations in Central Europe, considering various eco-evolutionary factors, such as population admixture, competition and long-distance dispersal. Our results do not support direct competition but rather the cohabitation of pastoralists and farmers in Central Europe, with limited gene flow between populations. They also suggest occasional long-distance migrations accompanying the expansion of pastoralists and a demographic decline in both populations following their initial contact. These results link recent archaeological and paleogenomic observations and move further the debate of genomic changes during the early Bronze Age.


Supplementary material 2: Model choice cross validation
This section presents the confusion matrices for the various model choice analyses. Lines of each table represent the scenario that have generated the pseudo-observed data while columns represent the relative probability estimated for each tested scenario. 100 replicates have been used in all cases. Confusion matrix of comparison between 8 scenarios ("Control", "F", "P", "F&P" with and without direct competition), without sampling layer optimization.  Table S2.4. Two scenarios with sampling layer optimization. Confusion matrix of comparison between the 2 scenarios without competition ("Control", "F&P") with sampling layer optimization.

Supplementary material 3: Sampling layer optimization
After the first ABC analysis, which consisted in comparing the relative probability of 8 scenarios where each sample could be sampled either in the pastoralists or in the farmers layer, we used the ABC analysis to fix some population layer when it was possible (see Material and Method) to optimize further estimation by maximizing the exploration of the parameter space. Here, we show that fixing the layer of four populations (those for which the results were unambiguous) improves the estimation.  The envelope represents 90% of the simulated data for both scenarios "control" and "F&P" with and without sampling layer optimization. The cross represents the observed data and is always inside the envelope. If the cross was outside the envelope, it would indicate a poor fit of the model. Estimation of the sampling layer based on scenario "F&P" without competition. The R "abc" package was used with loclinear method and tolerance 0.05. The mention "F or P" is returned when HDI99 lower value is lower than 0.5 and HDI99 upper value is greater than 0.5, meaning that the ABC estimation procedure does not allow to fix one layer. When both lower and upper HDI99 are lower than 0.5, value "F" (for farmers) is returned, while when both are greater than 0.5, value "P" (for pastoralists) is returned. The population Esperstedt_MN was not part of this layer optimization procedure as sampling in the "P" layer was never possible.

Supplementary material 4: Posterior distribution and cross validation
The next figures ( Figures S4.1-12) show the results of the parameter estimation for scenario "F&P" without competition, which is the most likely given the model choice procedure and with four fixed sample layers to optimize the results (see Supplementary material 3). On the left panel of each figure, the parameter estimation is represented as probability densities (prior before demographic filter: dotted line; prior after demographic filter: dashed line; posterior distribution: light grey). On the right panel of each figure, the cross-validation result represents the ability of the ABC analysis to retrieve a parameter value, by considering it as pseudo-observation and performing the ABC analysis. 100 replicates are done using the "abc" R package. Those results have been obtained under the best model (i.e. F&P, without competition and with the layer optimization).

Supplementary material 5: Long Distance Dispersal
We simulated the expansion of pastoralists to Central Europe under the scenario "control" without competition, with different rates of LDD: 0.0, 0.005, 0.01, 0.014, 0.2, 0.3, 0.4, 0.5. In that situation, we found that the absence of LDD make it impossible to sample the 6 early populations in the pastoralists layer, while the 7 th has less than 25% chance to not be sampled ( Figure S5: A). When simulating the scenario "F&P" without competition, we observed the same result in the pastoralist layer, but in addition, we also found that some populations taken in the farmer layer has chances not to be sampled ( Figure S5: B). This is because sampling is occurring at the time of the demographic drop in some populations and the number of individuals in that deme is not sufficient for sampling to occur. Those results show that a minimum number of LDD is necessary to reach all the sampling demes in Central Europe at the time where they have been sampled, under the most probable conditions (scenario "F&P" without competition).

Supplementary material 6: SPLATCHE3 proportions and f4-ratio
Even if they are expected to be closely linked, the genetic contribution from the pastoralist layer in SPLATCHE3 and the estimated genomic proportions from steppe ancestry associated with the population with Yamnaya culture from Haak et al. 1 are not computed in the same way. The former are exact virtual proportions, while the latter are estimated proportions from observed molecular data.
The purpose of this supplementary information is to show the relation between ancestry proportions calculated by the program SPLATCHE3 and the admixture proportions estimated based on f4-ratio statistics in the dataset of Haak et al. that we used as observed statistics in our ABC analysis. For this, we thus simulated two different models derived from our study (see below) and for each model we computed the correlation between SPLATCHE3 proportions and f4-ratios computed with the qpAdm function from Admixtools using the R package admixr 2 . Haak et al. used an unbiased estimation of those f4-ratios that we are not able to exactly reproduce with SPLATCHE3 because it requires to simulate several worldwide outgroup populations. Nevertheless, we assume that if f4ratios are correlated to SPLATCHE3 proportions, unbiased f4-ratios will be too as they are supposed to be closer to SPLATCHE3 proportions by definition.

Models investigated
We tested two different models to check the relation of SPLATCHE3 proportions and f4-ratios under different parameter context, one is based on the parameter prior distributions (Modelprior) and one is based on the posterior distribution of the model estimated as the best one in our study (Modelposterior).
Modelprior is based on the prior demographic distributions from our simulation framework (Table  3). It uses the middle of the prior range of each demographic parameter as input (including LDD, rF, rP, SDD, EDD, mF, mP). 10 individuals per population are then sampled under the "control" scenario without competition. The admixture rate ranges from 0.006 to 0.03 to explore whether the correlation between SPLATCHE3 proportions and f4-ratio is unbiased across its whole range.
Modelposterior is based on the posterior distribution estimated from our simulations of the scenario "F&P" without competition, which is estimated to be the most likely in our study. Each simulation is made of a random combination of demographic parameters taken in their posterior distributions estimated in our study, ranging from HDI90 lower to upper values ( Table 3). The number of individuals per population is identical to the dataset analyzed in our study and taken from Haak et al. 1 . Here the idea is to verify that the correlation between SPLATCHE3 proportions and f4-ratio is good in the range of estimated parameter values.
For both models, we tested for the correlation between steppe-ancestry proportions calculated with SPLATCHE3, and admixture proportions estimated with f4-ratios using 1,000 simulated SNP for 100 simulations per model case. The sampling dates are fixed to their average value to avoid more variability than the one already produced by LDD. Figure S6.1 shows a scheme of the hypothesized phylogenetic tree of the virtual populations under study. The f4-ratio is based on 5 populations, hereafter named X, A, B, C and O.

F4-ratio parameters
"X" represents the test populations under study. In the simulations, those populations are sampled in both population layers to show the differences.
"A" represents the farming population from Central Europe having no pastoralists ancestry, from which the f4-ratio statistics, associated to the test populations (X), is compared to. This population is sampled just before the arrival of pastoralists in Central Europe. The sampling date is the same than the real sample "Esperstedt_MN" (generation 193).
"B" represents the ancestral population from Central Europe associated to farmers. This population is sampled just after the moment first pastoralists appear in the simulation: 5,550 BP (generation 178).
"C" represents the ancestral population associated to pastoralists. This population is sampled in Southern Russia just after the moment first pastoralists appear in the simulation: 5,550 BP (generation 178).
"O" is the outgroup population, represented by the first farmers from the Middle East just after the beginning of the simulation: 9,950 BP (generation 2).
The resulting order of the populations that must be taken into account in the f4-statistics is thus: X, A, B, C, O. We can extract the pastoralists ancestry proportion from the resulting f4-ratio using the following equation: pastoralists ancestry proportion = 1 -f4-ratio

Pearson correlation
The Pearson correlation is a parametric test that gives a measure of the association between two continuous variables (here the f4-ratio and the SPLATCHE ancestry). A minimum of 8 values per variable are necessary to get a reliable statistical power. The Pearson correlation is calculated considering all population values (both tested measures are averaged per population).

Results
The correlation between SPLATCHE3 and f4-ratio values was tested for Modelprior with different γ values taken within its prior range. In addition, it was also tested for Modelposterior. The mean Pearson correlation is calculated with at least 17 Pearson correlation between f4-ratio and SPLATCHE ancestry (each one made of 1,000 SNP) per γ value (Table S6.1) and with 28 simulations with demographic parameters drawn randomly from the posterior distribution of the best scenario "F&P" (Table S6.1). The number of simulations is not always the same due to demographic filtering (i.e. simulations for which some population samples were not possible). The results show that the correlation decreases linearly with increasing γ (Figure S6.3) and is higher with a larger sample size (Modelprior against Modelposterior).
Moreover, the Figure S6    . "F" denotes a population sampled in the farmers layer, while "P" denotes a population sampled in the pastoralists layer. Ancestry refers to the steppe pastoralist ancestry.

Discussion
Unbiased admixture proportions estimated by Haak et al. 1 and used as observed data in our study are expected to be equivalent to SPLATCHE3 proportions. Our results show that there is indeed a significant correlation between both values (> 75% in all explored cases) and that they can thus be compared.
We observed that the Pearson correlation is higher with lower admixture rate γ than with larger one ( Figure S6.3). The γ posterior distribution that we estimated using the ABC analysis is estimated to 1%, corresponding to an expected Pearson correlation between SPLATCHE3 and f4-ratios of about 89% (Table S6.1 and Figure S6.3).
The correlation is better with Modelprior compared to Modelposterior because of the smaller variance in SPLATCHE3 proportions in the former due to larger sample size (10 genomes per sample instead of 1 to 10, see Table 1). Indeed, the variance of SPLATCHE3 proportions is much larger than the variance of f4-ratios due to the stochastic nature of the simulations, which is the reason why 100,000 simulations are used per scenario for the ABC analyses in our study.
As shown by Figure S6.4 and Figure S6.5, the f4-ratio tends to give admixture proportions slightly larger than the steppe-ancestry proportion obtained with SPLATCHE3. This bias increases with the values of γ (higher with larger γ) and decreases with the sample size.
Haak et al. improved their estimate of admixture proportion using 15 outgroup populations from different parts of the world, leading to "unbiased" estimates. Due to the complexity of simulating those outgroup populations, we did not reproduce their method exactly. Instead, we compared the steppe-ancestry proportions calculated with SPLATCHE3 to simple estimates of admixture proportions calculated with f4-ratio 3 , on which the estimate of Haak et al. is based, and we assume that simple f4-ratio and "unbiased" estimates based on those ratios are linked. In theory, "unbiased" population admixture based on f4-ratio should be closer to SPLATCHE3 steppe-ancestry proportion than simple f4-ratio.
Our comparison shows that if there is a bias, it should be in the direction of a slight underestimation of SPLATCHE3 steppe-ancestry proportions compare to classical f4-ratio. In any case, this bias should be comparable for all simulated scenario as they explore identical ranges of parameter values, so giving confidence to the model comparison. The parameter estimation could be slightly biased in favour of parameter values that increases the amount of pastoralist steppe-ancestry proportions in the sample data to compensate for the slight values increasing the proportion of pastoralist in SPLATCHE3. It means potentially a slight overestimation of the admixture rate and of the difference of carrying capacity in favor of the farmers over the hunter-gatherers and a slight underestimation of the amount of long-distance dispersal, which tend to limit genetic introgression from the local population to the invasive population.

Supplementary material 7: Representativeness of the data analysed
To check the robustness of the paleogenomic pattern investigated and the representativeness of the data analysed in relation to a larger dataset available from various publications, we retrieved from the AADR database (Allen Ancient DNA Resource, version 20 January 2021, https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-presentday-and-ancient-dna-data), all indexes (palaeogenomes) corresponding roughly to the period and area of interest, meaning from Central Europe between latitudes 46° to 53° and longitudes 8° to 22°, dated to the period from 5,300 BP to 2,800 BP (using the "mean BP" value).
By filtering with a mean coverage ≥ 0.05, we extracted 307 indexes in 105 GroupID (population samples). We filtered the dataset for duplicated, closely related and contaminated genomes by removing: i) indexes with same "Master_ID", keeping the one with the highest number of SNPs, ii) the "Master_ID" whose "GroupID" contains a "brother, sister, father, mother, son, daughter, sibling" content in the "GroupID" name, and iii) indexes for which the "GroupID" contains a "contam" or "dup" (contaminated or duplicated) in the "GroupID" name. There are 267 paleogenomes from 67 populations that remain after this filtering, with 1150639 SNPs (1102826 polymorphic). We then added a date to each population, considered later as "target population", calculated as the average of the "average BP date" of the genomes comprising the GroupID.
We then performed a qpAdm analysis using the program "Admixtools 2" (Maier & Patterson, https://uqrmaie1.github.io/admixtools/index.html), similarly to what was done to estimate the steppe ancestry proportions we are using in our study. To estimate the "steppe ancestry proportions", the qpAdm analysis was done independently on each retained population samples ("GroupID") from AADR. We choose the "outgroup populations" based on those in Haak et al. 1  We finally compared the datasets analysed in our study (34 paleogenomes in 10 populations, thereafter called "Haak et al 2015") with the filtered "AADR" dataset of 267 genomes from 67 populations. We divided the period of interest in 5 equal intervals of 500 years (the maximum time step to have at least one population sample from each dataset in each interval, see Figure S7.1). We then calculated the correlation between the average steppe ancestry computed in the two datasets at each time interval. We computed a one-tail Pearson correlation test by using the R package pwr version 1.3, R version 4.0.5. We found a strong positive correlation between the two datasets, statistically significant at the 5% threshold (r = 0.826, P value = 0.043). Note that 24 paleogenomes are common to both datasets.
Our analysis shows that the dataset analysed in our study (34 paleogenomes) is robust to the current state of paleogenomic data revealing a general pattern of a sudden increase in steppe ancestry in Central Europe during the Bronze Age.

Supplementary material 8: Conversion of LDD parameters
In the current study, we based the prior values of the long-distance dispersal (LDD) parameters on those estimated by Alves et al. 4 , but we needed to convert some of them due to the difference in deme size. Indeed, here we used demes of 100*100km while Alves et al. used demes of 150*150km. Below are the different LDD parameters in SPLATCHE3: • α: "Gamma Distribution shape parameter mode".
• µKm: average distance of LDD in km.
• σKm: standard deviation of the LDD distance in km.
α and µ are retrieved directly from Alves et al. (α150=1.209,µ150=5.357) but the other parameters are calculated based on these those two but adapted to demes of 100*100km.
We use the following formulas to retrieve β, σ, µKm, parameters from α and µ.

• µKm=µdist*(dist); σKm=σdist*(dist); ModeKm=Modedist*(dist)
We convert the average number of demes of 150*150 km into distances in km: After conversion, the Gamma Distribution corresponding to demes of 150 km by 150 km or 100 km by 100 km have identical distribution densities, only the scale changes ( Figure S8.1). The number of demes that can be migrated considering the gamma distribution for demes of 100 km by 100 km has grown by 150%. Thus, new formulas can be retrieved from those results, to facilitate the conversion: • distA: deme length before correction, • distB: deme length after correction • param: parameter to be corrected, being Mode, µ, σ or 1/β=θ. The parameter α, that controls the shape of the distribution, is not changed after correction ( Figure S8.1).
• paramdistB = paramdistA * distA/distB Figure S8.1: Distribution of the long-distance migration probability (represented with a Gamma density) before and after conversion. The shape is the same, only the range of the number of demes is 1.5 times higher after conversion (due to 1.5 times smaller demes side length).

Supplementary material 10: Random choice of 100'000 simulations
This section provides the same analysis as in the main text with a different random subset of 100,000 simulations to show the robustness of model choice and parameter estimation to this subset procedure. The absolute difference between the two subsets is always smaller than 0.04 and the probability order of scenarios is the same (Table 10.1). The estimators for the well estimated parameters (γ, LDD and sDD) is very close with less than 11% relative difference (Table 10.2, while this relative difference increases for the parameters with less precise estimators such as NmF and NmP (14.5% -18.4%).
Table S10.1 Posterior probability for the comparison of 8 scenarios ("Control", "F", "P", and "F&P" with and without direct competition) without sampling layer optimization. To be compared to the Table 2 in the main text. The difference is computed with the posterior probabilities given in the main text using a different random subset.

Supplementary material 11: Choice of the number of simulated SNP
To be sure that the number of simulated SNP is high enough to get robust estimates of genomic proportions from source populations (pastoralist of farmer) with SPLATCHE3, we repeated the most probable Scenario "F&P" by simulating 1,000 SNP. Then we compared the resulting genetic proportions of pastoralist with the one obtained with a subset of only 50 SNP ( Figure S11.1). It shows that, overall, the steppe-ancestry proportions calculated with 50 SNP give the same results than if it was calculated with 1,000 SNP. We found that the only population showing slightly different steppe-ancestry proportions is the one dated to 226 generations, corresponding to the "Alberstedt" sample. In that case, the sampling occurs right after the demographic drop (starting at 222 generations) so the combination between smaller population sizes and randomness of LDD events increase the variance of ancestry proportions. However, we judged that the gain in computational time due to simulation of only 50 SNP instead of 1,000 SNP (>20x), allowing to better explore the parameter space, worth this slight approximation in rare cases.

Supplementary material 12: Choice of tolerance rate for ABC
We also investigated variation of the choice of the tolerance rate for the ABC analyses. We tested different values of tolerance rates between 0.001 (100 retained simulations) and 0.2 (20,000 retained simulations) when comparing the four scenarios without competition (Figure S12.1 and S12.2). We found an instability of the results with low tolerance rate (< 0.01), due to the too low number of retained simulations ( Figure S12.1) and a diminution of the goodness of fit measure by the GOF p-values with larger number of tolerance rates (> 0.01). Consequently, we decided to use a tolerance rate = 0.01 (1,000 retained simulations) for all analyses presented in the main text.