Holocene Selection for Variants Associated With General Cognitive Ability: Comparing Ancient and Modern Genomes

Abstract

Human populations living during the Holocene underwent considerable microevolutionary change. It has been theorized that the transition of Holocene populations into agrarianism and urbanization brought about culture-gene co-evolution that favored via directional selection genetic variants associated with higher general cognitive ability (GCA). To examine whether GCA might have risen during the Holocene, we compare a sample of 99 ancient Eurasian genomes (ranging from 4.56 to 1.21 kyr BP) with a sample of 503 modern European genomes (F_st = 0.013), using three different cognitive polygenic scores (130 SNP, 9 SNP and 11 SNP). Significant differences favoring the modern genomes were found for all three polygenic scores (odds ratios = 0.92, p = 001; .81, p = 037; and .81, p = .02 respectively). These polygenic scores also outperformed the majority of scores assembled from random SNPs generated via a Monte Carlo model (between 76.4% and 84.6%). Furthermore, an indication of increasing positive allele count over 3.25 kyr was found using a subsample of 66 ancient genomes (r = 0.22, p_one-tailed = .04). These observations are consistent with the expectation that GCA rose during the Holocene.

The Holocene (11.7 kyr BP to the present) was a time of considerable microevolutionary change among human populations (Hawks et al., 2007). Based on various converging lines of molecular evidence, it has been estimated that the rate of adaptive evolution among these populations may have been 100 times greater than during the preceding Pleistocene (Cochran & Harpending, 2009; Frost, 2011; Hawks et al., 2007). Novel adaptations that are known to have arisen and spread during this period include lactase persistence (McCracken, 1971) and alterations in hemoglobin, permitting enhanced tolerance to diminished oxygen levels among populations living at high altitude (Bigham et al., 2010).

Other significant adaptations that may have been favored by selection during this period are the learning and problem-solving mechanisms that give rise to general cognitive ability (GCA). At the level of individual differences in cognitive performance, GCA is characterized by positive correlations among different measures of cognitive abilities (Jensen, 1998). GCA is highly heritable, with behavior genetic twin studies indicating that between 22% and 88% of its variance, depending on age, may be due to the action of additive genetic effects (Bouchard, 2004). Consistent with this, studies employing genome-wide complex trait analysis genomic-relatedness-based restricted maximum-likelihood (GCTA GREML) have found that a significant amount of the variance in GCA can be directly attributed to large numbers of variants with individually small, but cumulatively large (additive) effects on the phenotype (Davies et al., 2011).

GCA is not just present in human populations, but is present (at lower levels) also in a variety of animal taxa, including chimpanzees and other primates, raccoons, mice, rats, and ravens (Galsworthy et al., 2014). The existence of these species’ differences furthermore indicates that GCA has been under directional selection, particularly in the primate clade. This seems to be especially true in the case of abilities associated with tool use, which in primates are among the most strongly associated with GCA, in addition to being the most heritable. They are also associated with the strongest signals of recent directional selection (Fernandes et al., 2014; Woodley of Menie et al., 2015). This is consistent with the theoretical expectation (Geary, 2005) that GCA is an adaptation to coping with evolutionarily novel challenges—these being challenges that occur infrequently throughout the phylogeny of a lineage—thus cannot be solved with recourse to specialized cognitive systems adapted to coping with recurrent problems (such as those associated with cheater detection or language acquisition). Instead, such problems require generalized and open-ended problem-solving systems, such as learning and working memory, in order to tailor solutions to them. The ability to innovate a solution to such a problem (via the development of a tool) is a key manifestation of the action of these generalized problem-solving systems (Geary, 2005).

