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.Previous studies inferred phenotypes in these source populations, using either a few single locior polygenic scores based on genome-wide association studies,and investigated the strength and timing of natural selection on lactase persistence or height, among others.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 biobankand 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.

The recent and putatively ongoing nature of the inferred selective pressure on the six traits shown in Figure 4 A 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 4 B. These include some loci previously shown to be selected in West Eurasians (rs4988235 at MCM6/LCT,pigmentation-related SNPs at HERC2/OCA2, TYRP1, TYR, TPCN2,and rs653178 at ATXN2) 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 tissuesfor ALDH2, an aldehyde dehydrogenase known for its role in the alcohol metabolism.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.

We independently asked whether the phenotype-associated regions above also exhibit signs of recent natural selection. We applied CLUESto 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 4 A; 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.

(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.

(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

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.In other words, the local admixture imbalances we detected so far are not necessarily transferred to the whole population.

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 ).

(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 reporteffect 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.

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 reasonsare still visible in covA (see Figures S1 G–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.

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 2 B). We used the loadings from a principal component (PC) analysis on whole genome covAs ( Figure 2 C) 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 2 A) highlight the overall excess or lack of certain ancestries in relation with a given phenotype but are largely intertwined, ICs ( Figure 2 B) 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.

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.

(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.

(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 = 2 .

We defined three sets of candidate regions for any given trait by considering windows of 5kb, 50kb, or 500kb centered around GWAS cataloghits 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.

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.Nevertheless, our dataset captures a genetic basis for most of them, as confirmed by the trait heritability measured in our sample ( Figure 1 ).

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.

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 frequencyand because of varying grades of similarity among the ancestries for historical demographic reasons (see covA joint distributions in Figures S1 C–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 S1 G). Conversely, as shown by median covA for the Siberian component in Figure S1 I, 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 S1 H and S1J).

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 S1 A). See STAR Methods for further discussion ofproperties and simulation details, including the definition of ancestral specificity.

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.

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 (andrespectively):

We identified 27 complex traits of interest, based on information availability in the Estonian Biobank(EstBB) and genome-wide association studies (GWAS) catalog.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.

Discussion

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.

5 Olalde I.

Allentoft M.E.

Sánchez-Quinto F.

Santpere G.

Chiang C.W.K.

DeGiorgio M.

Prado-Martinez J.

Rodríguez J.A.

Rasmussen S.

Quilez J.

et al. Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European. , 23 Key F.M.

Fu Q.

Romagné F.

Lachmann M.

Andrés A.M. Human adaptation and population differentiation in the light of ancient genomes. 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.

6 Mathieson I.

Lazaridis I.

Rohland N.

Mallick S.

Patterson N.

Roodenberg S.A.

Harney E.

Stewardson K.

Fernandes D.

Novak M.

et al. Genome-wide patterns of selection in 230 ancient Eurasians. , 8 Cox S.L.

Ruff C.B.

Maier R.M.

Mathieson I. Genetic contributions to variation in human stature in prehistoric Europe. 7 Saag L.

Vasilyev S.V.

Varul L.

Kosorukova N.V.

Gerasimov D.V.

Oshibkina S.V.

Griffith S.J.

Solnik A.

Saag L.

D’Atanasio E.

et al. Genetic ancestry changes in Stone to Bronze Age transition in the East European plain. , 30 Buckley M.T.

Racimo F.

Allentoft M.E.

Jensen M.K.

Jonsson A.

Huang H.

Hormozdiari F.

Sikora M.

Marnetto D.

Eskin E.

et al. Selection in Europeans on Fatty Acid Desaturases Associated with Dietary Changes. , 31 Mathieson S.

Mathieson I. FADS1 and the Timing of Human Adaptation to Agriculture. 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 adaptationwith 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 ( h 2 = 0.087 ± 0.009 ) , brown hair color ( h 2 = 0.079 ± 0.009 ) , and even heart rate ( h 2 = 0.044 ± 0.009 ) , 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.