Skin Pigmentation-Associated SNPs Show a Signal of Positive Selection in Ancient West Eurasians.

We obtained skin pigmentation-associated SNPs from the UK Biobank GWAS for skin color released by the Neale Lab (14). We analyzed the evolution of these variants using two datasets of publicly available ancient DNA (Fig. 1 A–C and SI Appendix, Fig. S1). The first (“capture-shotgun”) consists of data from 1,158 individuals dating from 45,000–715 y BP genotyped at ∼1.2 million variants (7, 16, 18192021222324252627282930313233343536373839404142434445464748495051–52). The second (“shotgun”) is a subset of 248 individuals with genome-wide shotgun sequence data (16, 18192021222324–25, 28, 32, 35, 38, 40414243–44, 47, 49, 51, 53545556575859–60). This smaller dataset allows us to capture variants that may not be well tagged by the genotyped SNPs in the capture-shotgun dataset. Weighting variants by their GWAS-estimated effect sizes, we show that the polygenic score—the weighted proportion of dark pigmentation alleles—decreased significantly over the past 40,000 y (Fig. 1 A and B, P < 1 × 10⁻⁴ based on a genomic null distribution and accounting for changes in ancestry; SI Appendix, Fig. S2). We recapitulate this significant decrease using a different set of UK Biobank GWAS summary statistics calculated using linear mixed models (61) (SI Appendix, Fig. S3) and identifying independent signals based on predefined approximately independent linkage disequilibrium (LD) blocks (62) (SI Appendix, Fig. S4), rather than LD clumping. The genetic scores of the oldest individuals in the dataset fall within the range of present-day West African populations, showing that Early Upper Paleolithic (∼50–20,000 y BP) individuals, such as Ust’Ishim, carried few of the light skin pigmentation alleles that are common in present-day Europe.

Fig. 1.

Genetic scores for skin pigmentation over time. The solid lines indicate fitted mean scores with gray 95% confidence intervals. Weighted scores based on UK Biobank SNPs and unweighted scores based on manually curated SNPs for samples in the (A and D) shotgun dataset, (B and E) capture-shotgun dataset, and (C and F) capture-shotgun dataset dated within the past 15,000 y. Scores of labeled samples are plotted but not included in regressions. Averages for 1000 Genomes superpopulations are plotted at time 0. Area of points scale with number of SNPs genotyped in the individual.

Over the past 15,000 y, which covers most of the data, the polygenic score decreased (P < 1 × 10⁻⁴) at a similar rate to the estimated rate over the 40,000-y period (5.27 × 10⁻⁵ per year vs. 3.47 × 10⁻⁵ per year; P = 0.04). Stratifying individuals by ancestry (hunter-gatherer, Early Farmer, and Steppe) using ADMIXTURE revealed baseline differences in genetic score across groups (Fig. 1C). Mesolithic hunter-gatherers carried fewer light pigmentation alleles than Early Farmer or Steppe ancestry populations. This difference (0.091 score units) was similar to the difference between Mesolithic hunter-gatherers and Early Upper Paleolithic populations (0.097 score units). We tested for evidence of change in score over time within groups while controlling for changes in ancestry, finding no significant change within hunter-gatherer or Early Farmer populations (P = 0.32; P = 0.29). Steppe ancestry populations, however, exhibited a significant decrease (P = 0.007), with a steeper slope than that of Early Farmer or hunter-gatherer (P = 0.004 and P = 0.08, respectively).

While the UK Biobank is well powered to detect variation that is segregating in present-day Europeans, it may not detect variants that are no longer segregating. We therefore manually curated a separate list of 18 SNPs identified as associated with skin pigmentation from seven studies of populations representing diverse ancestries (Dataset S1E). Using this set of SNPs to generate unweighted scores, we observed similar results as with the UK Biobank ascertained SNPs (Fig. 1 D–F). We find a significant decrease in genetic score over both 40,000 and 15,000 y in the capture-shotgun dataset (P < 1 × 10⁻⁴ for both), although the decrease in the shotgun dataset was not significant (P = 0.17), probably reflecting a lack of power from the small number of ancient samples and SNPs.