The Holocene is believed to have been an environment of evolutionary relevance for GCA as it was characterized by major transitions among human populations, especially in Eurasia, away from a hunter-gatherer mode of subsistence, towards a sedentary agriculture-based one and beyond that to urbanization (Cochran and Harpending, 2009; Hawks et al., 2007). This brought with it many evolutionarily novel problems, such as having to cope with increased population densities (Cochran and Harpending, 2009; Hawks et al., 2007). Innovations that played a major role in facilitating this transition would have included the domestication of cultivars and animals and the development of novel tools for raising the productivity of land (such as the plough; Cochran and Harpending, 2009). Cultural innovations such as monotheism, monarchy, aristocracy, feudalism, and currency-based economics arose in response to the need for coping with the hierarchical power distribution characteristic of large, static populations (Cochran and Harpending, 2009). Those populations that were successful in using innovations to solve novel problems would furthermore have had an advantage in warfare, being better able to innovate weapons and tactics, allowing them to replace less successful populations. Population growth would have increased the chances of rare GCA-enhancing mutations arising, which would have favored the aggregate fitness of those populations, permitting them to expand to the greatest extent (Cochran and Harpending, 2009; Darlington, 1970; Hawks et al., 2007). Thus, positive culture-gene co-evolutionary feedback occurring during the Holocene might have rapidly increased GCA (Cochran and Harpending, 2009). Consistent with this, there is evidence that contemporary populations, which historically made the transition from a hunter-gatherer to an agrarian subsistence paradigm, have higher frequencies of working memory-enhancing genetic variants, compared to populations that never or only recently made this transition (Piffer, 2013).

High-resolution genome-wide association studies (GWAS) have identified a number of single nucleotide polymorphisms (SNPs) primarily associated with the development of the central nervous system that predict variance in both GCA and educational attainment—with which GCA shares approximately 60% of its linkage-pruned genetic variance (Davies et al., 2016; Domingue et al., 2015; Okbay et al., 2016; Rietveld et al., 2013, 2014; Selzam et al., 2017). These SNPs can be concatenated into cognitive polygenic scores (henceforth POLY_COG). The term ‘cognitive’ is being used here as a broader term than GCA, as it includes not only the cognitive ability measured by standardized tests but also other sources of variance that contribute to educational attainment. Polygenic scores are more or less normally distributed and can be used as a molecular index of GCA and educational attainment in standard regression-type analyses.

Recent advances in the recovery of ancient DNA have led to the accumulation of a wealth of genomic data on ancient human and other hominid populations living in the Holocene and prior epochs. These data permit comparisons to be made between ancient and modern human populations using POLY_COG, thus permitting the question of whether GCA rose during the Holocene to be tested directly. This analysis will here be conducted via the comparison of an ancient Eurasian genome sample from between 4.56 and 1.21 kyr BP with a modern, ancestrally matched European one.

Materials and Methods

Modern Genomes

Previous studies comparing ancient and modern genomes have utilized the Phase3 1000 Genomes dataset (1000 Genomes Project Consortium, 2015) as a modern reference set, on the basis that it is the most globally representative set of genomes currently available (Allentoft et al., 2015). The 2,504 samples in the Phase3 release are sourced from 26 populations, which can be categorized into five super-populations by continent. These are East Asian (EAS), South Asian (SAS), African (AFR), European (EUR), and American (AMR). These genomes were obtained in the form of VCF files, via the online 1000 Genomes FTP. The EUR sample is comprised of 503 individuals from five different countries (Finland, Great Britain, Italy, Spain, and the United States—specifically, north-western Europeans from Utah).

Ancient Genomes

