Estimating recent migration and population-size surfaces

In many species a fundamental feature of genetic diversity is that genetic similarity decays with geographic distance; however, this relationship is often complex, and may vary across space and time. Methods to uncover and visualize such relationships have widespread use for analyses in molecular ecology, conservation genetics, evolutionary genetics, and human genetics. While several frameworks exist, a promising approach is to infer maps of how migration rates vary across geographic space. Such maps could, in principle, be estimated across time to reveal the full complexity of population histories. Here, we take a step in this direction: we present a method to infer maps of population sizes and migration rates associated with different time periods from a matrix of genetic similarity between every pair of individuals. Specifically, genetic similarity is measured by counting the number of long segments of haplotype sharing (also known as identity-by-descent tracts). By varying the length of these segments we obtain parameter estimates associated with different time periods. Using simulations, we show that the method can reveal time-varying migration rates and population sizes, including changes that are not detectable when using a similar method that ignores haplotypic structure. We apply the method to a dataset of contemporary European individuals (POPRES), and provide an integrated analysis of recent population structure and growth over the last ∼3,000 years in Europe.

Taking the limit ∆t → 0, we arrive at the density The random walk approximation to the coalescent 20 Here, we introduce an approximation, 21 P (X i (t) = κ, X j (t) = κ) ≈ P (X i (t) = κ)P (X j (t) = κ).
The intuition is that the probability that lineage i and j coalesce before time t is extremely 22 small such that the two lineages approximately behave like two independently moving par-23 ticles. Each lineage can be modeled by a random walk with transition matrix M . These 24 assumptions were also made in the context of continuous spatial diffusion models for hap-25 lotype sharing Baharian et al. (2016); Ringbauer et al. (2017) , and even further back, as a 26 general approximation to the two-dimensional continuous-space coalescent process (Barton 27 et al., 2002;Wilkins, 2004;Blum et al., 2004;Novembre and Slatkin, 2009;Robledo-Arnuncio 28 and Rousset, 2010). 29 This approximation implies that where lineages i, j are initially sampled in deme α, β. Or equivalently in matrix form, where Q = diag(q 1 , ..., q d ).

32
Corollary 1.1.1 Let time slice k be defined by the interval t k−1 < t < t k , M k denote the migration rate matrix in time slice k, and Q k = diag(q k 1 , ..., q k d ) where q k α denotes the coalescent rate in deme α at time slice k. Let T i,j denote the coalescent time between lineage i, j sampled in demes α, β, then under the independence assumption, for t ∈ (t K−1 , t K ), Expected number of lPSC segments given the demography Θ We obtain the desired result by substituting (S11) and (S12) into (S10) and canceling like-44 terms.

45
Expected age of a segment 46 We choose PSC segment lengths based on their expected age which is derived below.

47
length L 1 centiMorgans and L 2 centiMorgans is approximately 300 population size (N ) is sufficiently large. 50 We choose to work in units of basepairs, and will convert back to units of morgans at the 51 end. We convert L 1 into units of base-pairs with the transformation: µ = L 1 100r and similarly 52 ν = L 2 100r .

53
Let us denote T |l, N as the random coalescent time of a PSC segment that is at least length l under a single-deme demography model with population size N . The expected coalescent time of an PSC segment longer than µ base-pairs can be expressed as where f L (l|t) = 4r 2 t 2 le −2trl denotes the probability density that a PSC segment is of length 54 l given it has a common ancestor event at time t, f T (t|N ) denotes the probability density 55 that a coalescent event occurs at time t under the demography model with population size 56 N .

57
Next, we expand a key term in equation (S13) and assume, Putting everything together, (S16) We can remove the dependence of N by taking lim N →∞ as done similarly in Baharian et al.

