Vasili Pankratov
email: vasilipankratov@gmail.com

Francesco Montinaro

Alena Kushniarevich

Georgi Hudjashov

Flora Jay

Lauri Saag

Rodrigo Flores

Davide Marnetto

Marten Seppel

Mart Kals

Urmo Võsa

Cristian Taccioli

Märt Möls

Lili Milani

Anto Aasa

Daniel John Lawson

Tõnu Esko

Reedik Mägi

Luca Pagani

Andres Metspalu

Mait Metspalu

Differences in local population history at the finest level: the case of the Estonian population Running title: Genetic structure of Estonia

Several recent studies detected fine-scale genetic structure in human populations. Hence, groups conventionally treated as single populations harbour significant variation in terms of allele frequencies and patterns of haplotype sharing. It has been shown that these findings should be considered when performing studies of genetic associations and natural selection, especially when dealing with polygenic phenotypes. However, there is little understanding of the practical effects of such genetic structure on demography reconstructions and selection scans when focusing on recent population history. Here we tested the impact of population structure on such inferences using high-coverage (~30X) genome sequences of 2,305

Estonians. We show that different regions of Estonia differ in both effective population size dynamics and signatures of natural selection. By analyzing identity-by-descent segments we also reveal that some Estonian regions exhibit evidence of a bottleneck 10-15 generations ago reflecting sequential episodes of wars, plague, and famine, although this signal is virtually undetected when treating Estonia as a single population. Besides that, we provide a framework for relating effective population size estimated from genetic data to actual census size and validate it on the Estonian population. Our results suggest that the history of human populations within the last few millennia can be highly region-specific and cannot be properly studied without taking local genetic structure into account. Our approach to estimating the census population size may be widely used both to cross-check estimates based on historical sources as well as to get insight into times and/or regions with no other information available.

Main text

As more and more datasets including genetic data from hundreds and thousands of individuals become available it becomes apparent that most if not all human populations exhibit at least some degree of geography-driven genetic structure even at small scales (for some examples see [START_REF] Leslie | The fine-scale genetic structure of the British population[END_REF][START_REF] Martin | Haplotype Sharing Provides Insights into Fine-Scale Population History and Disease in Finland[END_REF][START_REF] Bycroft | Patterns of genetic differentiation and the footprints of historical migrations in the Iberian Peninsula[END_REF][START_REF] Raveane | Population structure of modern-day Italians reveals patterns of ancient and archaic ancestries in Southern Europe[END_REF][START_REF] Pierre | The genetic history of France[END_REF] ). Many recent publications have shown the confounding effect of such population structure on studies of genetic associations and natural selection signals, mainly in the case of polygenic phenotypes [START_REF] Berg | Reduced signal for polygenic adaptation of height in UK Biobank[END_REF][7][8][9] . Here we study the fine-scale genetic structure of the Estonian population and the local differences in recent demographic history and action of natural selection between genetically defined Estonian subgroups to gain a deeper understanding of the forces shaping this population structure and the consequences it has for population genetics analyses. In doing so we make use of high-coverage whole genome sequences from more than 2300 Estonian individuals generated within a different study 10 .

Our exploratory principal component analysis (PCA) (Figure 1) shows the presence of genetic structure within Estonia with the main differentiation between South-East and North-East of the country in agreement with previous studies [START_REF] Martin | Haplotype Sharing Provides Insights into Fine-Scale Population History and Disease in Finland[END_REF]11,12 . To zoom-in into the fine-scale structure in Estonia, we used total genetic length of shared IBD segments detected with IBDseq 13 as input for the fineSTRUCTURE 14 clustering algorithm (Methods). We applied this approach to a subset of 468 individuals sampled in rural areas at the age of 50 or more, as this cohort is expected to be the least affected by recent migrations (Figure 2, SI1:2.3). We refer to this subset as "R50+" throughout the text (Methods).

IBD-based analysis (Figure 2) reinforces previous observations [START_REF] Martin | Haplotype Sharing Provides Insights into Fine-Scale Population History and Disease in Finland[END_REF]11,12 , including the strong differentiation between South-East and the rest of Estonia, and provides a deeper insight into Estonian genetic structure, showing that most of the revealed clusters are highly geographically localized. The sharing matrix provides additional details. First, off-diagonal sharing also reflects geography with clusters from the same area tending to have higher inter-cluster sharing. Second, intra-cluster sharing substantially varies among clusters, implying differences in effective population size (Ne), which is also supported by the results of homozygosity-by-descent analysis (Figure 3).