Ancestry and Selection Have Different Effects on the Histories of Each Variant.

Next, we investigated the evolution of each variant independently. We separated the effects of ancestry and selection in the data by regressing presence or absence of the alternate allele for each pigmentation-associated SNP on both date and ancestry, inferred using principal component analysis (PCA) (63). Significant effects of ancestry on frequency can be interpreted as changes in allele frequency that can be explained by changes in ancestry. Significant changes in frequency over time, after accounting for ancestry, can be interpreted as evidence for selection occurring during the analyzed time period.

The SNP rs16891982 at the SLC45A2 locus shows the strongest evidence for selection across analyses (Fig. 2), but we note this does not necessarily imply positive selection over the entire time transect. We observed little evidence for demographic transitions in driving the increase in light allele frequency over time for this SNP. In addition, we found robust signals of selection for the light allele of rs1126809 near TYR, rs12913832 and rs1635168 near HERC2, and rs7109255 near GRM5 (Fig. 2 B and C). In some cases—for example at rs6120849 near EDEM2 and rs7870409 near RALGPS1—light alleles decreased in frequency more than can be explained by changes in ancestry (Fig. 2 B and C). The dark allele at rs6120849 (EDEM2) is also significantly associated with several other traits including increased sitting and standing height in the UK Biobank, increased blood creatinine, and increased heel bone mineral density (64). The dark allele of rs7870409 (RALGPS1) is associated with increased arm and trunk fat percentage (65). These observations raise the possibility of pleiotropic effects constraining the evolution of these loci or of spurious pleiotropy due to uncorrected population stratification in the GWAS.

Fig. 2.

Plot of P values for ancestry and time for individual skin pigmentation SNPs. The dashed lines indicate fifth percentile of P values. UK Biobank SNPs using the shotgun dataset are plotted in A and capture-shotgun dataset in B and C. Manually curated SNPs using the shotgun dataset are plotted in D and capture-shotgun dataset in E and F. SNPs are colored according to whether the light allele is increasing over time (red) or decreasing (blue), with saturation determined by the magnitude of change.

On the other hand, the changes in frequencies of several variants appear to have been driven largely by changes in ancestry. For example, rs4778123 near OCA2, rs2153271 near BNC2, and rs3758833 near CTSC show evidence that their changes in frequency were consistent with changes in genome-wide ancestry (Fig. 2 B and C). At SLC24A5, rs2675345 shows evidence of change with both ancestry and time, suggesting that even after the spread of the light allele from Anatolia into Western Europe in the Neolithic (7) selection continued to occur post-admixture. Again, not all of these cases involve an increase in the frequency of the light pigmentation allele over time. The light allele of rs4778123 (OCA2) was at high frequency in hunter-gatherers but lower in later populations (SI Appendix, Fig. S6B). From the manually curated set of SNPs (Fig. 2 D–F), rs12203592 near IRF4 also displays a marked effect of ancestry with higher light allele frequency in hunter-gatherers (SI Appendix, Fig. S6E). While rs12203592 was not present in the UK Biobank summary statistics, another SNP at the IRF4 locus, rs3778607, was present but with a smaller ancestry effect (SI Appendix, Fig. S7).

Lack of PBS Selection Signal across Pigmentation SNPs in Source Populations of Europeans.

As we described in the previous section, the frequencies of several skin pigmentation SNPs appear to have been driven by changes in ancestry. This observation suggests that their frequencies have diverged between the European source populations and then subsequently been driven by admixture. Frequency differences across the source populations might reflect the action of either genetic drift or selection. To examine this question, we looked for evidence of selection at skin pigmentation SNPs in each ancestral population using the population branch statistic (PBS) test (66). We ran ADMIXTURE with K = 3 on the capture-shotgun dataset (SI Appendix, Fig. S8) and used the inferred source population allele frequencies to calculate PBS for each source in 20-SNP nonoverlapping windows, with the 1000 Genomes CHB (Han Chinese) and YRI (Yoruba) populations as outgroups.