59
(2016), (S17) Now that we have derived the expected age of PSC segment longer than µ, it is quite simple to expand the equation for PSC segments between µ and ν base-pairs, We transform back to units of centimorgans: let L 1 = 100rµ and L 2 = 100rν be in units of 61 centiMograns, we get the desired result where N (t) is the number of steps taken by time t, and Z i is a random variable representing 71 the direction and magnitude taken at step i. Since X(t) is a sum of iid variables, a form of 72 the central limit theorem applies here and X(t) converges to the normal distribution (Rényi, .

74
In a random walk on a triangular grid, a particle can move in one of the 6 directions and, where I 2 is the identity matrix. We do not show the additional steps of using the law of Here, we discuss how we interpret the diffusion constant m(∆x) 2 . Let the distance d = X(t) = x 2 1 + x 2 2 , then Thus the square root of the diffusion constant can be written as suggests an interpretation as the distance traveled by an individual after one generation, and 85 sometimes is called the "dispersal" distance or the "root mean square dispersal distance".

86
In these units, dispersal is more directly comparable to empirical estimates of dispersal from

Diversity rates versus coalescent rates
tance model and within-deme "diversity rates" to approximate expected pairwise coalescent 93 times, in which, whereÊ[T α,β ] is the resistance distance approximation to the expected coalescent time between deme α and deme β, e qα is the "diversity rate" in deme α, and R α,β is the resistance distance between demes α, β (Petkova et al., 2016). The diversity rates have no simple expression in terms of population-genetic parameters under the multi-deme coalescent model. As an alternative, diversity rates can be interpreted as reflecting average within deme heterozygosity since e q = E[T w ] ∝ H α where the heterozygosity for deme α (H α ) is defined as, where D i,j is the average number of differences between (haploid) individuals i and j.

96
MAPS models the recombination process using rates estimated from a recombination rate 97 map. In this model, population sizes and migration rates can be inferred separately rather 98 than as a joint parameter. Intuitively, the recombination rate serves an independent clock 99 to calibrate estimates.

100
More formally, a statement of identifiability is a statement regarding the likelihood.

101
MAPS models the expected number of lPSC segments shared between pairs of (haploid) 102 individuals, and can be computed with an integral. The integral can be broken up into a 103 product of two functions: a function describing the decay of PSC segments as a function 104 of time ("recombination rate clock"), and the coalescent time probability density f T i,j (t).

105
The migration rates and population sizes only appear in f T i,j (t), and cannot be factored into 106 parameters involving combinations of the migration rates and population sizes. tessellation for the coalescent rates is T q = (l q , q, c q , µ q ).

113
The location of each (unordered) Voronoi cell is distributed uniformly across the habitat, where U denotes the uniform distribution. The number of cells (a-priori) are drawn from a 115 negative binomial distribution, The effects of each Voronoi cell is normally distributed with variance ω 2 .
The probability of a particular (unordered) cell configuration is, We assume, We set 1.5 as the upper bound for log10(ω m ). For the normal distribution, 95% of the density 120 is within two standard deviations of the mean. As a result, by setting the upper bound for 121 for log10(ω m ) to be 1.5, we are constraining m so that the probability that it is within 3 122 orders of magnitude from the mean is 0.95 a priori, and similarly we set 1 as the upper 123 bound for log10(ω q ) to restrict the population sizes so to be within 2 orders of magnitude 124 from the mean with probability 0.95 a priori. 125 We assume, 126 µ m ∼ U (−10, 4) (S33) µ q ∼ U (−10, 4).

(S34)
We place a uniform prior on the log of the mean rates to reflect that we are uncertain about 127 the order of magnitude. Here, the data is highly informative of the mean, as a result, we 128 can allow the support of the prior to vary by many orders of magnitude.
We can write the cell specific migration rates as, Written this way, we can make separate updates to the cell effects (e i ), the mean (µ) and zero, we add MH joint random-walk updates to µ and e i as follows, where ∼ N (0, 1). The intuition here is that a constant is subtracted from the mean, 143 which then gets added to the individual cell effects e i which ensures that the mean cell effect 144 is approximately 0.

145
The steps above are applied to both the migration rates and coalescent rates.

146
Updating the number of cells 147 The number of cells change the dimension of the likelihood, and a result, we must use 148 a Reversible Jump MCMC step so that the ratio of densities in the Metropolis-Hastings 149 acceptance ratio is well-defined (Green, 1995). We choose to update the number of cells with 150 a birth-death update (Stephens, 2000). Fortunately, in such a case, the updates reduce to 151 standard Metropolis-Hastings because the dimension matching constant (i.e. the "Jacobian") 152 equals one (Petkova et al., 2016;Stephens, 2000). See equations S31 and S32 in Petkova where x denotes the current state of the MCMC, x the proposed state, e c+1 is the proposed cell effect drawn from a standard normal, l() is the likelihood function, and p() is the prior. Conversely, in a death-update, we randomly choose one cell uniformly to kill. In this case, the acceptance ratio for a death proposal (going from c + 1 cells to c cells) is α(x, x ) = min(1, p(birth) p(death) l(x )p(x )N (e c |0, 1) l(x)p(x) 1 c+1 ).