Data on the ancient Holocene genomes were obtained from the European Nucleotide Archive (Leinonen et al., 2011) in the form of BAM files (accession number PRJEB9021). These genomes, which have been utilized in previous research to analyze the spread of variants associated with skin and eye color, among other traits (Allentoft et al., 2015), were aged between 4.56 and 1.21 kyr (M = 3.44 kyr, SD = 0.62 kyr), covering the early Bronze to early Iron Ages. Their genome-wide coverage is relatively low (0.01–7.4X average depth, overall average = 0.7X; Allentoft et al., 2015). There are 102 genomes in total; however, two samples (RISE507 and RISE508) were sourced from the same individual (Allentoft et al., 2015). The total allele count was therefore calculated employing the full sample, alternately using each version of the genome—the results of the two analyses were then averaged. The majority of the ancient genomes were sourced from regions that are presently considered to be parts of Europe, the remainder being sourced from Central Asia. Late Bronze Age European and Central Asian gene pools resemble present-day Eurasian genetic structure (Allentoft et al., 2015). Consistent with this, the F_st values range from 0 to 0.08, when present-day 1000 Genomes EUR are compared with the various populations comprising the Ancient samples, indicating little to modest genetic differentiation (little differentiation corresponds to an F_st range of 0 to 0.05, and modest to an F_st range of 0.05 to 0.15; Hartl & Clark, 1989). These values are lower than the distance between present-day Europeans and East Asians (F_st = 0.11; Allentoft et al., 2015). The two ancient genomes belonging to the Siberian Okunevo culture (RISE515 and RISE516) were somewhat exceptional however, exhibiting modest differentiation relative to the EUR sample when compared with the sample size and coverage weighted average F_st value computed for the additional 92 ancient genomes for which these data were available (F_st = 0.074 vs. 0.013 for the remainder of the sample). On this basis, they were removed from the final ancient sample in order to reduce the genetic differentiation between the ancient and modern samples. The final sample of ancient genomes was therefore comprised of 99 genomes, sourced from sites located in present-day Armenia (8%), Czech Republic (6%), Denmark (6%), Estonia (1%), Germany (10%), Hungary (10%), Italy (3%), Kazakhstan (1%), Lithuania (1%), Montenegro (2%), Poland (7%), Russia (36%), and Sweden (8%).

Reference Genome

Paired-end reads are aligned to the reference human genome obtained (in the form of a FASTA file) from the UCSC database (GRCh37/hg19).

Read Realignment and Base Recalibration

After removing duplicate reads, the reads were realigned around the known indels from the 1000 Genome sample using the GenomeAnalysisTKLite-2.3-9 (GATK) toolkit. The known indels set was obtained from the GATK resource page (McKenna et al., 2010). After performing realignment, the base recalibration step is performed (DePristo et al., 2011; Van der Auwera et al., 2013). Known variant position is taken into account to recalibrate the quality score.

Variant Calling and Computing POLY_COG

After performing realignment, the GenomeAnalysis- TKLite-2.3-9 toolkit UnifiedGenotyper was used to identify the SNPs and short indels. The positive allele counts (positive here refers to the sign on the β value derived from GWAS regressions—indicating that they increase the trait of interest) and allele frequencies of the target SNPs were calculated from the Ancient Genome VCF files and also from the 1000 Genomes VCF files using a custom-made PERL script (available on request). Variant calling was used to construct the first POLY_COG score from a list of 162 SNPs that were found to have high genome wide significance (p ≤ 5 × 10⁻⁸) as ‘hits’ for educational attainment in a recent large-scale, meta-analytic GWAS study (Okbay et al., 2016). After controlling for ascertainment bias by excluding those variants that were completely absent from the ancient genomes (due to the ancient samples having approximately 30% of the 1000 Genomes coverage) the resultant POLY_COG score consisted of 130 ‘hits’ that were present to some degree in both the ancient and modern genomes.

A second POLY_COG score was constructed using data from another GWAS study of educational attainment and related phenotypes (Davies et al., 2016). In this study, 1,115 SNPs reached GWAS significance, of which 15 were independent (linkage pruned) signals. Four SNPs were absent from the ancient genomes; thus, a POLY_COG was calculated using the remaining 11 SNPs. A third POLY_COG score was constructed utilizing a technique developed in Piffer (2015, 2017; cf. Woodley of Menie, Piffer et al., 2016), employing only the nine GWAS hits from Okbay et al. (2016) that were in close linkage disequilibrium (linkage cut-off r ≥ 0.8, 500 kb linkage window) with ‘hits’ predicting educational attainment across two other large GWAS studies (Davies et al., 2016; Rietveld et al., 2013). Linkage was determined using the NIH LDLink program (Machiela & Chanock, 2015) with the 1000 Genomes Phase3 CEU population as a reference group. These are presented in Table 1.

TABLE 1 Nine Linked Replicated SNPs With Linkage Disequilibrium (D′) and R ² Coefficients