Few of the skin pigmentation loci show extreme PBS values (Fig. 3), and even fewer show evidence of divergent frequencies across ancient groups (SI Appendix, Fig. S10). On the Early Farmer and Steppe ancestry branches, but not the hunter-gatherer branch, the SLC24A5 locus exhibited the strongest signal (Fig. 3 B and C), indicating selection at the locus in both Early Farmer and Steppe source populations. The RABGAP1 locus also exhibited an elevated signal of selection in Early Farmer and Steppe, which remained when using a 40-SNP window size (SI Appendix, Fig. S11 B and C). In all three ancestry groups, we observed an elevated PBS at the OCA2 locus, with the most extreme value on the hunter-gatherer branch. However, unlike the SLC24A5 variant, neither of these SNPs (rs644490 and rs9920172) showed a substantial effect of ancestry from the individual SNP regression analysis. As a whole, the PBS values of the 170 skin pigmentation SNPs do not significantly deviate from the genome-wide distribution of PBS for any of the ancient source populations (SI Appendix, Fig. S12).

Fig. 3.

Joint PBS distributions for 20-SNP windows across the genome for GBR-CHB-YRI on the x axis and X-CHB-YRI on the y axis, with X being (A) hunter-gatherer, (B) Early Farmer, and (C) Steppe. The boxes represent windows centered around UK Biobank skin pigmentation SNPs and are colored redder according to how much higher the light allele frequency is in X compared to YRI. Nearest genes are labeled. The orange and red lines represent the top 1 and 0.1 percentiles.

Signals of Selection Are Restricted to the SNPs with the Largest Effect Sizes.

Since our results suggest that pigmentation-associated variants exhibit different changes in allele frequency, we further examined the signal of selection found in Fig. 1. In general, SNPs with larger GWAS effect sizes changed more in frequency over time than those with small effect sizes (P = 2.8 × 10⁻¹⁸; adjusted R² = 0.36) (SI Appendix, Fig. S13A). To understand the contribution of these large GWAS effect size SNPs, we iteratively removed SNPs from the polygenic score calculations (Fig. 4A). Removing the SNP rs16891982 at SLC45A2, we still found a significant decrease of score over time (P < 1 × 10⁻⁴), but the estimated rate of change (β_time) decreased by 39%. Removing the top two SNPs at the SLC45A2 and SLC24A5 locus further attenuated the signal with a decrease in β_time of 58% (P = 3.96 × 10⁻⁴). We removed the five SNPs with the largest GWAS effect sizes and found the effect of time attenuated (P = 0.026) with a decrease in β_time of 78%. Removing the top 10 SNPs attenuated the signal even more (P = 0.050) with an 81% decrease and removal of 15 or 20 SNPs abolished the signal (P = 0.147, P = 0.52).

Fig. 4.

(A) Regressions of genetic score based on UK Biobank SNPs using the capture-shotgun dataset over date, with scores using all SNPs and iteratively removing top GWAS effect size SNPs. (B) Q_x empirical −log₁₀(P values) using UK Biobank skin pigmentation-associated variants with all 1000 Genomes populations. Different numbers of SNPs were used to calculate Q_x, which were ordered by GWAS-estimated effect size. (C) Variance explained by UK Biobank SNPs with nearest gene labeled between UK and YRI populations (y axis) and within the UK population (x axis).

Using present-day populations from the 1000 Genomes Project (67), we tested for polygenic selection with the Q_x statistic (10), which measures overdispersion of trait-associated SNPs relative to the genome-wide expectation. We performed the test with different numbers of SNPs ordered by largest to smallest absolute GWAS effect size (Fig. 4B). Using all 170 SNPs in the test, we failed to observe a signal of polygenic selection (P = 0.27). However, when we restricted to the top five largest GWAS effect size SNPs, we found a significant signal (P < 2.5 × 10⁻⁷). We observed significant Q_x statistics at P < 0.05 up to the top 30 SNPs, but it appears this signal is driven by the top SNPs of largest effect. After removing the top 5 SNPs and testing the remaining 25 top SNPs, the signal of polygenic selection disappeared (P = 0.59). Examining the predicted variance explained per SNP (Fig. 4C and SI Appendix, Fig. S14), we predict that, of the UK Biobank SNPs, two SNPs at SLC45A2 and SLC24A5 make the largest contribution to the difference between West African and European populations. Most UK Biobank skin pigmentation-associated SNPs are predicted to contribute relatively little to the between-population variance.

