Estimating mobility using sparse data: Application to human genetic variation

Geographic, Temporal, and Trait Distances.

The geographic distance between all sample pairs was calculated in kilometers using the Haversine Eq. (51) to account for the curvature of the Earth as follows:Gij=2r arcsin(sin((φj−φi)/2)2+cos(φi)cos(φj)sin((λi−λj)/2)2),[1]where G is the distance in kilometers between individuals i and j; φ_i and φ_j are the latitude coordinates of individuals i and j, respectively; λ_i and λ_j are the longitude coordinates of individuals i and j, respectively; and r is the radius of the earth in kilometers.

Temporal distances between samples were calculated as time in years between sample pairs. Previously reported date ranges based on stratigraphy or direct radiocarbon dating were used for all individuals (SI Appendix, Table S1). In all analyses, sample dates were randomly drawn from a uniform distribution corresponding to the upper and lower bounds of a time period for a given specimen.

Genetic distances were calculated as pairwise proportion of alleles that are not identical by state (pairwise heterozygosity), using the function ibs.dist from the Bioconductor package SNPstats v.1.18.0 (52) in the R statistical analysis environment v3.2.2 (53).

The S_max Estimator.

To consider the full range of scaling factors on a finite interval, we choose to represent S as the tangent of an angle α between 0° and 90°, where α = 0 corresponds to S = 0 (geographic variation alone explains the observed trait distances between entities) and α = π/2 corresponds to S = ∞ (temporal variation alone explains the observed trait distances between entities). Formally, the time–space product matrix (D) was calculated as follows:Dij=Gij2+(STij)2,[2]where i and j are the specimens considered, D is the time–space product matrix, G is the geographical distance matrix, T the temporal distance matrix (given by the difference in age of the samples), and S is the scaling factor [S = tan(α)].

To find the scaling factor, S_max, that maximizes the correlation between the trait distance matrix and D, the time–space product matrix, we calculated the Pearson correlation coefficient between these matrices for 200 (500 for the simulated data) scaling factor values (Fig. 1). The scaling factor value in the time–space product matrix that produced the strongest correlation with the trait distance matrix is recorded as S_max, the mobility estimator.

Simulation Tests.

The reliability and the robustness of the S_max statistic in recovering information about past mobility was explored using a spatiotemporally explicit simulation model. The simulation world consists of a grid of 8,000 by 8,000 demes. Each simulation starts with one entity placed in a randomly chosen deme and lasts 20,000 generations. The model simulated exponential population growth to a carrying capacity of 10,000 entities, followed by a stochastic birth–death process (54) mobility, and trait mutation. We generated spatiotemporal trait variation data under different mobility parameter values using the same S_max estimation protocols as described above for each dataset. A total of 10,000 independent replicates of the simulations and analyses were generated, and the utility of the S_max statistic in recovering information about mobility was assessed by correlation.

The migratory process was modeled as Gaussian random walks: In each generation, each entity moves independently in the x and y directions by distances picked randomly from a normal distribution with mean = 0 and SD = σ_mig. This corresponds to the average distance moved in a single step (d_mig) of π/2 σ_mig = 1.2533 σ_mig. Thus, d_mig is the parameter of interest. We choose 1,000 random values of d_mig from a uniform distribution with a range of 1–100. We modeled drift as a Moran-type birth–death process (54). At each generation, each entity undergoes binary fission with probability P = 0.1, creating a duplicate of itself at the same location. The two entities subsequently move and evolve independent of each other. When the number of entities reaches or exceeds the carrying capacity (10,000), excess entities are deleted at random among all entities present in that generation. Mutation was modeled as a one-dimensional Gaussian random walk for each trait (N_traits = 50). Each trait was assigned an initial value of 1,000, and new (mutated) values were picked from a random normal distribution with mean equal to the current value and SD fixed at 0.05.

Following a burn-in period of 10,000 generations, entities were sampled from simulations with a probability of 0.00001 at each generation. The x and y coordinates, time of sampling in generations, and the values for the 50 traits were recorded for all sampled entities.

Pairwise trait distances between all sampled entities in each of the simulated datasets were calculated using the Euclidean distance formula as follows:Mij=∑k=1n(dik−djk)2,[3]

where M_ij is the distance between the two entities i and j; d_ik and d_jk are the values of the trait k for individuals i and j respectively; and n is the number of recorded traits.

Out of 10,000 simulations, 9,866 (98.66%) resulted in an EC greater than zero. To match the simulated data with the empirical data, we filtered the simulated data based on the measured EC values and removed all simulations that produced an EC value smaller than 0.001. This resulted in 9,155 simulations being used in the correlation analysis.

To assess the reliability of the S_max statistic in recovering information about mobility, R² values were calculated for the correlation between the simulated d_mig values and their corresponding S_max values.

Human Mobility in Late Pleistocene and Holocene.