By focusing on only the smaller numbers of linked replicates, the third POLY_COG score was developed as a means of negating the likely high numbers of false positive SNPs present in larger POLY_COG scores. This approach is based on maximizing the signal of directional selection by focusing on the GWAS hits with the highest significance (or replication rate across databases), as opposed to maximizing the amount of variance explained, which instead is the aim of the classical within-population GWAS design. While it is true that the effect of selection on a phenotype is proportional to the additive genetic variance, the selection coefficient for a given SNP is proportional to the amount of variance explained by that SNP. Since the rate of allele frequency change is proportional to the selection coefficient, and because both GWAS data and evolutionary processes are notoriously very noisy, it is important to maximize signal by limiting the analysis to the SNPs with higher GWAS significance, which are also more likely to account for phenotypic variance. This is the same reason why detecting selection on phenotypes that are influenced by a few genes with large effects can be accomplished with traditional, less sensitive tools (e.g., F_st index) than detecting a signal of polygenic selection (Berg & Coop, 2014).

Analysis

All statistical analyses were conducted using R v. 3.3.2 (R Core Team, 2016). In order to test the hypothesis that POLY_COG had different proportions of the positive and negative alleles between the two populations against the null hypothesis that there would be no differences, we used Fisher's Exact and chi-squared tests, since the data are binary (positive vs. negative allele/ancient vs. modern populations). At the same time, the contingency table dealt with unequal representation of genomes for each SNP within the ancient group (i.e., some SNPs were called if they were present in only one individual) by incorporating this information (each SNP is assigned a weight proportional to the number of genomes in which the SNP was called); odds ratios (OR) computed utilizing Fisher's Exact Test were preferred over the chi-squared test because it allows exact estimation of the likelihood of the null hypothesis, rather than relying on an approximation. An alternative to Fisher's Exact test is the log-likelihood ratio G-test, which is employed here as an additional test of significance. This was implemented using the R package DescTools (Signorell, 2016). Additional analyses were also performed on the two smaller POLY_COG (11 and 9 SNPs) using two unconditional tests (considered more powerful alternatives); specifically, Barnard's (1945, 1947) and Boschloo's (1970) Exact tests, utilizing the R package Exact (Calhoun, 2016).

In constructing POLY_COG for conventional GWAS-type analyses, it is common practice to weight each SNP by its associated GWAS regression β value, as these polygenic scores are typically the sum of very large numbers of SNPs with highly heterogeneous effect sizes as predictors of the phenotype of interest (e.g., Okbay et al., 2016). These regression weights were not found to significantly correlate with the frequency differences between modern and ancient populations for the largest of our POLY_COG, however (130 SNP; r = -0.128, p = .07), likely because selecting only the SNPs with the highest p values introduced range restriction among the regression weights. Therefore, given that the SNPs comprising our POLY_COG could be considered to be roughly ‘equally good’ as predictors of the trait of interest, the inclusion of regression β in the analysis as weights is therefore unlikely to add additional power, which in turn justifies the use of a more straightforward and also parsimonious binarization approach in comparing the ancient and modern genomes.

A second analysis was carried out on a subsample of 66 ancient genomes for which sample age was known (age range: 1.21–4.46 kyr). The correlation was computed between individual positive allele counts derived using the 130 SNP POLY_COG with genome age (ascertained by utilizing radiocarbon dating), in order to test for change over time. Since on average the ancient samples had approximately 30% of the 1000 Genomes coverage, carrying out this analysis using the smaller POLY_COG would add too much noise and produce unreliable estimates for individual genomes. Two additional genomes for which age data were present (RISE413 and RISE42) were identified as having significantly outlying values of positive allele count based on the outlier-labeling rule (Hoaglin et al., 1986) because there was only one target SNP present per genome. For the analysis only genomes whose coverage was limited to two or more of the 130 SNPs were retained. Removing these two outlying genomes also reduced the skew of the regression residual from 1.16 to 0.698, making the remaining dataset more amenable to parametric regression.

Random SNP Polygenic Scores

As an additional robustness test, a method developed by Piffer (2017) was employed, in which the Monte Carlo model was run using randomly drawn and non-overlapping sets of SNPs matched for minor allele frequency (MAF). The expectation is that the three POLY_COG should outperform the majority of random polygenic scores in terms of revealing a frequency difference, favoring modern populations, consistent with theory. The random polygenic scores were computed for non-overlapping sets of SNPs for the modern and ancient populations and for each allele (A, B). A 2 × 2 matrix was computed for each comparison (Allele A/Ancient population, Allele A/Modern population, Allele B/Ancient population, Allele B/Modern population). As with the POLY_COG, OR for each matrix were then computed. Finally, the proportion of ORs with a weaker effect (i.e., higher OR value) than what was observed for the three POLY_COG was computed. Frequency matching was carried out using SNPSnap (Pers et al., 2015), by entering the nine linked replicate SNPs, which it must be recalled are likely the closest SNPs to the true causal variants, and setting LD at r ² < 0.1 (for EUR). This resulted in a total of 7,369 SNPs being retrieved.