In order to understand how gene flow barriers and/or differences in local population density shaped the IBD-sharing pattern in the R50+ dataset, we inferred migration surfaces using MAPS 15 . We used two windows of IBD segments length (in centimorgans, cM), 2-6 cM and more than 6 cM, which under a simplistic model of infinite population size have mean segment ages of 50 and 12.5 generations respectively 15 . Results for the two length bins generally agree with each other, suggesting higher levels of gene flow in the North along with a barrier separating South-East Estonia (SI1:2.4). A second barrier, separating the islands, especially Hiiumaa, from the mainland is also evident. This observation suggests that the population ancestral to modern South-East Estonians was partially isolated from the rest of the country at least since 50 generations ago. Interestingly, this genetic differentiation is consistent with linguistic data suggesting that the deepest split within the Finnic languages separates Southern Estonian from the other branches of the phylum that includes Northern Estonian 16 .

As local differences in admixture with external populations may have played a role in creating the observed genetic structure within Estonia we looked at patterns of haplotype sharing between R50+ Estonians and different non-Estonian populations (Table SI2:3

.1-I).
Here we used a conventional CHROMOPAINTER/fineSTRUCTURE/GLOBETROTTER (CP/FS/GT) approach 17 (Methods). Figure 4 shows the results of non-negative least squares (NNLS) [START_REF] Leslie | The fine-scale genetic structure of the British population[END_REF] , modelling each individual from the R50+ dataset as a result of admixture between non-Estonian groups revealed by CP/FS (Figure 4, SI1:3.1, Table SI2

-3.1-IV).
Admixture signals in Figure 4 show clear geographic patterns that match known historical evidence of external migration to Estonia, including Swedish settlements on the western coast and islands in 14-15th centuries and Finnish immigration to North-East Estonia in the 17th century 18 . Comparing NNLS results between clusters from Figure 2 we found that some of them, such as NE_1 and NE_2, stand out in terms of sharing with external groups but most of the clusters have overlapping distributions of NNLS scores (SI1:3.1). A similar pattern is observed in IBD-sharing patterns (SI1:3.2). These results suggest that admixture with non-Estonian groups can only partially explain the fine genetic structure observed in Figure 2.

We show that, despite the small territory it occupies, the Estonian population exhibits a readily detectable genetic structure, reflected in patterns of IBD segments sharing (Figure 2) and allele frequencies (Figure 1, Table SI2 -2.3-III, Table SI2-2.3-IV). Next, we sought to explore whether this differentiation has any effect on the reconstruction of demographic processes, namely whether there are region-specific differences in effective population size dynamics and action of natural selection. We hence applied IBDNe, which estimates effective population size (Ne) in past generations 19 , and SDS (Singleton Density Score), a tool for detecting signatures of natural selection 20 , as both methods give insight into very recent time periods, when regional differences in population history may be anticipated. For both analyses, we used the entire dataset of 2,305 samples, for which clusters were inferred using the same approach as for the R50+ subset (Figure 5).

We ran IBDNe 19 on the four most distinct clusters from Figure 5, representing four regions of Estonia: North-West, North-East, South-West and South-East and observed rather distinct Ne trajectories (Figure 6a, SI1:4.2). In particular, all clusters (except for eSE_5) show evidence of an effective population size decline between 10 and 20 generations ago, which is not detected when the entire dataset is analyzed (Figure 6a). Overall, these results suggest that population dynamics are region-specific and hence population-wide result may depend on the sampling scheme. For a deeper understanding of this phenomenon and the effects of other factors on IBDNe results, we applied IBDNe to genetic data, simulated under various demographic scenarios (SI1:4.1). Furthermore, the same approach has been applied to genotype data from the UK population, where regional differences in Ne dynamics are observed as well (SI1:4. Given our understanding of confounders of the observed regional Ne patterns, we exploited the fine-grained temporal resolution enabled by IBDNe to infer changes in actual census sizes (Nc) of the ancestors of contemporary Estonians, adapting previous theoretical work 21 to empirical case of human populations (Methods). We applied equation [3] (Methods) to the Estonian-wide Ne trajectory inferred using the Est1527 subset, which excludes clusters that can be considered as outliers in terms of external admixture and/or Ne trajectory (SI1:4.4).

We then compared the inferred Nc with available historical estimates (Figure 6b) showing a remarkable match between the two with the exception of the last three generations, for which IBDNe estimates are extrapolated from preceding time points 19 . This match may be attributed to i) our success in adequately controlling for events of recent gene flow and population structure; ii) the relatively recent time intervals considered, which limits the range of spatial interaction among the ancestors of contemporary Estonians. However, note that the pronounced fluctuations in Nc reported by historians between 1500 and 1700 are only very roughly approximated by the Ne-derived curve which, as expected 22 , provides only relatively long-term harmonic average of Ne. Nevertheless, we suggest that after controlling for confounders such as population structure and admixture and keeping in mind all the assumptions implied by the biological notion of Ne, our approach could be used to convert Ne to human Nc at any time interval for which historical records are missing, including the ones provided by PSMC 23 , which are beyond the scope of the current paper.