Capture-Shotgun Ancient DNA Dataset.

We downloaded version 37.2 of the publicly available datasets of 2,107 ancient and 5,637 present-day (Human Origins dataset) samples from the Reich Lab website: https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers. Ancient individuals were treated as pseudohaploid (i.e., carry only one allele) because most samples are low coverage. For many of these samples, 1,233,013 sites were genotyped using an in-solution capture method (88), while those with whole-genome shotgun sequence data were genotyped at the same set of sites. For our regression models, we included ancient individuals from West Eurasia, defined here as the region west of the Ural Mountains (longitude, <60° E) and north of 35° N (SI Appendix, Fig. S1). We excluded samples with less than 0.1× coverage. We removed duplicate and closely related samples (e.g., first-degree relatives), by calculating pairwise identity by state (IBS) between ancient individuals. For each individual, we identified a corresponding individual with the highest IBS. Based on the distribution of highest IBS values, we identified the pairs with IBS Z score >7 as duplicates and retained the individual with fewest missing sites. We also manually removed related samples that were explicitly annotated, selecting the highest coverage representative. Finally, we incorporated additional ancient samples from the Iberian Peninsula from two recent papers (36, 45), applying the same criteria for coverage and relatedness. Ultimately, we compiled a list of 1,158 ancient pseudohaploid individuals in our analyses covering a time span of the last 40,000 y. Only 11 individuals were dated to >15,000 y BP, and 1,147 individuals lived during the last 15,000 y.

Shotgun Ancient DNA Dataset.

We organized a dataset of samples that were shotgun sequenced without enrichment using the capture array (23, 25262728293031–32, 35, 39, 42, 45, 47484950–51, 54, 56, 58, 60616263646566–67). To generate the first 10 PCs from sites on the 1240K capture array, we combined shotgun samples with available genotype data from v37.2 of the Reich laboratory dataset and 15 individuals that we ourselves pulled down from BAM files. We manually removed duplicate samples in the shotgun dataset, preserving the sample with higher coverage. We checked for first-degree relatives in the dataset by considering familial annotations but found none. The final shotgun dataset included 249 individuals.

Skin Pigmentation SNP Curation.

We obtained summary statistics for UK Biobank GWAS for skin color (data field 1717) from the publicly available release by the Neale Lab (version 3, Manifest Release 20180731) (14). The GWAS measured self-reported skin color as a categorical variable (very fair, fair, light olive, dark olive, brown, black). To identify genome-wide significant and independent SNPs, we performed clumping using PLINK v1.90b6.6 (89) with 1000 Genomes GBR as an LD reference panel (--clump-p1 5 × 10⁻⁸ --clump-r2 0.05 --clump-kb 250), and followed up with clumping based on physical distance to exclude SNPs within 100 kb of each other. We made two separate lists of UK Biobank SNPs for the shotgun and capture-shotgun datasets because the capture-shotgun dataset was restricted to the 1240K array sites. For the capture-shotgun dataset, we intersected all UK Biobank SNPs with 1240K array SNPs before identifying 170 independent and genome-wide significant SNPs (Dataset S1A), which were used in Figs. 1 B and C, 2 B and C, 3, and 4. For the shotgun dataset, we identified 242 independent and genome-wide significant SNPs (Dataset S1B), which were used in Figs. 1A and 2A. For SI Appendix, Fig. S4, we identified 93 SNPs from the Neale Lab GWAS (Dataset S1F) by picking the top SNPs that are genome-wide significant from nearly independent LD blocks in the human genome specified beforehand (62). For SI Appendix, Fig. S3, we identified 176 SNPs (Dataset S1G) using the same clumping parameters with UK Biobank GWAS summary statistics calculated with fastGWA (61).

