Enrichment of low-frequency functional variants revealed by whole-genome sequencing of multiple isolated European populations

The genetic features of isolated populations can boost power in complex-trait association studies, and an in-depth understanding of how their genetic variation has been shaped by their demographic history can help leverage these advantageous characteristics. Here, we perform a comprehensive investigation using 3,059 newly generated low-depth whole-genome sequences from eight European isolates and two matched general populations, together with published data from the 1000 Genomes Project and UK10K. Sequencing data give deeper and richer insights into population demography and genetic characteristics than genotype-chip data, distinguishing related populations more effectively and allowing their functional variants to be studied more fully. We demonstrate relaxation of purifying selection in the isolates, leading to enrichment of rare and low-frequency functional variants, using novel statistics, DVxy and SVxy. We also develop an isolation-index (Isx) that predicts the overall level of such key genetic characteristics and can thus help guide population choice in future complex-trait association studies.

P opulation variation in disease susceptibility has been shaped by environment, demography and evolutionary history. Isolated populations (isolates) have generally experienced bottlenecks and strong genetic drift, so by chance some deleterious rare variants have increased in frequency while some neutral rare variation is lost, both helpful characteristics for the discovery of novel rare variant signals underpinning complex traits [1][2][3] . Studies to date have focused on individual isolates and have identified several disease-associated signals [4][5][6][7][8][9][10][11][12] . However, isolates differ in the time when they became isolated, their initial population size, the level of gene flow from outside and other historical demographic factors, and consequently also differ in their power for association studies 2 . We thus generate and analyse low-depth (4 Â -10 Â ) whole-genome sequences (WGS) from eight cohorts drawn from isolated European populations and compare each isolate with the closest non-isolated (general) population, for which we also generate or access WGS data. We then investigate empirically how these historical differences influence the population-genetic properties of isolates, and frame these insights in terms of their consequences for study design in complex trait association studies.