We then questioned whether natural selection could have also acted differently within the Estonian population. In doing so, we applied singleton density score (SDS) 20 to the entire dataset of 2,305 samples as well as to two regional subsets, South-East Estonia (SE, consisting of 1,029 samples belonging to clusters eSE_1 -eSE_5 in Figure 7) and the remaining 1,276 samples from the rest of the country (nonSE) (Methods, SI1:5.1).

First, we inspected the genome-wide distribution of positive SDS scores in the three datasets (Figure 7) for any evidence of recent selection acting at individual loci.

Unlike other studies that used SDS 20,24 we don't observe any hits with very low p-value (possible reasons are discussed in SI1:5.3). However, we see that the distribution of SDS scores differs between the three datasets (Figure 7, Table SI2 -5.3-I). Whereas one genomewide significant hit (rs75386033 and rs79907158 on chromosome 6) is detected in the SE, nonSE and the entire dataset had many more hits with p-values in the range between 5x10 -8 and 1x10 SI1:5.3,SI1:5.4,. Whereas most of the top SDS signals do not overlap between the three datasets analyzed, one region on chromosome 10 corresponding to the WDFY4 gene appears in both SE and nonSE (Figure 7, Table SI2 -5.3-I). It has been shown that WDFY4 is involved in immune response toward viral and tumor antigens 25 as well as in autoimmune diseases [26][27][28] . Functional annotation of variants with positive SDS scores 29 coupled with enrichment test 30,31 did not reveal any annotation category to be specific for a particular subset studied (S1:5.3, S1:5.4, Tables SI2-5.3-IV-V, Tables SI2-5.4-IV-VII). Likewise, alternative enrichment test employing the GWAS catalog showed that similar phenotype categories are present in the three tested datasets (S1:5.3, S1:5.4, Table SI2:5.3-III, Tables SI2:5.4-II-III).

