If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Institute of Genomics, University of Tartu, Tartu, 51010, Estonia
Department of Neurosciences “Rita Levi Montalcini,” University of Turin, Torino 10126, Italy
Institute of Genomics, University of Tartu, Tartu, 51010, Estonia
Department of Epidemiology, University of Groningen, Groningen 9700 RB, the Netherlands
Institute of Genomics, University of Tartu, Tartu, 51010, Estonia
Research Department of Genetics, Evolution and Environment, University College London, London WC1E 6BT, UK
Department of Biology, University of Padova, Padova 35131, Italy
Institute of Neurosciences, National Research Council (CNR), Padova 35121, Italy
Faculty of Health and Medical Sciences, University of Surrey, Guildford GU2 7YH, UK
Ancient groups differentially contributed to complex traits in contemporary Europeans
•
In contemporary Estonians 11 out of 27 traits show association with some ancestry
•
Hunter-Gatherer and Yamnaya ancestries divergently influence cholesterol levels
•
Post-admixture selection is not necessary to have trait-ancestry associations
Summary
The contemporary European genetic makeup formed in the last 8,000 years when local Western Hunter-Gatherers (WHGs) mixed with incoming Anatolian Neolithic farmers and Pontic Steppe pastoralists.
This encounter combined genetic variants with distinct evolutionary histories and, together with new environmental challenges faced by the post-Neolithic Europeans, unlocked novel adaptations.
However, how ancient populations contributed to present-day phenotypic variation is poorly understood. Here, we investigate how the unique tiling of genetic variants inherited from different ancestral components drives the complex traits landscape of contemporary Europeans and quantify selection patterns associated with these components. Using matching individual-level genotype and phenotype data for 27 traits in the Estonian biobank
and genotype data directly from the ancient source populations, we quantify the contributions from each ancestry to present-day phenotypic variation in each complex trait. We find substantial differences in ancestry for eye and hair color, body mass index, waist/hip circumferences, and their ratio, height, cholesterol levels, caffeine intake, heart rate, and age at menarche. Furthermore, we find evidence for recent positive selection linked to four of these traits and, in addition, sleep patterns and blood pressure. Our results show that these ancient components were differentiated enough to contribute ancestry-specific signatures to the complex trait variability displayed by contemporary Europeans.
EstBB contains matching genotype and phenotype information for individuals from a relatively homogeneous population, which contains all three ancestry components found in Europe, with the proportion of remnant Hunter-Gatherer ancestry among the highest in Europe and an additional minor (<5%) Siberian component associated with Iron Age movements.
In order to associate specific ancient European ancestry components with predicted phenotypes, we introduce covA, a measure of the relative similarity between any contemporary individual and the ancestries that contributed to its genetic makeup. For each sample in the EstBB, we computed its covA with each of the ancestral source populations, focusing on genomic regions potentially connected to each trait. We then used covA as a predictor to model traits, also in comparison with the same statistic computed for the whole genome. Finally, we test if those regions associated with the genetic contribution from a specific ancestry experienced a post-admixture selective pressure on top of the observed local unbalance in contributing ancestries.
covA measures similarity with ancestral groups
Here we introduce covA, the covariance between allele frequency (p) in a contemporary individual i (i.e., its allele dosage) and a given ancestral population j with respect to the contemporary and ancient average frequencies (
and
respectively):
(1)
The covA statistic is expected to be high when the allele frequencies of the individual i and the ancestry j are similar in comparison with the differences within the contemporary population and across the ancestries that contributed to its genetic makeup. Furthermore, covA can be computed averaging over the contribution of multiple single nucleotide polymorphisms (SNPs), either across the whole genome or for specific regions of interest.
In order to test the potential of covA to distinguish between genetic contributions from different ancestries, we simulated polygenic traits in a modern population composed of three ancestral groups and verified that when predicting simulated traits, covA estimated coefficient correlates well with their ancestral specificity (Pearson’s correlation coefficient ρ = 0.919-0.937, Figure S1A). See STAR Methods for further discussion of
properties and simulation details, including the definition of ancestral specificity.
For each EstBB individual, we computed genome-wide covA between the individual and each of the ancestries among Western Hunter-Gatherers (WHG), Neolithic Farmers from Anatolia (Anatolia_N), Yamnaya Pastoralists from the Pontic Steppes (Yamnaya), and Siberian (Siberia). We defined these ancestry groups based on genetic and chronological proximity to a set of identified focal individuals (see STAR Methods and Table S1 for a list of the ancient genomes assigned to each group). As expected, covAs calculated on the different ancestries are strongly interdependent because they include as term the average ancestral frequency
and because of varying grades of similarity among the ancestries for historical demographic reasons (see covA joint distributions in Figures S1C–S1F). In particular, covA tends to be negatively correlated between different ancestral components with the exception of Yamnaya and WHG, reflecting complex demographic relationships between the two, involving WHG-like Eastern Hunter Gatherer ancestry presence in Yamnaya.
Furthermore, although the Estonian population is considered relatively genetically uniform, some geographic differences exist with the southeastern inland counties having higher haplotype sharing with Latvians, Lithuanians, and Russians compared with the rest of the country, as recently shown in Pankratov et al., 2020.
This result is also mirrored in our analyses with median covA for WHG being higher in southeastern inland counties (see Figure S1G). Conversely, as shown by median covA for the Siberian component in Figure S1I, the Siberian ancestry seems to be more abundant in northeast Estonia, consistently with Finnish ancestry shown by Pankratov and colleagues.
Yamnaya and Anatolia_N covAs are instead more evenly distributed (Figures S1H and S1J).
Phenotype-associated genomic regions show specific ancestry similarity patterns
We examined 27 complex traits (31 if considering separate classes of pigmentation) for which we had sufficient records in the Estonian Biobank (see Table S2). We corrected and adjusted them for confounding covariates, including sex, age, genotyping platform, and others as specified in Table S2. Because our analysis relies on SNPs overlapping between ancient and contemporary genetic data, a portion of the genetic influence over these traits—especially when conveyed by rare alleles—might elude our experimental setting.
All traits analyzed and their estimated heritability after covariate adjustment. Bars indicate standard errors of the estimate. Numbers in parentheses indicate the number of unrelated samples for which phenotypic information was available for each trait.
See Figure S2 for partial heritability estimates on candidate regions.
hits for corresponding trait categories (see STAR Methods and Table S3). As shown in Figure S2, these genomic regions harbor a higher heritability intensity (
Mb) than the whole genome, supporting their suitability as candidate regions for the traits of interest.
Next, we used covAs computed on the candidate regions as a predictor to model traits and asked whether they showed significantly different regression coefficients when compared to 50 size-matching random genomic sets; this was found true in 11 out of 27 traits (double-sided Z-test, Benjamini-Hochberg false discovery rate (FDR) = 0.05), see Z-scores in Figure 2. This analysis has the advantage of automatically controlling for virtually all potential confounders that apply to the genome in its entirety—e.g., social, economic, and cultural statuses—thus allowing us to not include any such covariates in the model. In addition, this analysis pinpoints genetic signals that are likely to be functionally connected to the trait. Among others, blood cholesterol levels are shown to be positively correlated with similarity to Yamnaya in cholesterol-associated regions with respect to the rest of the genome, while the opposite is true for WHG.
Figure 2Ancestry-trait association on candidate regions
(A) Z-scores of covA coefficients, the color refers to the ancestry tested.
(B) Z-scores of coefficients associated with covA independent components (IC) computed with whole genome-based covA PC loadings. Each color is associated with one of the three ICs. For each trait we show the Z-score of the standardized coefficient associated with candidate regions against a distribution of 50 random genomic regions of matching size. Candidate regions are determined around GWAS hits for appropriate traits as windows with three different widths: 5 (small dot), 50 (medium dot), and 500 (large dot) kilobases. Pastel dots are deemed not significant at Benjamini-Hochberg FDR = 0.05, p value from double-sided Z-test; asterisks mark traits to be considered significant according to (B); dotted lines correspond to absolute Z-scores
.
(C) Loading matrix for genome-wide covAs and their PCs, used to transform covAs into their ICs. The three genome-wide (GW) PCs accounted for 0.498, 0.327, and 0.175 covAs variance, respectively. PCs and relative ICs can be interpreted as axes defined by 2 or 3 covAs.
See Table S4 for an interpretation of significant cases.
Since covA exhibits a high correlation across ancestries, we avoided implementing a model with largely multicollinear predictors including covA for all ancestries and instead adopted separate models for each ancestry, complementing them with a regression on covA independent components (ICs) (Figure 2B). We used the loadings from a principal component (PC) analysis on whole genome covAs (Figure 2C) to transform region-specific covAs into ICs. This, though not returning actual PCs in each candidate region, drastically reduces the collinearity (highest variance inflation factor = 1.62 in hair color 50kb candidate regions), while allowing simpler interpretation and, crucially, cross-region comparisons required for Z-scores computation. While covAs (Figure 2A) highlight the overall excess or lack of certain ancestries in relation with a given phenotype but are largely intertwined, ICs (Figure 2B) can be interpreted as independent axes defined by 2 or 3 covAs. We therefore adopted ICs to discriminate significant ancestry-trait associations, as they are independent variables in a comprehensive predictive model. Significant results, interpreted in light of the ICs, are summarized in Table S4 and discussed below. Among others, this analysis confirms the association between cholesterol levels and the Yamnaya-WHG axis previously mentioned.
Comparison with genome-wide ancestry similarity
We followed up the association between phenotypes and local excess or lack of a given ancestry and explored whether a similar pattern held at whole genome level by computing genome-wide covAs. Here, being unable to correct for environmental confounders with a Z-score approach and avoiding genotype-based PCs as covariates in order not to hinder potential genome-wide signals, we run the risk of obtaining spurious ancestry-trait associations. This is due to uneven ancestry similarity across Estonia concurrent with geographically associated socio-economic differences that can potentially confound genotype-phenotype associations. Although the confounding effect of population structure is minimized by the inclusion of a relatively uniform population, small differences related to historical reasons
are still visible in covA (see Figures S1G–S1J). Therefore, we include a city/countryside residency covariate in the models, defined as 1 for people living in the wealthiest and most populous county (Harju county) and 0 otherwise, and a covariate for educational attainment, which is a good proxy for family socioeconomic status.
This control allows us to suggest a significant influence of genome-wide ancestry on 16 traits out of 27, as shown in Figure S3, even when geographical and social stratification is present (coefficient p value significant at Benjamini-Hochberg FDR = 0.05). Again, covA-based PCs were used to interpret significant results.
Interestingly, we do not always observe concordance between the region-specific and genome-wide results, as shown in Figure 3, pointing to the fact that region specific trends are not always sufficient to drive genome-wide signals to significance or might even arise in a contrasting genomic background. This is especially true for less polygenic traits (e.g., pigmentation) but also for more polygenic ones, as indicated by height association with WHG. On the other hand, we also find genome-wide ancestry-trait connections that are not exacerbated in candidate regions, thus losing Z-score significance. This can occur for a single ancestry (e.g., Anatolia_N or Siberia and height) or cause the loss of trait associations altogether, as for alcohol consumption, depression, sleep duration, social jetlag, diopters, pulse pressure, and creatinine levels. Finally, we observed that genome-wide covAs for WHG and Yamnaya tend to be linked to most phenotypes in a similar fashion, in contrast with results found in candidate regions where the two ancestries behave in a more independent manner (Figure S3).
Figure 3Comparison between genome-wide and region-specific ancestry-trait associations
(A and B) The y axis represent Z-scores of covA coefficients, for covA computed on candidate regions (CR) of 5 kilobases as in Figure 2. X-axes represent genome-wide (GW) covA estimated coefficients: we report
effect sizes for continuous traits in (A) and Odds Ratios for categorical traits in (B). Independent models are run for different covAs. Colors label the ancestry tested, while inner and outer color intensity represents significance of CR covA Z-score and GW covA coefficients, respectively. Pastel colors indicate not significant results at Benjamini-Hochberg FDR = 0.05 (double-sided Z-test p value for CR covA Z-score or double-sided coefficient p value for GW covA coefficients). Labels indicate selected outlying ancestry-trait associations.
GW covA coefficients are also reported in Figure S3.
Selection signatures at candidate regions with ancestry-trait association
So far, we only explored associations between a given trait and a local or genome-wide excess of a given ancestry. The observed local admixture imbalance points to a role of that ancient contribution in explaining a given phenotype. However, these results alone do not show whether after the admixture event the incoming genetic material also underwent a selective sweep within the recipient population, altering population-wide allele frequencies as investigated in Mathieson et al., 2015.
to the list of GWAS hits used as index for our candidate regions to obtain per-SNP evidence of recent (up to 500 generations ago) natural selection and to see which phenotypes show enrichment in SNPs with strong selection signals compared to a random set of GWAS hits. Out of the genomic regions responsible for ancestry-trait association shown in Figure 2, pigmentation-related SNPs (eye and hair color) showed extremely high CLUES logLR values (Figure 4A; Figure S4) in accordance with previous results,
as well as SNPs related to BMI and cholesterol, pointing to ongoing or recent selection at these loci. Diastolic blood pressure (DBP) and sleep-related SNPs also showed the same extreme signature, but the candidate regions encompassing them did not reach significance in ancestry-trait association.
(A) CLUES log likelihood ratios (logLR) values distribution for GWAS hits for six selected phenotypes. For each phenotype, at most 100 top SNPs with highest logLR values and the corresponding ranks from the random GWAS hits distribution are shown. Gray dots show mean values for each rank in the background distribution while the whiskers show the 5-95 percentile range. The logLR values for tested SNPs are shown in red or blue depending on whether the value lies above the 95th percentile of the values from the background distribution with a given rank. Number of tested SNPs for each phenotype are shown in panel titles. Sleep indicates SNPs connected to all sleep traits as indicated in Table S3. See also Figure S4.
(B) Maximum likelihood estimates of derived allele frequency trajectories for top three SNPs with highest logLR values for each phenotype. When more than one SNPs come from the same locus, only the top-scoring SNP is shown.
The recent and putatively ongoing nature of the inferred selective pressure on the six traits shown in Figure 4A is further exemplified by the steep increase in derived allele frequencies over time inferred for the top three SNPs of each trait and shown in Figure 4B. These include some loci previously shown to be selected in West Eurasians (rs4988235 at MCM6/LCT,
) and some others yet to be explored. In particular, rs17630235, associated with BMI and DBP, is an expression quantitative trait locus (QTL) in several epithelial tissues
Although this selective signal might be due to rs17630235 proximity with ATXN2, it is tempting to speculate about the changed role of ALDH2 in a post-neolithic society, which made available several fermentable substrates. Other selected SNPs include rs74555583 and rs11539148, both associated with sleep patterns (chronotype). Most notably, the latter is a missense variant in the catalytic domain of QARS1, for which also functions as splicing QTL.
QARS1 itself encodes an enzyme involved in the glutaminyl-tRNA synthesis and, when mutated, leads to microcephaly, cerebral-cerebellar atrophy, and seizures.
Here we combined existing knowledge on genotype-phenotype associations and the availability of ancient genomes to assess the impact of ancient migrations on the phenotypic landscape of contemporary Europeans. We leveraged on traits measured in living individuals, complementing previous works where phenotypes were inferred for ancient genomes instead. As a whole, the most affected traits include pigmentation and anthropometric traits together with blood cholesterol levels, caffeine consumption, heart rate, and age at menarche.
Importantly, while our genome-wide results highlight an overall excess of an ancestry in the carriers of a given phenotype, this is not necessarily mirrored at the genetic loci for which the genotype-phenotype association is ascertained in the literature. A genome-wide excess can completely explain a regional signal, leading to non-significant Z-scores, and even indicate a different direction for the same ancestry. While the first scenario can be due to the extreme polygenicity of a trait, possibly coupled with an inaccurate tagging of the actual functional regions by the GWAS catalog hits, the second might indicate an incomplete correction of non-genetic factors in the genome-wide analysis. Indeed, it is possible that place of residence and educational attainment alone cannot fully account for confounding environmental effects such as socioeconomic status. Conversely, candidate region Z-scores are disentangled from background confounders and virtually free from collinearity issues when they also agree with the relevant ICs. In this light, we here chose to report and discuss results showing region-specific significance for covAs and matching ICs (as reported in Table S4), thus refraining from making inferences on traits such as eye pigmentation in Yamnaya, among others.
WHG ancestry in present day individuals is linked to lower cholesterol levels, higher BMI, and putatively contributed brown hair and light eye color to the contemporary Estonian population. This last association has been previously described based on the HERC/OCA2 haplotypes found in ancient WHG samples.
In addition, loci associated with these features also appear to have undergone selection in Estonians. Other region-specific associations for this ancestry include decreased hip circumference and increased caffeine consumption and heart rate.
An enriched Yamnaya ancestry is linked to a strong build, with tall stature (in agreement with previous studies
) and increased hip and waist circumferences, both at genome-wide and region-specific levels, but also to black hairs and high-cholesterol concentrations when focusing on candidate regions. The associations of Yamnaya and WHG ancestries to respectively higher and lower cholesterol levels, together with the observed signatures of selection at loci connected to cholesterol and BMI, add a new component to our understanding of post-neolithic dietary adaptation
with potential implications to disease risk and outcomes in present-day populations.
Anatolia_N enrichment in trait-related genomic regions is connected with a reduced BMI-corrected waist-to-hip ratio, reduced BMI, light (but not green) eyes and fair hair, increased age at menarche, and reduced heart rate. Notably, covA(i,Anatolia_N) has a substantial weight only in IC2, the single IC that reaches significance when predicting heart rate, suggesting a prominent role for this ancestry in determining this trait.
Finally, the Siberian ancestry is connected with dark hair pigmentation, higher heart rate, lower caffeine consumption, and most prominently, green eye color and lower age at menarche. Importantly, while the results connected to the Siberian ancestry are not of broad applicability to all European populations, covA(i,Siberia) and relative ICs received effect-sizes with mixed significance in all the previous traits except for age at menarche and pigmentation, suggesting that other ancestries might have a larger impact. In other words, we do not find other phenotypes that can be best explained by similarity with Siberia, implying that the presence of this ancestry in the Estonian genome does not significantly affect the inference based on the other, pan-European ancient components.
A general caveat about significance levels observed in this study is that as we refrain from reducing interdependent traits by arbitrary choices, even testing multiple alternatives of the same trait, we expose ourselves to inflated false negatives. We deemed it best to acknowledge and control this risk by avoiding overly stringent multiple testing corrections as Bonferroni and adopting the Benjamini-Hochberg procedure to control FDR instead. In addition, as highly significant traits tend to have higher heritability, it is likely that our analysis might not have enough statistical power for poorly heritable traits. Nevertheless, as we are able to highlight ancestry-trait associations for caffeine consumption
, brown hair color
, and even heart rate
, this condition should be limited only to the very few traits exhibiting lower heritabilities.
Importantly, our inferences are applicable to contemporary individuals of European ancestry, where the phenotypes were collected. Conversely, using them to extrapolate features of ancient populations, although tempting, should be done with caution due to the interaction of their genetic legacy with a radically different lifestyle and environment. Furthermore, when seeking a biological interpretation of our results, it should be kept in mind that certain narrowly defined, contemporary phenotypes such as caffeine consumption may point to broader biological pathways.
Taken together, our results show that the ancient components that form the contemporary European landscape were differentiated enough at a functional level to contribute ancestry-specific signatures on the phenotypic variability displayed by contemporary individuals, regardless of which target population one may examine. In particular, when looking at Estonians, for 11 out of 27 traits surveyed here we could confirm a significant relationship between presence of a given ancestry in genetic regions associated with a given phenotype and how this is expressed by contemporary individuals. While showing that both autochthonous (WHG) and incoming groups contributed genetic material that shapes the phenotype landscape observed today, we also demonstrated that a subset of these loci further underwent positive selection in the last 500 generations. Although not determining whether the selected alleles (and phenotypes) were predominantly contributed by the autochthonous or incoming groups, by connecting genotypic ancestry and complex traits measured in a large dataset, our results reveal both neutral and adaptive consequences of the post-neolithic admixture events on the European phenotype landscape.
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Davide Marnetto ([email protected]).
Materials availability
This study did not generate new unique reagents.
Experimental model and subject details
We relied on a set of 50,341 individuals from the Estonian Biobank
with coupled genotypic and phenotypic data (contemporary Estonian sample set). Among them, 49,799 samples were extracted for the following analyses as described in ”Method details” section, including 32,709 adult (age
18) females and 17,090 adult males, with median age 43. The Estonian Biobank managed recruitment, inclusion and exclusion criteria, and obtained informed consent for all participants in the study. Data were accessed upon approval of the Estonian Committee on Bioethics and Human Research: Approval Number 285/T-13 obtained on 17/09/2018 by the University of Tartu Ethics Committee.
Method details
Sample selection and ancient European grouping
After removing second-degree relatives (pi-hat
) from the contemporary Estonian sample set (n = 50,341) we obtained a subset of 37,952 individuals and used it as a scaffold to perform a PC Analysis (PCA) with Eigensoft-6.1.4. Other individuals were projected on the same PCA space. Outliers identified in this process (with parameters numoutlieriter: 5 numoutlierevec: 10 outliersigmathreshold: 6) were discarded. Samples that on the first round of genome-wide covAs were more distant than 8 Interquartile Ranges (IQR) from the upper or lower quartile against any of the ancestries were also discarded, resulting in 49,799 individuals included in our sample set. For each trait of interest we first removed individuals with missing data for traits and covariates and subsequently discarded second-degree relatives.
To define ancestral European groups we started from the Allen Ancient DNA Resource (AADR) V44.3 merged with present-day individuals typed on the Human Origins array (see Data Availability section). From this set we defined a manually curated core set for each ancestral group, then performed a PCA on a space defined by modern Eurasian and North African individuals west of Iran (included), where the ancient samples were projected. We expanded these core sets to other individuals from AADR dataset using multi-dimensional ellypses with diameters equal to 3 core set SDs. We used 4 dimensions: the annotated dating and the first 3 PCs generated above. With this process we selected 90 WHG, 92 Anatolia_N, 74 Yamnaya S1. In addition, from the ones available from the same dataset, we took 7 samples as representative of the broader Siberian ancestry, assuming any Siberian individual would be equidistant to the other ancestral European groups: S_Even-3.DG, S_Even-1.DG, S_Even-2.DG, Bur1.SG, Bur2.SG, Kor1.SG and Kor2.SG. 957,869 SNPs remained in our dataset after merging the contemporary and ancient sets.
covA definition
covA is the covariance in allele frequency (p) within a contemporary individual i (i.e., its allele dosage) with the ancestral group of interest j, computed respectively against the allele frequency
of the contemporary population C and the average frequency
in all the A ancient groups:
where j is one of the four ancestral groups, the statistics are different but strongly correlated (data not shown): this is expected when the
outgroup population is an outlier to all populations, contemporary and ancestral, considered in covA, as in
. Indeed
has a strong relationship with f-statistics,
i.e.,
where C is the contemporary population (Estonians in our case) and A is an ideal population with
. Nevertheless, as opposed to f-statistics, which include allele frequencies in groups that portray actual populations,
includes
, an average allele frequency which only serves as balanced comparison for the ancestries under analysis. In relation with our aim, this constitutes an advantage of covA, which does not take into account drift or selection occurred in the branch that connects the outgroup population with the internal node shared by the other populations under analysis. Related to this, see our simulation in Figures S1A and S1B. See also covA distributions in Estonia in Figures S1C–S1F.
(date of download: 20/11/2020) and then extracted for each trait a set of hits connected to it filtering on the reported trait (”TRAIT/DISEASE” field) or selecting the appropriate trait in the Experimental Factor Ontology (EFO) field, as specified in Table S3. Then we took windows of 5, 50 and 500 Kbs centered on the selected hits and merged them where overlapping, obtaining three sets of candidate regions for each trait.
to simulate a contemporary population of 5000 individuals with the following composition: 50% EUR1, 30% EUR2, 20% EUR3. These are simulated groups which separated 10kya and admixed 3kya. Each genome consists of 1000 independent loci of 50kb. We then simulated complex traits as following: we extract randomly one causal SNP in the central quartiles of each window, among those with MAF
, and assign standard-distributed effect sizes. We forced an ancestry-trait interaction by enforcing uniformly distributed correlations between the effect sizes and the causal allele frequency specificity, i.e., the allele frequency difference between the contributing ancestry and the other ancestries. As we modeled a completely symmetrical admixture among the ancestral groups, for each simulated trait we could directly measure its shift in the ancestry j as a proxy for the true ancestry-trait association. This directional shift is described as average trait value for that ancestral group standardized over all ancestral individuals
. See also Figures S1A and S1B.
Phenotypes treatment and heritability
Continuous traits were treated as specified in Table S2 and regressed against the covariates according to the same table. Individuals with traits or covariates more distant than 4 IQRs from the upper or lower quartile were considered as outliers and discarded. After adjusting traits as described, their heritability was computed using LDAK 5.0.
First we computed a kinship matrix with the LDAK-Thin Model: we thinned down SNPs on the non-related sample set defined above with parameters–window-prune 0.98–window-kb 100, then used–calc-kins-direct with the resulting weights and–power 0.25. Finally we estimated heritability using REML solver.
Testing for ancestry-trait association with covA and covA-based PCs
We fitted each standardized trait
with a model including one standardized covA for each ancestry j and estimated its coefficient:
. We adopted a logistic regression for categorical traits, which were transformed to {0, 1} where 1 stands for the specified category and 0 for all the others. covAs were computed on the corresponding candidate region for each trait. In addition, each trait was regressed against three PC-transformed covAs:
. Notably, we transformed all covAs using the loadings obtained from a PC analysis run on whole genome covAs, thus obtaining components that were largely independent, yet not strictly principal. These Independent Components (ICs) were standardized and included together as predictors. To evaluate association we used coefficient Z-scores computed against the same parameter extracted from 50 random genomic sets with matching size (double-sided Z-test). Sample sizes are reported in Figure 1 and Table S2.
In the genome-wide analysis, we adopted similar steps, performing individual regressions for all the covAs and coupling this with a model including all covAs PCs, but socioeconomic variables were added as covariates in all models as described in the result section. Note that in this analysis ICs are not needed anymore, but actual PCs are used. Then, the standardized coefficient (β or effect size), or the Odds Ratio (OR) were directly used to assess ancestry-trait association for continuous and categorical traits respectively (double-sided coefficient p value from glm R function). This analysis was restricted to samples for which socioeconomic covariates were defined, i.e., 38,996 samples (including relatives): the actual sample size for this analysis is therefore less than reported in Figure 1 and Table S2.
Significance was evaluated at Benjamini-Hochberg False Discovery Rate = 0.05 on the results from IC/PC regressions. See comments on the implications of this definition in the ”Discussion” section.
Testing for signals of positive selection
In order to test individual SNPs for signatures of positive selection we utilized the Relate/CLUES pipeline.
This was applied on a curated subset of 1800 unrelated samples; further details on its application are described in ”Relate/CLUES application” section. CLUES was run once for each of the 14,712 unique GWAS hits for traits analyzed here with a derived allele frequency (DAF) above 1% and passing the 1000 Genomes strict mask (SNP size for each trait given as number in Figures 4 and S4). To obtain an expected distribution we randomly sampled 10,000 GWAS hits from the GWAS catalog meeting the same conditions and ran CLUES for positions not present among the 14,712 SNPs. Next, for each phenotype we compared its distribution of the logLR values to that of random GWAS hits. We took 1000 random subsets (with replacement) from the 10,000 logLR values each of the same length as the number of GWAS hits for a given phenotype and ranked the logLR values from lowest to highest within each subset. In this way we obtained 1000 values for each logLR rank from 1 to N where N is the number of SNPs analyzed for a given phenotype. For each rank we calculated the mean and the
and
percentiles. Finally, we rank SNPs within each trait and compare each logLR value to the mean and
percentiles range for the corresponding rank of the background distribution. As we are interested in deviations in the higher ranks we focus on the top 100 ranks for each phenotype. Such an approach is conservative as we are testing not against presumably neutral SNPs but against random GWAS hits that are shown to be enriched in signals on natural selection compared to random SNPs in the genome.
mutation rate of 1.25x10-8 and effective population size of 30000. We then removed: a) 115 related samples; b) PCA outliers (2.5% in both directions along 1st and 2nd PCs both for Estonian-only PCA and when projecting Estonian samples on a “European” PC space); c) singleton count outliers (top/bottom 2.5%); d) samples with pairwise total IBD sharing
cM with more than one other sample; d) 57 random samples. These filtering steps lead to a set of 1800 samples and are described in Pankratov et al., 2020.
We estimated coalescence rate through time (EstimatePopulationSize.sh module) with a random subset of 100 samples (for which subtrees were extracted with RelateExtract SubTreesForSubpopulation), mutation rate 1.25x10-8, generation time 28 years, number of iterations 5, tree dropping threshold 0.5 and time bins defined as 10x years ago where x changes from 2 to 7 with an increment of 0.1. We ran CLUES
using Relate trees for the 1800 individuals as input. We first sampled branch length (SampleBranchLengths.sh module) for the tree containing a SNP of interest (200 samples). We removed 100 of 200 trees as burn-in and selected every 5th tree for importance sampling by CLUES. We focused on the time window between 0 and 500 generations ago.
Data and code availability
This paper analyzes existing, publicly available data. The accession numbers for these datasets are listed in the key resources table. Data from Estonian Biobank are under managed access and subject to approval of the Estonian Committee on Bioethics and Human Research. All original code has been deposited at https://bitbucket.org/dmarnetto/anccontributions and is publicly available as of the date of publication.
Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.
Acknowledgments
This work is supported by the European Union through the European Regional Development Fund, project no. 2014-2020.4.01.16-0024, MOBTT53 (D.M., K.P., L.M., and L.P.); MOBEC008 (V.P., M.Mondal, M.Metspalu, and A.E.); 2014-2020.4.01.16-0030 (F.M. and M.Metspalu); 2014-2020.4.01.15-0012 (M.Metspalu); through the Horizon 2020 research and innovation program grant no. 810645 (V.P., M.Mondal, M.Metspalu, and A.E.); and through the Horizon 2020 MSCA Initial Training Network, grant no. 765937 (R.C.). L.S. and M.Metspalu are supported by the Estonian Research Council through PUT PRG243. S.M. is supported by the STARS@UNIPD 2019 Consolidator Grant for the project CircadianCare. We would like to thank Paolo Provero and Fernando Racimo for insightful discussion and comments. Most of the analysis was run on the High-Performance Computing Center at the University of Tartu. The Estonian Biobank Research Team includes Mari Nelis, Lili Milani, Tōnu Esko, Andres Metspalu, and Reedik Mägi.
Author contributions
D.M. and L.P. conceived and designed the study; A.E. contributed in the statistical design; D.M. and V.P. performed data analysis; M.Mondal, F.M., K.P., L.V., L.M., and L.P. contributed to data analysis; S.M. and R.C. provided analyses and expertise about sleep traits; F.M., L.S., L.L., and M.Metspalu contributed with ancient genetics expertise; D.M. and L.P. drafted the manuscript; all authors reviewed and approved the submitted paper.
Table S1: Ancestral reference samples, related to STAR Methods
The table contains for each ancient sample used in the study: ID, ancestral group assignation, and if the sample is part of the core set for that group.
Table S2: Traits analyzed, related to STAR Methods
The table contains the traits analyzed and how they were encoded, corrected, and adjusted for covariates. In addition, for each trait we include sample size and case number for categorical traits.
References
Lazaridis I.
Patterson N.
Mittnik A.
Renaud G.
Mallick S.
Kirsanow K.
Sudmant P.H.
Schraiber J.G.
Castellano S.
Lipson M.
et al.
Ancient human genomes suggest three ancestral populations for present-day Europeans.
00108-7&id=gr1.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00108-7&id=gr2.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00108-7&id=gr3.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00108-7&id=gr4.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}