We also manually curated a list of skin pigmentation-associated SNPs from the literature (Dataset S1C). We identified 12 suitable studies (12, 13, 15, 80, 90919293949596–97), 9 of which were GWAS conducted in diverse populations (Dataset S1D). We did not include SNPs from all of the considered papers for our analysis because some studies reported associations that were not genome-wide significant (P < 5 × 10⁻⁸). However, we included subsignificant SNPs from a GWAS in Europeans (95) with evidence of replication (P < 5 × 10⁻⁸) in the UK Biobank GWAS. For a GWAS in African populations (12), we selected a SNP for each LD block, which was defined by the study authors, with the greatest support based on multiple lines of evidence (e.g., P value and functional annotations). For a GWAS in South Asians (94), we picked SNPs that were 200 kb away from each other and based the selection of SNPs that were in physical proximity by considering P value. The other studies presented a set of independent SNPs that we considered for inclusion. Lastly, we manually clumped the SNPs curated across studies, selecting SNPs 200 kb apart with more statistical and/or functional evidence than nearby SNPs. Specific details about the choices made on which SNPs to exclude can be found in the “Notes” column for highlighted SNPs in Dataset S1C. We ended with a final list of 18 SNPs for the manually curated set (Dataset S1E), 10 of which are present on the 1240K capture array. We also made pseudohaploid calls for each of the 18 SNPs in the shotgun dataset, picking a random allele from the reads at each site as for the overall dataset.

ADMIXTURE Analysis on Capture-Shotgun Data.

We performed unsupervised ADMIXTURE (98) on the dataset of ancient individuals with K = 3, which we found to produce the lowest cross-validation error for 2 < K < 9. The three identified clusters could easily be identified as corresponding to hunter-gatherer, Early Farmer, and Steppe (also referred to as Yamnaya) ancestry (SI Appendix, Fig. S8).

Time Series Analysis of Genetic Scores.

We calculated genetic scores for each individual from present-day 1000 Genomes populations (67) and ancient individuals. We computed the scores in two ways. First, we weighted by the GWAS-estimated effect sizes. Because of variable coverage across ancient samples, not all SNPs were present in a given sample. To account for missing information in the creation of weighted scores, we devised a weighted proportion in which we divided the realized score over the maximum possible score given the SNPs present in the sample: $$ Scor{e}_{weighted}={{\displaystyle \sum }}_{i=1}^m{d}_i{\beta }_i/{{\displaystyle \sum }}_{i=1}^m{\beta }_i$$, where m is the number of skin pigmentation SNPs genotyped in the individual, d_i is the presence of the dark allele at the ith SNP, and β_i is the GWAS-estimated effect size (Fig. 1 A–C). For the manually curated list of SNPs where we did not have comparable effect size estimates, we computed an unweighted score, effectively assuming that all variants had the same effect size. This score is the proportion of dark alleles an individual carries out of the SNPs used in the construction of the score: $$ Scor{e}_{unweighted}=\left(1/m\right){\displaystyle {\sum }_{i=1}^m{d}_i}$$, where m is the number of SNPs genotyped in the individual and d_i is the presence of the dark allele at the ith SNP (Fig. 1 D–F).

We used logistic regression to examine the association between weighted and unweighted genetic score and time for all ancient samples. We fitted separate models for ancient samples in the shotgun and capture-shotgun datasets. We included ancestry as a covariate in the model by including the first 10 PCs. PCA was performed using smartpca v16000 (63) to generate PCs from present day populations and we projected ancient individuals onto these axes of variation. For the shotgun dataset and the capture-shotgun dataset that included all samples (40,000 y BP), we used 1000 Genomes samples. For the capture-shotgun dataset that restricted to samples later than 15,000 BP, we used West Eurasians from the Human Origins dataset (16).

We model the score of each individual as the proportion of successes in a binomial sample of size m_i. That is, $$ m_i{S}_i\sim Binomial\left(m_i,p_i\right)$$, where, for individual i, m_i is the number of SNPs genotyped, S_i is the (weighted or unweighted) score, the probability of success p_i is given by the following:

and date_i, lat_i, and lon_i are the dates, latitude, and longitude of each individual. For the weighted score, m_iS_i is noninteger, but the likelihood can be computed in the same way.