On the other hand, frequency differences of rs75386033 derived allele T (10.3% in SE vs 6.1% in nonSE, Weir and Cockerham 32 (https://www.ncbi.nlm.nih.gov/gene/2911). These SNPs themselves are not known to be associated with any phenotypes, however, there are some indications that variant rs362870 which is in high linkage disequilibrium with rs75386033 and rs79907158, might be a cis-eQTL for the EPM2A gene (SI1:5, Table SI2:5-II), suggesting a plausible biological effect behind the frequency change. EPM2A gene is associated with Lafora disease which is a form of progressive myoclonus epilepsy [33][34][35] . This gene codes for a protein called laforin which is involved in regulating glycogen synthesis and potentially prevents glycogen accumulation in neurons [33][34][35] .

Given the lack of information on the phenotypic effect of this GRM1 allele and its modestly strong SDS signal, it is unclear whether the raise in frequency happened due to actual selection or because of random genetic drift especially given the fact that South-East Estonians exhibit signals of long-lasting low Ne and further differentiation into smaller subclusters (Figures 2 and6). Nevertheless, differential SDS signals between the entire Estonian dataset and its subsets including GRM1, WDFY4, suggest that recent selection, restricted to regional subpopulations, may remain undetected if population-wide datasets are treated as a single entity.

In conclusion, here we describe a dataset of more than 2300 high-coverage Estonian genomes from a population genetics perspective making it one of the smallest populations to date with such high-resolution data available. We show that the Estonian population, despite occupying a small area with no strong geographic barriers, is genetically structured and exhibits rather pronounced interregional differences with respect to recent admixture with neighbouring populations, population dynamics and potential action of natural selection.

These observations together with results of other studies suggest that population stratification could be ubiquitous in human populations, and should be taken into account in any large-scale genetic study including reconstructions of recent population history. We also show that we are able to accurately link effective population size to actual census size based on some simple assumptions about human population age structure and reproduction patterns. We envisage future studies exploiting this framework to ultimately unlock the potential of genomic data to provide a reliable estimate of past human census size, hence informing other historical sciences such as the study of cultural evolution, history and archaeology.

METHODS

Data reporting

No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Whole Genome Sequencing data

We used whole genome sequences of the Estonian Biobank participants reported in Kals et al. 10 . We applied exactly the same criteria for filtering individuals (sequencing quality control filters, match between WGS and chip genotype, total number of SNVs, self-reported Estonian ethnicity etc., see Kals et al. 10 for details) except for relatedness and singleton count (see below). For manipulations with vcf files bcftools-1.8 36 was used unless specified otherwise. Additionally to sample filtering applied by Kals et al. 10 , we removed seven samples with missing call over 3% as well as related individuals. To do so we used PLINK-1.9 37 and KING-2.1.6 38 to estimate relatedness coefficient and removed one individual from each pair with this value equal to 0.0442 or higher, corresponding to third degree relatives.

This resulted in a dataset consisting of 2,305 individuals that was used for all downstream analyses.

For analyses that require phased and/or imputed data (CHROMOPAINTER, SDS) phasing and imputation was done using Eagle v2.

CHROMOPAINTER/fineSTRUCTURE/GLOBETROTTER

To study genetic similarities between Estonians and other European populations we used the CHROMOPAINTER/fineSTRUCTURE (CP/FS) pipeline 17 . It involves a chromosome "painting" procedure which represents each chromosome of an individual (the recipient) as a mixture of chunks received (copied) from every other individual in the dataset (donor). The number of chunks copied by a recipient from each of the donors makes a "copying vector" which are used in the FS algorithm to group individuals into populations.

Initial chromosome painting parameters were estimated using 30% of the phased dataset of 1068 Estonian and non-Estonian samples (Table SI2-3.1-I). FS was run for 15 million MCMC iterations in two parallel runs to assess convergence. The tree-building step was performed using the approach from Leslie et al. [START_REF] Leslie | The fine-scale genetic structure of the British population[END_REF] and the run with the highest observed posterior likelihood was used to cluster the samples into genetic groups. Inferred FS groups were further manually inspected and clustered into the higher-order FS populations (S1:3.2).

These FS groups were used as surrogate populations to infer admixture with GLOBETROTTER and estimate ancestry profile with NNLS.

Next, GLOBETROTTER (GT) 17 was used to detect signals and dates of admixture for the Estonian groups defined using the approach described above. Unlike many other methods, GT allows the structure of unsampled source populations which were involved in the admixture event(s) to be assessed by modelling them as a mixture of sampled surrogate populations. GT inference was performed using a "regional" approach 17,42 . Estonian clusters were only allowed to copy from external surrogates, but not from other Estonians.

CHROMOPAINTER parameters were estimated for each Estonian target group individually and the average over all target populations was used to prepare input copying vectors for GT. Two separate runs, with and without standardization by "NULL" individual, were performed and consistency between runs was checked. To assess whether unbalanced surrogate population sample size could have biased our GT inference, we performed five additional GT runs by down sample both target and surrogate populations to 20 individuals.

Finally, given complex admixture signal in Estonia, we implemented non-negative leastsquares (NNLS) method [START_REF] Leslie | The fine-scale genetic structure of the British population[END_REF] . This allowed us to assign relative ancestral proportions to each individual in the R50+ panel using the non-Estonian surrogate groups identified by FS as sources. NNLS values for CP/FS Estonian groups were extracted from GT output while for individual samples these were calculated with an in-house R script. Obtained results were then summarized across Estonian parishes as well as across IBD/FS clusters.

Detecting segments identical-by-descent (IBD segments)

To

IBD segments between Estonians and non-Estonian individuals were detected by applying

refined IBD version 12Jul18.a0b 43 with default parameters except for length (window=40.0, length=1.0, trim=0.15, lod=3.0) to the same dataset that was used for CP/FS/GT, as in this case the dataset is highly structured. This was followed by applying the merge-ibd utility version 12Jul18.a0b to merge together segments separated by at most 1 cM long gaps and no more than 2 positions with genotypes discordant with IBD.

Both for IBDseq and refined IBD/ibd-merge results segments shorter than 2 cM were discarded, as longer segments are detected with higher reliability.

MAPS

In order to evaluate the extent of gene flow across the whole country together with local population densities, we estimated migration surfaces using MAPS 15 , which harnesses a similarity matrix summarizing the total number of IBD segments shared in a given population. In doing so, we used the IBD segments shared among pairs of individuals inferred with IBDseq as described in the previous section. Subsequently we have classified the shared genetic fragments as "short" (between 2 and 6 cM) and "long" (more than 6 cM), and performed two independent MAPS runs for each length class to assess convergence.

Estonian territory was modeled as having a total of 200 demes. Each run had 5 million iterations thinned every 10,000 and preceded by a burn-in of 2 million discarded cycles. The obtained migration surfaces were subsequently plotted using the plotmaps R package 15 . We repeated the whole procedure after removing samples belonging to clusters from Figure 2 with mean sharing above 60 cM to assess their effect on MAPS results.

IBD-based fineSTRUCTURE (IBD/FS)

In order to exploit patterns of genetic similarity between samples that arose very recently and get insight into fine genetic structure of the Estonian population, we used total genetic length of IBD segments longer than 2 cM as a measure of genetic similarity between pairs of individuals. We refer to this measure as "IBD-sharing". Next, to obtain natural genetic grouping of the samples we used a matrix of IBD-sharing as input for fineSTRUCTURE v2.0.7 14 . Although our approach is different from the original CHROMOPAINTER/fineSTRUCTURE method 14 , it is very similar in its idea to the approach used in 3 , and, put loosely, treats each cM shared between a pair of individuals as a CHROMOPAINTER chunk copied by the recipient from the donor. The fineSTRUCTURE algorithm already has an inbuilt method of compensating for the fact that the units used to measure similarity/relatedness between samples (either chunks in the classical approach or cM in our approach) don't represent fully independent pieces of information by transforming the raw value into an effective one by applying a c-factor. The c-factor was calculated using the fs combine command with the -C option applied to matrices of IBD-sharing for each individual chromosome. For more details see Supplementary Information SI1:2.1 and SI1:2.2. When running fineSTRUCTURE for both R50+ and the entire dataset the first 2,000,000 MCMC iterations were removed as burn-in and subsequently MCMC was run for additional 2,000,000 MCMC iterations sampling every 10,000th run. When building the tree we used the approach described in Leslie et al., 2015 1 and corresponding to the "1" value of the -T option in fineSTRUCTURE v2, which, put informally, maximizes the concordance between samples' final cluster assignment and its' assignment in individual MCMC runs. To validate this approach we applied it to the simulated data used in Lawson et al., 2012 14 , and calculated the same measure of correlation between real and inferred cluster assignment of the samples for different number of chromosomal regions used to detect IBD segments (SI1:2.2).

We applied this approach to the R50+ dataset (468 samples) and the entire dataset (2,305 samples). In both cases fineSTRUCTURE was run twice to assess convergence (SI1:2.3, Tables SI2-2.3-I and SI2-2.3-II). In the case of the R50+ dataset to reduce the number of clusters revealed by the fineSTRUCTURE algorithm we have hierarchically joined together clusters with short terminal branches by cutting the tree at such a level so as to avoid having clusters consisting of less than 5 samples. In the case of the entire dataset clusters referred to throughout the text were obtained by cutting the tree at a level chosen after visual inspection (SI1:2.3).

Fst calculations

Fst between Estonian clusters was calculated using smartpca from the EIGENSOFT package v7.2.0 41 after LD-pruning (r 2 > 0.4, windows of 1,000 SNPs) and removing sites with MAF < 0.05 and missing rate > 0.1.

Per-site Weir and Cockerham 32 Fst estimator between SE and nonSE subsets was calculated using vcfttools 40 after filtering sites for LD, MAF and missing rate the same way as described above.

Geographic data visualization

Geographic coordinates of the corresponding birth town/parish were assigned to each sample with birth place information available (2,168 out of 2,305 samples). For MAPS these coordinates were used directly. When plotting IBD/FS and NNLS results for the R50+ panel, coordinates of the samples were changed manually to avoid over-plotting. When plotting samples from the entire dataset jittering was used for the same purpose. Shp objects used to plot maps of Estonia with parish and county borders were retrieved from the Estonian reproducing. We therefore adapted a more general approach from 22 which incorporates the inbreeding coefficient (Fis), relative fraction of males (m) and excess in variance of reproductive success compared to the Poisson distribution (DV):

N b (t ) = (1+ Fis ) 4 ×( 1 (1-m)× m + DV )× N e (t )
[1]

Formula [1] yields the Nb, the number of breeding individuals (individuals capable of reproducing) at time t under the assumption of absence of gene flow and population structure, non-overlapping generations and equal variance of number of offspring between sexes. It is dependent on parameters such as m, Fis and DV that cannot be reliably estimated for each time bin. We therefore explored a range of plausible scenarios described by different values of m, Fis and DV based on the following assumptions: i) Fis calculated on chromosome 1 for contemporary human populations from the 1000 genomes dataset 53 as well as for Estonians ranges from -0.016 to 0.004 (S1:4.4), leading to the conclusion that for most human populations the term (1+Fis) can be safely approximated to 1; ii) m, the relative fraction of reproducing males, must be comprised between 0.1 and 0.9, considering further polarizations of this parameter as implausible for our species; iii) DV, the difference between the expected and the observed variance in number of offspring per adult can be estimated to range between -1 and 3. The latter estimate was obtained by taking Poisson distributions constrained between 0 and 10 (considering 10 as the maximum number of surviving children per adult) with an average between 1 and 5, and by empirically inflating the most extreme bins (0 or 10 children per adult) 5-fold. Such an exercise yields DVs ranging between -0.2 and 2.5, which we conservatively rounded to -1 and 3, respectively. This range is also consistent with data from contemporary Estonians, available from the Estonian Biobank, and showing a DV of -0.76 based on 7,863 females born between 1900 and 1955 and in the age of menopause at the time of enrolment. When plugged together into formula [1], these estimates yield a minimum of 0.75 (with m=0.5 and DV=-1) and a maximum of 3.53 (with m=0.1 or 0.9 and DV=3). To provide a single point estimate of Nc we rewrite formula [1] as [2] using a geometric mean between 0.75 and 3.5 and thus making our estimate slightly more than 2-fold away from the provided range boundaries. Note, that although there are indications that in some human populations DV can be higher than 3 54 , such cases can be considered to be at the very extreme of human reproductive behaviour spectrum as even hypothetical "super-male" populations would have a sex-average DV of 1.8 given m equals to 0.5 [START_REF] Heyer | Sex-specific demographic behaviours that shape human genomic variation[END_REF] . Hence we suggest our approach to be applicable to many human populations provided that immigration and population structure can be properly accounted for. In addition, the range of DV can be changed to study populations with extreme inequality in reproductive success.