Results

A 2 × 2 contingency table was created for each of three alternative POLY_COG assigning allele counts by GWAS effect status (positive or negative effect) to ancient and modern genomes. OR of the association between positive or negative alleles and ancient or modern populations were computed and are reported in Table 1. As Fisher's Exact and chi-squared tests revealed nearly identical results, only the former are reported. In addition, the G-test (log-likelihood ratio) is reported as an alternative test of the significance of the difference. As a further robustness test, the difference in POLY_COG between the two sets of genomes was examined using the unconditional Barnard's and Boschloo's Exact tests. These statistics are extremely computationally intensive when large numbers of pair-wise comparisons are employed; therefore, they could only be used for the two smaller POLY_COG (9 and 11 SNPs). The results were identical to those obtained via the conditional techniques.

The Monte Carlo model was run with 6,740 random SNPs (after removal of the ones absent from the ancient genome in order to control for ascertainment bias). Hence, 748, 613, and 52 random SNP polygenic scores (corresponding to the 9, 11, and 130 SNPs POLY_COG, respectively) were computed. The POLY_COG outperformed the random polygenic scores (producing lower OR values) in the majority of instances (9 SNPs = 76.4 % [572 out of 749 draws]; 130 SNPs = 84.6% [44 out of 52 draws]; 11 SNPs = 77.8% [477 out of 613 draws]), which indicates a low probability that the results presented in Table 2 are due to chance.

TABLE 2 Two by Two Contingency Tables With Positive and Negative GWAS Effect Allele Counts for Ancient and Modern Genomes for Each of the Three POLY_COG, Along with Fisher's Exact Test OR and G-Test Log-Likelihood Ratios

Note: The non-integer allele counts result from averaging the genomes of two samples (RISE507 and RISE508) that came from the same individual. The samples had different allele counts likely resulting from different coverage.

These results indicate that modern European genomes have higher POLY_COG relative to ancient ones (sourced from Europe and Central Asia).

A second analysis was conducted to investigate the association between POLY_COG and sample year within the subset of 66 ancient genomes for which radiocarbon age estimates were available (from Allentoft et al., 2015). If selection is operating on these variants throughout the 3.25 kyr covered by this sample, then less ancient genomes should have higher positive allele counts relative to more ancient ones. The correlation of positive allele frequency with the sample year (scaled using the BCE/CE calendar era) was in the expected positive direction (samples from more recent years had higher positive allele counts); the result is furthermore statistically significant when a one-tailed test is used (r = 0.217, 95% CI [-0.026, 0.043], p _one-tailed = .04), which is theoretically justified on the basis that the direction of the effect was anticipated (Kimmel, 1957). The scatter plot is presented in Figure 1.

FIGURE 1 Scatter plot of 130 SNP POLY_COG positive allele counts for 66 ancient genomes as a function of year (scaled in terms of the BCE/CE calendar eras). Note: r = 0.217, p_one-tailed = .04, N = 66. Gray area around trend line corresponds to the 95% CI.

Coverage was not correlated with either genome age (r = 0.022) or positive allele count (r = -0.054). Additionally, a multiple linear regression was run, with positive allele count along with coverage and age as dependent and independent variables respectively. This did not alter the effect size compared to the sample without the outliers for models with and without coverage controlled (β = 0.209 vs. 0.211). The regression models were run with each duplicate genome (RISE508 and RISE507) separately, the resultant β values were then averaged.

Discussion