Results
Samples, sequencing and QC. The data set includes newly generated low-depth (4x-10x) WGS from eight cohorts drawn from isolated European populations: one each from Kuusamo in Finland (FIK) and Crete in Greece (GRM 13 ), four from Friuli-Venezia Giulia in Italy (IF1, IF2, IF3 and IF4 (ref. 14)), and one each from Val Borbera in Italy (IVB 15 ) and the Orkney Islands in the UK (UKO 16 ); and the closest non-isolated (general) population: Finland (FIG 9 ), Greece (GRG), together with publicly available data for Italy (ITG 17 ) and UK (UKG 18 ) ( Fig. 1a and Supplementary Table 1). We generated a superset of variants called in these cohorts and all 26 population samples in the 1000 Genomes Project Phase 3 (ref. 17), and performed multi-sample genotype calling across all 9,375 samples (3,059 from the current study, 2,353 from the 1000 Genomes Project Phase 3 release and 3,781 from UK10K). Both individual population and amalgamated genotype call data, which have greater than 99% concordance with genotyping data (Supplementary Table 2), are available to the scientific community (Data availability).
General description of the variants in the isolates. We identified approximately 12.2 million variants with minor allele frequency (MAF) r2% (rare), 5.5 million with MAF 42-r5% (low-frequency) and 8.3 million variants with MAF 45% (common) across the ten populations newly sequenced here (eight isolates, GRG and FIG). Of these, 10.5, 0.7 and 0.3%, respectively, are novel (Table 1 and Supplementary Table 3). As expected, most of the isolates have lower numbers of variant sites per genome than their closest general population ( Supplementary Fig.  1, Supplementary Table 5). We find B188,000-B513,000 variants that are common with MAF 45.6% in each isolate but with MAFr1.4% in its closest general population (Table 1); B30,000-122,000 of these per isolate have frequency r1.4% in all the general samples studied, among which B150-B700 in coding regions and B500-B2,800 genome-wide are deleterious (Supplementary Table 4). These common and low-frequency variants are thus useful markers for whole-genome association studies in these populations and some of them (if absent from the general population) could potentially lead to novel association signals. They include known examples such as rs76353203 (R19X) in APOC3 in GRM, which is associated with high-density lipoprotein and triglyceride levels 6 .
Population-genetic analyses in the isolates. Previous population-genetic studies of isolates have, with some exceptions 11,19 , been based on common variants found on genotyping arrays, and have illustrated general characteristics such as low genetic diversity and longer shared haplotypes 9,[13][14][15]19,20 . Rare variants discovered from sequencing are on average more recent in origin than common variants 21 and therefore more powerful for distinguishing closely related populations and more informative about recent demographic history. We find that isolates are, as expected, genetically close to their matched general population in principal component analyses (PCA), ADMIXTURE 22 and TreeMix 23 using common variants (Fig. 1b, Supplementary  Figs 2-5 and Supplementary Table 6), but PCA using rare and low-frequency variants, as found previously 24 , distinguishes them more clearly from the general population and also from other isolates, particularly among the Italian samples (Fig. 1c, Supplementary Fig. 2). The majority of sharing of variants present just twice across all samples of 36 individuals from each population (f 2 variants 21 ) takes place within the same population, and the isolates generally share more with their closest general population than with other populations. This latter trend, however, is not apparent for IF1-IF4, who show little sharing with any other population, pointing to a greater level of isolation and lower level of gene flow with their general population (Fig. 1d, upper triangle and Supplementary Fig. 7), which is confirmed by f3-statistics 25 comparing with a worldwide population panel of HGDP-CEPH samples using common SNPs ( Supplementary Fig. 6). f 3 -f 10 variant sharing demonstrates sharing by ITG and IVB with both Greek and UK populations (Fig. 1d, lower triangle and Supplementary Fig. 7), potentially indicative of their more ancient heritage.
Population demographic history. All populations studied here, both isolates and general, appear to have shared a comparable effective population size (Ne) history before 20 thousand years ago (KYA) based on the multiple sequentially Markovian coalescent method 26 (Supplementary Fig. 9). The isolates diverged from their general populations within the last B5,000 years based on LD estimations 27 (Supplementary Table 7 and Supplementary Fig. 8) and yet had sharp decreases in their population sizes in more recent times as estimated using inferred long segments of identity by descent (IBD) 28 (Fig. 1e,f and Supplementary  Fig. 10). Different isolates also split from their respective general populations at different times. For example, IF1-IF4 split from ITG B4-5 KYA, while most other isolates split from their general populations within the last B1,000 years (Supplementary Table 7).
The different demographic histories of different isolates should lead to different genetic characteristics. To summarize these features in a single quantitative measure that can be calculated from genotype data, as well as sequence data, we developed an isolation index (Isx) which combines information on the divergence time from the general population (Tdg), Ne and migration rate (M), such that early-divergence-time isolates with small Ne and low M have a high Isx value ( Fig. 2a and Supplementary Fig. 11 Fig. 12), the total length and number of runs of homozygosity (ROH) (Supplementary Fig. 13), inbreeding coefficient (F) (Supplementary Fig. 14) and length of LD (Supplementary Figs 15 and 16 and Supplementary Tables 9 and  10). All these characteristics are correlated, but the pairwise correlation coefficients show that Isx is a slightly better overall predictor of the other measures than any single existing measure ( Fig. 2c, Supplementary Fig. 17 and Supplementary Table 11); moreover, it is potentially more robust to confounding factors as it is calculated from three demographic parameters, while the others are all based on single measurements.   IF1  IF2  IF3  IF4  IVB ITG UKO   Purifying selection analyses. Several lines of evidence suggest relaxed purifying selection in the isolates due to their reduced Ne, although as expected we do not detect substantially increased genetic load per genome using the Rxy statistic 29 based on all of the variants in the genomes ( Fig. 3a and Supplementary  Table 12). First, we see different levels of enrichment of lowfrequency functional variants in isolates (Fig. 3b,c, Supplementary  Tables 13 and 14, Supplementary Fig. 18a) quantified by a new statistic, DVxy-coding, developed here (DV: drifted variants). DVxy-coding measures the ratio of functional coding variants (missense plus loss-of-function (LoF)) in isolates compared to the closest general population (and vice-versa), adjusted for the corresponding ratios of intergenic variants in order to correct for the effect of genetic drift. We applied this only to a subclass of DVs, defined as low-frequency (2-5%, the best choice according to the sample size we have) in any isolate, yet at least three-fold higher than in the closest general population (and vice versa). We find that DVxy-coding is 41 in all isolates and o1 in all general populations (Fig. 3c, Supplementary Fig. 18a and Supplementary  Table 13). We also calculated a similar DVxy-wg statistic by stratifying whole-genome variants according to their combined annotation dependent depletion (CADD) score (0-5, neutral variants; 5-10, mildly deleterious; 10-20, deleterious; and 420, highly deleterious; these cut-off choices balance the number of variants in each bin to allow us comparable statistical power among all bins, although the conclusions are robust to the particular cut-off values chosen and different bins ( Supplementary Figs 18b and 19)). The DVxy-wg values are differentiated for variants with CADD score of 10-20 and significantly so (assessed using the jack-knife bootstrap method) for ones with CADD scores 420, with DVxy-wg values 41 in all isolates and o1 in all general populations (Fig. 3b, Supplementary Fig. 18b and Supplementary Table 14). This demonstrates enrichment of low-frequency functional variants, both coding and genome-wide with CADD score 410, in the isolated populations. Moreover, both DVxy-coding and DVxy-wg values are correlated with Isx, suggesting that different isolation characteristics lead to different levels of enrichment of functional variants.
We also investigated the relaxation of purifying selection by assessing functional (missense) singleton variants (SV) pooled for all of the genes that have at least one singleton missense or synonymous variant in a pair of populations (one isolate and its general population), correcting with pooled synonymous variants (SVxy statistic,). We find a substantial deviation from 1 for functional singletons in all of the isolates (Fig. 3d and Supplementary Table 15), with SVxy values positively correlating with Isx ( Fig. 2c and Supplementary Fig. 20). We also find that the proportion of relaxed essential genes 30 with SVxy 41 in isolates is significantly higher than in the general population (Supplementary Table 15). Such rare and low-frequency drifted functional variants, measured by both SVxy and DVxy, are particularly relevant for boosting the power of association studies 6 .
Positive selection analyses. We do not find convincing evidence for positive selection in any isolate using deltaDAF 31 , PCAdapt 32 or singleton density score (SDS) 33 , although we do identify some highly differentiated variants (Supplementary Fig. 21 and Supplementary Tables 16 and 17), including in the proteincoding genes ALK, SPNS2, SLC39A11 and ACSS2, which can nevertheless be accounted for by drift. Interestingly, we also find six highly differentiated variants shared between different isolates from Italy, IF2, IF3 and IF4, but interpret them as likely to result from drift or positive selection for the ancestral allele in the ITG (Supplementary Table 17). We find that the SDS method has little power in our samples because of their small size, and failed to detect selection even at the lactose tolerance SNP in the UKO, a known strong signal of recent selection ( Supplementary Fig. 22).