N b (t ) =1.63 × N e(t )
The value estimated using [2] corresponds to the number of individuals in reproductive age.

It can be converted into total census size (Nc) of a human population at a given time point by dividing it by the estimated fraction of breeding individuals, which we here assume to be roughly 0.33. This is supported by actual data on the Estonian population from the "Statistics Estonia" database (http://andmebaas.stat.ee/Index.aspx?lang=en#) showing that the fraction of people between 20 and 40 years old was between 0.33 to 0.38 during the period between 1970 and 2018. Incorporating this idea into [2] results in equation [3].

N c(t ) =4.89 × N e(t) [3] which we used to obtain the curve in Figure 6B. Sources of historical estimates of Estonian population size used in that figure are provided in S4.2-II.

When using Ne as a proxy for actual population size one should keep in mind the potential effect of gene flow between populations. For example under a stepping stone model with constant population size and migration Ne estimated using samples from one deme is expected to increase when going back in time as more and more ancestors of sampled individuals would represent other demes [START_REF] Browning | Ancestry-specific recent effective population size in the Americas[END_REF] . In other words, coalescent-based Ne estimates reflect the number of ancestors of a sampled population, which may have lived in any location in space, rather than strictly the number of individuals in a given area at a given time point.

Singleton density score (SDS) selection scan

As SDS 20 does not handle missing data, imputed genomes of 2,301 unrelated individuals (four PCA outliers removed) were used. SDS 20 analysis was applied to three datasets separately, namely, the entire dataset and it's two subsets, Estonia SE and Estonia nonSE.

The latter two were defined based on the IBD/FS results (Figure 5): SE (individuals with South-East Estonian ancestry belonging to clusters eSE_1-eSE_5) and nonSE (individuals coming from the other parts of the country and belonging to other clusters). We did not apply SDS to finer subclusters of the dataset due to sample size issues. The number of individual genomes used in the analysis has a direct impact on the number of SNPs analyzed, the power to detect selection at any given SNP and the length of the terminal branches of the coalescent tree and hence the timing of the selection events that can be inferred 20 . Since in this study we are focusing on recent selection signals and differences between subpopulations that show only limited differentiation, we aimed at using samples of at least 1000 individuals. As the dataset needs to be homogeneous in terms of number of singletons per individual, this value was calculated with vcftools 0. the geometric mean between the two. In both panels on the y axis, "k" stands for "thousands" and "M" for millions. - a q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Bel Lat Lit cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM 10 cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM cM 20 30 40 a q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Haapsalu Jogeva Kohtla-Jarve c q q q q q q q q q q q q q qqqqq qq qqqq qq qq qq qqq qq qqq qqqq qqq qqq qqq qqq qqq qqq qqq qqq qqq q q qq q q q q q q q q q q q q q q 25 k 50 k 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0 1 1 0 0 1 2 0 0 1 3 0 0 1 4 0 0 1 5 0 0 1 6 0 0 1 7 0 0 1 8 0 0 1 9 0

Mor

  3).Based on our simulations and MAPS results, we propose that most of the differences in Ne dynamics between Estonian subpopulations may be attributed to different patterns of gene flow and external admixture. South-West and North-West Estonia are characterized by an of gene flow (SI1:2.4), leading to similar Ne trajectories that deviate only during the last 20 generations (Figure6, SI1:4.2) reflecting very recent differences in population size dynamics and/or migration. This also brings about the idea that the strong bottleneck in South-West could contribute to the observed population structure, in particular leading to differentiation of South-West and its subgroups. On the other hand, South-East Estonia has the most distinct Ne trajectory according to Figure6a, having a substantially lower long-term Ne compared to other regions. Together with MAPS results (SI1:2.4) this might suggest a recent expansion of a previously small-size eSE_5-like population and its admixture with other local subpopulations occupying South-East Estonia thus contributing to other eSE groups. This, in turn, results in a rather recent increase in relative proportion of individuals with eSE_5-like ancestry in the entire Estonian population affecting the Ne reconstructions for the entire dataset (SI1:4.2).

  Fst=0.0117 corresponding to the 0.999 percentile of the genome-wide distribution (Figure S5.4-I) together with its low standardized SDS p-value form strong evidence for a recent frequency increase of the rs75386033 T allele in South-East Estonia. Both rs75386033 and rs79907158 lie within an intron of the GRM1 gene,

  03 and positions in LD (r 2 > 0.4 within sliding windows of 200 positions). Results obtained in setting a were used to identify Estonian samples with extreme position in the PCA plot to be removed from the R50+ panel and from the dataset used for SDS (S1:1.2). In setting b we used 255,536 bi-allelic SNPs that overlap between the different datasets and passed LD-pruning (r 2 > 0.6 within sliding windows of 1,000 positions), MAF (<0.05) and no-call (>0.05) filters. We first calculated the principal components (PCs) based on all non-Estonian samples and then projected the Estonian individuals onto the first two PCs.