The Holocene appears to have been an environment of evolutionary relevance for GCA, as modern European genomes have higher POLY_COG relative to those sourced (predominantly) from Bronze Age Europe and Central Asia. A Monte Carlo model employing polygenic scores comprised of non-overlapping sets of random SNPs revealed that the three POLY_COG outperformed between 76.4% and 84.6% of random polygenic scores. It must be noted also that a major strength of our analysis is that it used real allele frequency data, as opposed to traditional Monte Carlo simulation-based approaches, which rely on computer-generated (virtual) data.

A second analysis revealed that within a subset of the ancient genome sample, there was a (one-tailed) significant trend toward increased positive allele counts over time. As was mentioned in the section Introduction, these SNPs are assumed to exert their effects on educational attainment mainly through GCA, but also through other traits relevant for educational attainment. While the increase in these variants over time is certainly consistent with the expectation of rising GCA, the possibility that their increase indicates a simultaneous rise in other factors that make unique contributions to educational attainment (such as ‘slow’ life history or ‘high-K’ social cognitive characteristics; Giosan, 2006) cannot be ruled out.

This increase in POLY_COG can realistically be accounted for in two different ways: (1) microevolutionary selection pressures on standing genetic variation, resulting in soft sweeps (Hermisson & Pennings, 2005; Orr & Betancourt, 2001), and (2) population expansion, replacement, and admixture due to migration.

The idea that selection takes many millennia to produce noticeable changes in a population is at odds with data indicating that within human populations, changes in the means of heritable traits such as height, weight, cholesterol levels, blood pressure, age at first birth, and age at menopause can occur due to selection over a relatively short period of time (Stearns et al., 2010). Selection operating over the course of millennia (as in the present case) would be expected to produce quite considerable microevolutionary change. As was discussed in the section Introduction, Holocene populations appear to have undergone accelerated adaptive microevolution relative to those living in the Pleistocene (Cochran & Harpending, 2009; Hawks et al., 2007). Increasing cultural complexity and technological sophistication among these populations may therefore have arisen in part from selection favoring GCA. Cultural and technological change can in turn create, via culture-gene co-evolutionary feedback, conditions favoring higher GCA (Cochran & Harpending, 2009; Piffer, 2013). This process likely continued until the Late Modern Era, where it has been noted that among Western populations living between the 15th and early 19th centuries, those with higher social status (which shares genetic variance with, and is therefore a proxy for GCA; Trzaskowski et al., 2014) typically produced the most surviving offspring. These in turn tended toward downward social mobility due to intense competition, replacing the reproductively unsuccessful low-status stratum and effectively ‘bootstrapping’ those populations via the application of high levels of skill to solving problems associated with production and industry, eventually leading to the Industrial Revolution in Europe (Clark, 2007, 2014). The millennia-long microevolutionary trend favoring higher GCA not only ceased, but likely went into reverse among European-derived populations living in the 19th century (Lynn, 1996; Lynn & Van Court, 2004), largely in response to factors such as the asymmetric use of birth control and prolonged exposure to education among those with high GCA (Lynn, 1996). Consistent with this, it has been found that various POLY_COG negatively predict reproductive success in contemporary Western populations (Beauchamp, 2016; Conley et al., 2016; Kong et al., 2017; Woodley of Menie, Schwartz et al., 2016). It is important to note that this recent microevolutionary trend (working in the opposite direction) has likely attenuated the difference in POLY_COG between the modern and ancient genomes noted in the present study.

Changes in allele frequencies can also occur via population expansion and replacement. The end of the Last Glacial Maximum initiated major cultural changes, leading to the development of agriculture, population growth, and the subsequent Neolithicization process (8–5 kyr BP), which involved the spread of genes as well as culture from the Middle East (Sokal et al., 1991). Another big cultural change associated with major population movements was the Bronze Age, whose culture started replacing the Neolithic farming cultures in temperate Eastern Europe 5 kyr BP (Allentoft et al., 2015). This Pontic-Caspian Steppe genetic component among contemporary Europeans is associated with the Corded Ware and Yamnaya peoples, and possibly accounts for the spread of Indo-European languages (Allentoft et al., 2015). Being potentially highly fitness enhancing, admixture among populations may have furthermore led to the spread of the causal variants that are tapped by POLY_COG. However, recent evidence suggests that the European and Central Asian gene pools toward the end of the Bronze Age closely mirror present-day West Eurasian genetic structure (Allentoft et al., 2015). Hence, population movements and replacements likely played a larger role in the evolution of GCA prior to the period covered by the present sample. Sicily and Sardinia, and the southern fringes of Europe in general are notable exceptions however, as they tend to cluster genetically with Neolithic farmers (Allentoft et al., 2015).