We considered genome-wide data comprising 354,199 SNPs typed in 329 West Eurasian (i.e., west of the Ural mountains) individuals (SI Appendix, Fig. S2) temporally ranging from ∼39,000–1,000 y before present (see SI Appendix, Fig. S1). We merged the overlapping SNPs typed in archaeological samples published in refs. 22⇓⇓⇓–26 and 55⇓⇓⇓⇓⇓–61 (see SI Appendix, Table S1 for list of samples and references) that met the geographic and temporal criteria described above. No additional bioinformatic processing of the data were carried out for this study.

The 329 individuals were assigned to one of following three groups based on their estimated age and subsistence strategy based on their archaeological context: pre-LGM hunter–gatherers n = 19 (temporally ranging from 39,000 y B.P. to 26,000 y B.P.), post-LGM hunter–gatherers n = 47 (temporally ranging from 19,000 y B.P. to 5,000 y B.P.), and Holocene farmers n = 263 (temporally ranging from 10,000 y B.P. to 500 y B.P.).

Sliding window analysis was performed on all individuals in the dataset postdating 16,000 y B.P. The S_max statistic was estimated for 121 overlapping 4,000-y windows, each differing by 100 y.

To take age uncertainty into account, we report the mean scaling factor angle from 10,000 replicates with sample dates randomly resampled from their age ranges. We estimated 95% confidence intervals through a jackknifing procedure in which a randomly chosen sample in each window was removed from analysis, and the 0.025 and 0.095 quantiles were calculated from the resulting distribution.

To estimate the IBD signal through time, we fitted a linear model of genetic distances as a function of geographic distances in each time window (with sample jackknifing and age resampling as before, using the lm function from the R package base version 3.2.2) (53) and reported the slope of the line.

Confidence Intervals and Robustness of S_max Estimator.

We tested the assumption that there is an IBD pattern by correlating the genetic (trait) distance matrices in all time bins and in all windows with the respective geographic distance matrices and the date-resampled temporal distance matrices and calculated the P values for these correlations. We find a positive and statistically significant IBD pattern in space in all windows (SI Appendix, Fig. S3 A and B, respectively and SI Appendix, Fig. S6). The isolation by temporal distance pattern is positive and significant for most windows, but some windows show negative correlations or are not significant. We find that these windows correspond to time periods where we observe low EC (SI Appendix, Fig. S3C) and also low P values for the EC (Fig. 4B).

To account for the uncertainty in sample ages, we calculated the scaling factor angle 10,000 times using dates resampled at random from a uniform distribution for each sample, as described above, and report the average of the scaling factor angle of the given distribution as a point estimate.

We also performed a leave-one-out analysis (10,000 replicates, combined with sample date resampling) to explore the combined effect of sampling and dating uncertainty and constructed approximate equal-tailed 95% confidence intervals for all groups and windows.

To assess the statistical significance of S_max estimates, we consider the EC—defined as the Pearson correlation coefficient between the trait difference matrix and the time–space product matrix when S = S_max, minus the Pearson correlation coefficient between the trait difference matrix and either the temporal or geographical distance matrix alone, whichever is higher.

To obtain a null-distribution of EC, we permuted trait data for individuals among the spatiotemporal sample locations 10,000 times and calculated EC for each permutation, as described above. Finally, we calculate the proportion of EC values from the permuted datasets that are equally high or higher than that obtained from the observed data. This permutation test permits assessment of how frequently the EC for the observed data are produced by chance alone or, alternatively, as the result of the method used for estimating the S_max statistic. The resultant P value is the probability of observing an equally high or higher EC value in permuted, supposedly signal-less data, and provides an indication of the information content of each dataset.

Simulated Scenario of Changing Migration Rate.

We modified our simulations to represent a population experiencing two changes in migration rate, resulting in three episodes of constant migration rate. We assumed a generation time of 25 y and chose the effective population size to be 2N_e = 10,000, standard figures in population genetic models of European populations (62). We next chose three levels of migration with relative magnitude on par with what was inferred from the empirical data: m₁ = 0.0002, m₂ = 0.01, and m₃ = 0.05. To ensure equilibrium conditions during the start of the sampling period, we discarded the first 10,000 steps of the simulation (using migration rate m₁). We then simulated a time period of 20,000 y, divided into three episodes with constant migration rate: m₁ for 25,000–15,000 y ago, m₂ for 15,000–10,000 y ago, and m₃ for the last 5,000 y of the simulation. This roughly corresponds to the time spans associated with Mesolithic hunter–gatherers, Neolithic farmers, and post-Neolithic cultures in our empirical dataset. From a population genetic point of view, whole genome data as used in the empirical estimates correspond to a large number of approximately independent replicates. Because our model does not include recombination, we accounted for this effect by increasing the sample size to 10,000 individuals. SI Appendix, Fig. S4 shows the migration rate estimation using the S_max statistic using a 4,000 y-wide sliding window.

R version 3.2.2 (53) was used for analyses throughout this manuscript. The correlations between temporal, geographic, and trait distance matrices were calculated using the mantel (method = “pearson”) function in R package Vegan version 2.3.0 (63). The permutation and bootstrap tests were performed using the function sample in the R package base version 3.2.2 (53).