Figure 5 .

 5 Figure 5. Genetic clustering of the entire Estonian dataset(2,305 samples). Samples were clustered using the fineSTRUCTURE clustering algorithm based on pairwise total genetic length of IBD segments as described in Methods. Obtained clusters were pulled together based on their position on the tree resulting in 12 higher order clusters shown here (SI1:2.3).A: Hierarchical relationships (tree) and average total length of IBD segments shared between clusters (heatmap). The length of the tree branches does not reflect any relationship between the clusters. Numbers in grey next to cluster names show the number of samples in each cluster. B: Geography of inferred clusters. Each dot within the contour of Estonia corresponds to one individual, while waffle plots show samples for 15 major Estonian towns with each dot corresponding to 5 individuals. This map was created in R (https://www.R-project.org/)44 using an shp object of the Administrative and settlement units provided by the Estonian Land Board, 2018.11.01(https://geoportaal.maaamet.ee/eng/Spatial-Data/Administrative-and-Settlement-Division-p312.html). See Methods for more details.

Figure 6 .

 6 Figure 6. Estonian effective population size dynamics. a, Effective population size estimates obtained by applying IBDNe19 to the entire dataset and to 4 clusters from Figure 5: eNW_1, eNE, eSW_2 and eSE_5. b, Comparison of historical and genetic estimates of Estonian population size. Historical estimates combine census data and reconstructions based on written or archaeological sources (S1:4.2-II). Genetic estimates are derived from IBDNe results, for which Est1527 subset was used (SI1:4.4-II) and refer to the broader population that contributed over time to the genomes of contemporary Estonians. When converting time points of the IBDNe curve into actual years we used the same logic as in the original publication19 and set generation 0 to correspond to the year when individuals in our sample had a mean age of 25 (1988). Generation time of 29 years was assumed. For year 1200 the minimum and maximum estimates are provided. In panel a shaded areas show 95% confidence intervals. In panel b shaded area corresponds to the range between the minimum and maximum genetic estimates of Nc (Methods), while the light blue line shows

Figure 7 .

 7 Figure 7. Genome-wide plots of positive standardized SDS scores for the entire dataset (a) as well as SE (b) and nonSE (c) subsets. Conditional suggestive (blue) and genome-wide (red) significance lines are drawn. Gene names are highlighted for intragenic variants with -log10(p) > 5. Datasets are described in the text and Supplementary Information SI1:5.1.

  3 39 on the dataset consisting of 2,420 samples to benefit from the presence of related individuals and subsequently relevant samples were extracted. All Estonian Biobank participants have signed a broad informed consent which allows research in the fields of genetic epidemiology, disease risk factors and population history. All work at Estonian Biobank is conducted according to the Estonian Human Gene Research Act. The original study generating the WGS data 10 was approved by the Research Ethics Committee of the University of Tartu (application number 234/T-12).

	The second dataset used for CHROMOPAINTER/fineSTRUCTURE/GLOBETROTTER
	except for R50+ Estonians included 425 samples from 27 populations genotyped on Illumina
	arrays and 175 samples from seven 1000 Genome Project populations (CHB, FIN, GBR,
	GIH, IBS, TSI, YRI) (Table SI2-3.1-I). Whole genome sequences were pruned to keep
	positions matching those overlapping between genotyping arrays resulting in approximately
	500K SNPs.
	Principal component analysis
	The 'Rural above 50 years old' (R50+) panel We ran principal component analysis (PCA) for the entire Estonian dataset in two settings: a)
	Non-Estonian samples
	To place the Estonian population genetic variation in Eurasian context we compiled two
	datasets containing the R50+ Estonian samples each and samples from various populations
	predominantly representing West Eurasia. The first dataset used for PCA contained 59
	samples from 17 populations sequenced on the Complete Genomics platform, 207 samples
	from 8 populations sequenced using Illumina technology and 255 samples from 14
	populations genotyped on Illumina arrays (Table SI2-1.2-I). Whole genome sequences were
	pruned to keep positions matching those overlapping between genotyping arrays resulting in
	approximately 450K SNPs.