The first theoretical scenario (microevolutionary change involving standing genetic variation) most likely therefore played the largest role in the POLY_COG increase during the time period covered by our samples.

It has been noted that concomitant in time with the apparent increase in POLY_COG observed presently is an apparent parallel selection trend favoring smaller brains (Hawks, 2011). These trends may at first appear to be contradictory, as GCA is positively associated with brain volume. It must be noted, however, that meta-analysis reveals that the phenotypic association is weak (ρ = .24; Pietschnig et al., 2015) with the genetic correlation being slightly stronger (r _g ~0.30; Okbay et al., 2016), indicating that the majority of the phenotypic and genotypic variance in brain volume is unrelated to variance in GCA. Furthermore, studies utilizing POLY_COG have found that it does not predict variation in brain volume, despite both POLY_COG and brain volume making independent contributions to GCA (Deary et al., 2016). These findings suggest that declining brain volume during the Holocene may have been a consequence of enhanced brain efficiency stemming from increased corticalization and neuronal connectivity, with more bioenergetically optimized brains simply requiring less mass to achieve greater processing power. It may therefore be this process that the increase in POLY_COG is tracking in our samples. A prediction stemming from Hawks (2011) is that variants that predict brain volume and not GCA should show the opposite trend in time to POLY_COG when examined in the context of the present samples.

A potential limitation of the present work is linkage disequilibrium (LD) decay. This takes place when a pair of SNPs on a chromosome in a population moves from linkage disequilibrium to linkage equilibrium over time, due to recombination events eventually occurring between every possible point on the chromosome (Bush & Moore, 2012). LD decay can be an issue when comparing polygenic scores computed using genetically different populations (i.e., either across time or space). Since most GWAS hits are actually tag SNPs, decay in LD implies that the causal SNPs will be less efficiently flagged by the tag SNPs, resulting in the tag SNPs coming to resemble a sample of random SNPs. Since the frequency of the average SNP allele is 50%, the GWAS hits will tend to converge towards an average frequency of 50%, with increasing LD decay. The implication of this for our analysis is that our estimates of the difference in POLY_COG between the reference population (modern European) and the ancient population are therefore conservative when the mean frequency is lower than 50%, while the opposite is true when the frequency is higher (the mean frequencies in the modern 1000 Genomes samples were the following: for the 130 SNP POLY_COG = 47%, 9 SNP POLY_COG = 44%, and the 11 SNP POLY_COG = 53%). Furthermore, because GCA likely involves around 10,000 causal variants (Hsu, 2014), and as there is relatively little genetic differentiation between the ancient and modern genomes, as indicated by the sample size and coverage weighted F_st value of 0.013, the results of simulations indicate that LD decay is ultimately expected to be minimal (Scutari et al., 2016) and should not be confounding the results.

The present effort constitutes a direct attempt to show the theoretically anticipated change in POLY_COG over time utilizing relatively simple methodology, which makes few modeling assumptions. The difference in POLY_COG between the ancient and modern genomes is indifferent to the use of alternate versions of POLY_COG and to the use of various statistical techniques for quantifying the significance of the group difference. The results of the Monte Carlo model indicate that the majority of frequency-matched random SNP polygenic scores underperform relative to the POLY_COG. Among a subsample of the ancient genomes, indications of increasing positive allele counts as a function of time are also present. It is acknowledged that there are limitations in the data employed—larger samples of ancient genomes would certainly increase confidence in the present finding. Furthermore, it is acknowledged that any of the discrete lines of evidence that have been brought to bear on the question can, in isolation, of course be criticized, with some lines being considered weaker than others. It must be noted, however, that multiple operationalizations of the variable of interest and multiple methods have been employed here, which, even if each individual line of evidence is to be considered ‘weak’, should nevertheless greatly increase confidence in our conclusion when the evidence is considered collectively as part of a larger nomological net work (Cronbach & Meehl, 1955).