Discussion
Isolated populations have special characteristics that can be leveraged to increase the power of association studies, as several previous studies have shown 19,34 . Nevertheless, only a small proportion of functional variants have increased in frequency in any one isolate, so multiple isolates must be investigated to reveal the full diversity of associated variants. Here, we probed an extended allele frequency spectrum of variants potentially underpinning human complex disease through the analysis of WGS data in multiple isolates matched to nearby non-isolated populations, capturing common, low-frequency and rare variants. We quantified different levels of isolation resulting from different demographic histories and have demonstrated that the Isx statistic, calculated even from SNP-chip data, reliably captures these relevant features. This study provides a systematic evaluation of the genetic characteristics of multiple European isolates and for the first time empirically demonstrates enrichment of rare functional variants across multiple isolates. With the advent of large-scale whole-genome sequencing, studies in isolates are poised to continue as major contributors to our understanding of complex disease etiology.   The sample sizes are very different across these collections, and we used three different standard-sized subsets of the samples for different analyses: (1) the whole data set; (2) the sample-size-matched data set, obtained either by randomly selecting samples from general population to match the isolated population (for example, we randomly select 377 from FIG to match FIK), or by randomly selecting a subset of the isolated population to match the general population (for example, we randomly select 108 IVB to match the general population ITG); (3) the minimum-sample-size data set of 36 individuals per population. By doing this, we maximize the use of the data for different analyses, and we specify which data set is used for each analysis. The sequencing depth is also different across different populations, within a 2.5-fold range (apart from GRG, in which variants were called differently, details in Supplementary Notes), and we allowed for these differences when interpreting the results.
Variant counts. We first re-annotated all variants using the Variant Effect Predictor annotation from Ensembl 76 with the '-pick' option, which gives one annotation per variant. We then performed variant counting at both the population and individual level, stratifying by functional categories and frequency bins. These counts were either plotted in figures or summarized as median values in tables. We carried out these analyses using both the sample-size-matched data set and the minimum-sample-size data set.
Population-genetic analyses. We used the whole data set for the analyses in this section, unless otherwise specified. PCAs were performed separately with common variants or rare variants using EIGENSTRAT v.501 (ref. 36). Shared ancestry between the populations studied here was evaluated using ADMIXTURE v1.22 (ref. 22). The relationships between the populations studied here, combined with worldwide populations from the HGDP-CEPH panel 37 , were also examined using ancestry graph analyses implemented in TreeMix v. 1.12 (ref. 23). We also used formal test of f3-statistics 25 to investigate population mixture in the history of the populations studied here, as well as worldwide populations from the HGDP-CEPH panel. Rare f 2 variants (with only two copies of the alternative allele in the minimum-sample-size data set) and moderately rare f 3-10 variants (3-10 copies of the alternative allele in the same data set) are particularly informative for investigating recent human history 21 . We investigated the sharing pattern of these two types of variant by summing all f 2 variants or any random two alleles of the f 3-10 variants shared by pairs of individuals. We plotted the results as a heat map using the image 1 function from the base R package (https://stat.ethz.ch/R-manual/ R-devel/library/graphics/html/image.html). Variants were aggregated by pair of individuals using the 'count' function of the plyr package, then arranged in matrix form and colourized using 'colorRampPalette' from the colorspace package (https://cran.r-project.org/web/packages/colorspace/index.html). ROH, inbreeding coefficient (F) as well as the length of LD-blocks were calculated in PLINK, and finally genome-wide F ST values between isolates and their general populations were calculated with the software 4P (ref. 38) using the minimum-sample-size data set.
Demographic inferences. LD-based 39-41 demographic inference was performed in the NeON R package 27 using the minimum-sample-size data set; the median and confidence interval were estimated using the 50th, 5th and 95th percentiles of the distribution of long-term Ne in each time interval. We used the multiple sequentially Markovian coalescent method 26 to infer demographic changes before 20,000 years ago using four individual sequences from each population. In order to account for some loss of heterozygous sites in the low-depth data, we used a slow mutation rate of 0.8 Â 10 À 8 mutations per nucleotide per generation and a longer generation time of 33 years. We then estimated more recent demographic changes (from the present to B9,000 years ago) using IBDNe 28 with the minimum-samplesize data set. We used IBDseq 42 to detect IBD segments in sequence data from chromosome 2 in all populations. We then used IBDNe with the default parameters and a minimum IBD segment length of 2 centiMorgan (cM) units. We assumed a generation time of 29 years.
Isolation index. In order to quantify the different isolation levels of different isolates, we developed an index that combines three demographic parameters: (a) Tdg, (b) Ne and (c) the level of private isolate ancestry (M). We call this estimate the Isolation index (Isx). It is defined as: Both Tdg and Ne were inferred from the LD-based method using the NeON R package 27 . M is difficult to estimate directly from SNP genotype data, so here we estimated the difference of shared ancestral components between an isolate and its general population from ADMIXTURE analysis. We ran ADMIXTURE with only one isolate and it closest general population using K ¼ 2. We then estimated the Rxy analysis. Rxy statistics 29 between each pair of populations (an isolate and its closest general population) for different functional categories were calculated using the matched-sample-size data for missense and LoF variants, including stop gain, splice donor and acceptor variants, using synonymous variants as controls (we did not use intragenic variants as control because of the ascertainment in the ITG which has high-depth exome sequences and low depth for the rest of the genome). We also calculated Rxy statistics for variants with CADD scores 43 greater than 10 and 20, using variants with CADD scores less than 5 as controls. The mean and s.d. for each Rxy value were obtained from 100 bootstraps.
DVxy analysis. A new statistic, DVxy, was developed to quantify the enrichment of low-frequency functional variants in the isolates using both the matched-samplesize and minimum-sample-size data sets. It calculates the proportion of functional variants in each isolate compared with its general population, correcting for genetic drift at the same time. We calculated DVxy specifically for the subset of variants with DAF 2-5% in the isolate, and at least three times lower in its closest general population, or vice-versa. We called these variants 'drifted variants' (DV). DVxy was calculated for both coding regions and whole genomes. For coding variants, we defined missense or missense plus LoF variants as functional variants. We counted the number of functional DVs and neutral (intergenic) DVs in each isolate (population x) and the corresponding general population (population y). The ratio between the fraction of DV variants from the isolated population (corrected by the count of intergenic variants) and the corresponding fraction of DV variants from its general population was defined as the DVxy statistic. If DVxy is equal to 1, there is no enrichment for the functional DVs in the isolate; less than 1 indicates depletion, and greater than 1 indicates enrichment.
DVxy coding¼ % DVx missense % DVx intergenic % DVy missense % DVy intergenic 0 For the whole genome, we used different CADD score cut-offs and bins. We calculated a DV statistic by stratifying the variants according to their CADD scores (0-5, neutral variants; 5-10, mildly deleterious; 10-20, deleterious; and greater than 20, highly deleterious) for each isolate and its closest general population. We finally calculated a ratio of the fraction of DV variants (from each class) between the isolate and its general population, and vice-versa. The following formula shows the DVxy-wg calculation for variants with CADD score between i and j in an isolate and its general population.
The 95% confidence interval for each calculation was obtained by randomly sampling data from 20 chromosomes 100 times.
SVxy analysis. We further investigated the relaxation of purifying selection in the isolated populations using SVs. Here, we also used the minimum-sample-size data set. Another new statistic, SVxy, was developed to measure the ratio of missense versus synonymous singletons per gene in each population, as well as the ratio of the sum of singletons in all genes which have at least one singleton in the pair of the populations (one isolate and one general population). We counted the number of missense singletons and synonymous singletons per gene in each population, and SVgene was calculated as: SVgene¼ SV missense count þ 1 ð Þ SV synonymous count þ 1 ð Þ SVgene41 indicates relaxation of purifying selection; SVgene ¼ 1 indicates neutrality; and SVgeneo1 indicates purifying selection. We then divided the gene list into essential genes 30 and non-essential genes (the rest), and calculated a statistic, G SV , for each population, defined as: G SV ¼ percentage of essential genes with SVgene41/percentage of non-essential genes with SVgene41 We finally calculated a statistic, SVxy, which is the ratio of SVpop of each isolate to SVpop of its general population. SVpop for each isolate and its general population was calculated using all genes which have at least one singleton in the pair of the populations and defined as SVpop ¼ S (SV missense counts)/S(SV synonymous counts).
We used the same annotation as in the variant counts. We calculated a confidence interval for each estimate using bootstrapping of 80% of the genes 100 times.
Correlation analyses. We calculated pair-wise correlation coefficients between the Isx values, population-genetic measurements ROH, F, F ST , and number and length of LD blocks, as well as the newly developed statistics DVxy and SVxy using the Pearson correlation in R.
Positive selection analyses. We calculated genome-wide pairwise derived allele frequency differences (deltaDAF) for each pair of populations (an isolate and its general population) as described previously 31 using the matched-sample-size data set. We also carried out PCAdapt analyses 32 for each pair of populations using the whole data set. Both analyses look for high derived allele frequency variants in the isolates, and will not be affected by sample size. Finally, we ran the SDS method 33 using the whole UKO and UKG data sets, which have the largest sample sizes for both isolate and its general population, and thus the greatest power for this method.
Data availability. Amalgamated genotype calls across all populations studied are available through the European Genome/Phenome Archive (EGAD00001002014) with Data Access Agreement described in Supplementary Information.