As information on parents' and grandparents' birthplace is mostly unavailable for the samples used here, we subsetted the 2,305 dataset for individuals born in rural areas and sampled at the age of 50 or older as we expect this cohort to be the least affected by recent migration and long-distance marriages, hence expecting it to preserve the original genetic structure. This resulted in a dataset of 474 individuals which we further pruned for PCA outliers (see below) and samples with more than 10,000 singletons (SI1:1.1-SI1:1.3). Persample number of singletons was estimated using vcftools-0.1.14 40 on the entire

(2,305

samples) non-imputed dataset. We ended up with a panel of 468 individuals, which we call "R50+". with only the 2,305 Estonians and b) combining the 2,305 Estonians with 521 non-Estonian samples from 18 European populations (Table

SI2

-1.2-I). In both cases smartPCA from EIGENSOFT-7.2.0 41 was used. In setting a we directly ran PCA on the dataset filtering for MAF below 0.01, no-call above 0.

  detect IBD segments in the Estonian dataset we applied IBDseq version r1206 (10) with

	default settings (errormax=0.001, errorprop=0.25, r2window=500, r2max=0.15, minalleles=2,
	lod=3.0) to the non-phased non-imputed dataset consisting of 2,305 Estonians. Choosing
	IBDseq over refined IBD 43 here is justified by working with samples coming from a relatively
	homogeneous population, which makes IBDseq frequency model applicable, while IBDseq
	has the advantage of not requiring phasing as well as having sequencing errors and rare
	alleles being explicitly accounted for. As IBDseq software reports only physical coordinates
	of a segment's start and end we interpolated segments' genetic length in cM using GRCh37
	recombination map (ftp://ftp.ncbi.nlm.nih.gov/hapmap/recombination/2011-01_phaseII_B37/)
	using R 44 . When working with the R50+ panel corresponding IBD segments were retrieved
	from the general output obtained on the 2,305 dataset. Homozygosity-by-descent segments
	were also inferred with IBDseq.

European Journal of Human Genetics. Published 25 July 2020. https://doi.org/10.1038/s41431-020-0699-4

Acknowledgements

This work was supported by the Estonian Research Council grants PRG243 (GH, RF, LS, LP, MM), PUT1339 (AK), PUTJD817 (MK), MOBTP108 (UV), IUT20-60 (AM, RM), IUT24-6 (AM, RM), PUT1660 (TE) and PRG184 (LM); by the European Union through the European Regional Development Fund Projects No. 2014-2020.4.01.16-0030 (FM, MM, VP), No. 2014-2020.4.01.15-0012 (AM, RM, TE, LM, GH, LS, MM), No. 2014-2020.4.01.16-0024 (DM, LP, VP), No. 2014-2020.4.01.16-0125 (RF), No. 2014-2020.4.01.16-0271 (RF) and 2014-2020.4.01.16-0125 (AM, RM, TE, LM); by NIH GIANT (AM, TE); by ERA-CVD grant

DATA AVAILABILITY

The sequencing data are available on demand. The procedure of applying for the access to the data can be found under the following link:

https://www.geenivaramu.ee/en/biobank.ee/data-access.

Land

Board website (Administrative and settlement units, 2018.11.01, https://geoportaal.maaamet.ee/eng/Spatial-Data/Administrative-and-Settlement-Division-p312.html). Geographic data were visualized in R 44 with the aid of the following packages:

sp 45,46 , sf 47 , rgdal 48 , rgeos 49 and ggplot2 50 .

IBDNe

In order to reconstruct recent Ne dynamics we used IBDNe version 07May18.6a4 19 with default settings (npairs=0, nits=1000, nboots=80, trimcm=0.

Genetic simulations

To simulate genetic data under various demographic scenarios to test the behavior of IBDNe we used mspms which is an ms-compatible version of msprime 51 . Commands used for simulation are provided in the Supplementary Information S1:4.1.

Estimating actual census size based on Ne

Several lines of evidence, based both on theoretical reasoning 52 and empirical comparisons 19 suggest that in industrial human societies census size (Nc) is roughly 3 fold the Ne assuming a panmictic and isolated population. However, application of this coefficient is limited to populations with specific reproductive characteristics, for

(MM) and PRESICE4Q (LM).

Data analysis was performed on the High-Performance Computing Center of University of Tartu.

The authors would like to thank the Genome Aggregation Database (gnomAD) and the groups that provided exome and genome variant data to this resource. A full list of contributing groups can be found at https://gnomad.broadinstitute.org/about.

The authors would like to thank Bayazit Yunusbayev for fruitful discussion and valuable advice.

Competing interests

Authors declare no competing interests.