For the stratified analysis, we divided the capture-shotgun dataset into mutually exclusive ancestry groups. We categorized individuals for this analysis as hunter-gatherer if they carried over 60% of the ADMIXTURE-estimated hunter-gatherer component. For placement into the Early Farmer group, we required that an individual carry over 60% of that component. Individuals we categorized as Steppe had over 30% of the ADMIXTURE component and were dated to less than 5,000 y BP. Individuals that are more recent than 5,000 y BP and have at least 30% of the Steppe ADMIXTURE component were classified as Steppe even if they had more than 60% of the Early Farmer component. Based on these cutoffs, there were 102 hunter-gatherer, 499 Early Farmer, and 478 Steppe ancestry assigned individuals. We used the same logistic regression model as above for the stratified analysis to control for ancestry.

Because the fitted model parameters are overdispersed (SI Appendix, Fig. S2), we computed P values from a genome-wide empirical null distribution for β_date. We made this null distribution by running the regression model described above on scores from sets of random, frequency-matched (±1%) SNPs across the genome, maintaining the same effect sizes for the weighted score. We matched the derived allele frequency based on the EUR superpopulation frequencies from 1000 Genomes. We reported P values based on 10,000 random samples.

Time Series Analysis of Individual SNP Allele Frequencies.

We performed logistic regression for each individual SNP using date of the sample and ancestry as covariates. That is, the probability that the haplotype sampled from individual carries the derived allele is p_i, where

We compared the full model above to a nested model with no principal components by performing a likelihood ratio test [in R we use anova(nested model, full model, test = ‘Chisq’)] to obtain a P value for the ancestry term which encompasses β_PC1…β_PC10. To obtain a P value for β_date, we compare a nested model without β_date to the full model.

Q_x Polygenic Selection Test.

We used the Q_x test for directional, polygenic selection on skin pigmentation (10). For this test, we used 170 skin pigmentation SNPs obtained from the UK Biobank in all 1000 Genomes populations. We restricted sites to those on the 1240K and constructed a covariance matrix from a total of one million SNPs. To calculate empirical P values, we sampled a total of 500,000 null genetic values matching skin pigmentation SNPs based on a ±2% frequency on the minor allele. However, we report one million runs for the top seven SNPs and two million runs for the top five and six to better distinguish the deviation of the tested SNP sets from the expectation under drift.

Variance in Skin Pigmentation Explained by Individual SNPs.

We estimated the variance explained in the phenotype by a SNP within the UK population using the following:

where β is the UK Biobank GWAS-estimated effect size and f is the allele frequency in the GBR population from 1000 Genomes. We calculated the proportion of variance explained as follows:

where N is the GWAS sample size (356,530) and σ is the SE of the estimated β (99). To estimate the variance explained between GBR and YRI, we calculated the between population variance as follows:

where f₁ and f₂ are the allele frequencies of GBR and YRI, β is the UK Biobank GWAS-estimated effect size, and f_avg = 0.5(f₁ + f₂).

Population Branch Statistic.

We calculated the PBS (66) for different sets of ancient groups and present-day groups. For the ancient populations, we used inferred source population allele frequencies from ADMIXTURE (the P matrix), whereas for present-day groups, we used the observed allele frequencies estimated from 1000 Genomes YRI, CHB, and GBR. Using these frequencies, we estimated pairwise Hudson’s F_ST between groups for the PBS test. The estimator is calculated as follows:

where for population i, n_i is the sample size and p_i is the allele frequency of the sample. We excluded SNPs that had a missing rate of >90% in the capture-shotgun dataset. Because of variable coverage in ancient genomes, not all sites were used in the calculation. F_ST was calculated using 20- and 40-SNP nonoverlapping windows throughout the genome. F_ST values were transformed to calculate genetic divergence between populations as T = −log(1 – F_ST). We calculated PBS for ancient group X (i.e., hunter-gatherers, Early Farmers, or Steppe) using 1000 Genomes Han Chinese (CHB) and Yoruba (YRI) as follows:

We identified nearby genes by examining the National Center for Biotechnology Information RefSeq track for assembly GRCh37/hg19 (downloaded from https://genome.ucsc.edu/cgi-bin/hgTables) for genes that overlap with windows with high PBS values. We note that the gene names are used simply as labels for each locus and do not necessarily represent the causal genes.