Searching for ancient balanced polymorphisms shared between Neanderthals and Modern Humans

Abstract Hominin evolution is characterized by adaptive solutions often rooted in behavioral and cognitive changes. If balancing selection had an important and long-lasting impact on the evolution of these traits, it can be hypothesized that genes associated with them should carry an excess of shared polymorphisms (trans- SNPs) across recent Homo species. In this study, we investigate the role of balancing selection in human evolution using available exomes from modern (Homo sapiens) and archaic humans (H. neanderthalensis and Denisovan) for an excess of trans-SNP in two gene sets: one associated with the immune system (IMMS) and another one with behavioral system (BEHS). We identified a significant excess of trans-SNPs in IMMS (N=547), of which six of these located within genes previously associated with schizophrenia. No excess of trans-SNPs was found in BEHS, but five genes in this system harbor potential signals for balancing selection and are associated with psychiatric or neurodevelopmental disorders. Our approach evidenced recent Homo trans-SNPs that have been previously implicated in psychiatric diseases such as schizophrenia, suggesting that a genetic repertoire common to the immune and behavioral systems could have been maintained by balancing selection starting before the split between archaic and modern humans.


Introduction
Many of the fundamental processes at the core of complex human behavior and cognitive abilities, including sensory processing, recognition of self and others, emotions, motivation, learning, memory, attention, vocalization and speech processing, executive function, as well as neural development, may be characterized by at least some degree of heritability (Vallender, 2011;Gilissen et al., 2014;Vissers et al., 2015;Johnson et al., 2016).Assuming that genes are partly responsible for these phenotypes, genetic variation between and within species can be expected to give rise to the large behavioral repertoire observed in nature (Bendesky and Bargmann, 2011).Furthermore, the diversity of these traits can be expected to be shaped by the fundamental forces of microevolution, including genetic drift, directional natural selection (either positive or negative), and balancing selection.
While directional selection tends to reduce variability close to the selected site (Lachance and Tishkoff, 2013), balancing selection results in the persistence of variation in the population or species, even in the face of loss due to drift, leading to an excess of polymorphisms with intermediate frequencies (Nielsen et al., 2009) and increasing genetic diversity around the site of selection (Charlesworth, 2006;Fijarczyk and Babik, 2015).Balancing selection can result from different processes, such as heterozygote advantage (overdominance), negative frequency-dependent selection, and heterogeneity in selective advantage across time or space -all possibly acting in changing environments requiring a fast rate of adaptation (Boon et al., 2007;Wolf et al., 2007;Bendesky and Bargmann, 2011;Pruitt and Riechert, 2011;Schaschl et al., 2015;Taub and Page, 2016).
Send correspondence to Maria Cátira Bortolini.Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Caixa Postal 15053, 92501-970 Porto Alegre, RS, Brazil.E-mail: maria.bortolini@ufrgs.brWhile many cases of directional selection have been reported in the literature, only a handful of examples of balancing selection have been described.This may be due to several factors.On the one hand, balancing selection may be a transient state, leaving marks so subtle that their detection may be difficult using current tests, leading to a large number of false negatives (Fijarczyk and Babik, 2015).On the other hand, most methods rely on the fact that balancing selection will lead to a decreased inter-population diversity.In populations that diverged not long ago, or that are experiencing some level of admixture, however, this pattern will be found at most neutral loci, thereby leading to a large number of false positives.One way of circumventing this issue is by comparison of different species, considering ancient balancing selection, in which most loci should indicate moderate to high divergence, and thus the number of false positives is expected to be smaller.
One of the best-studied targets of balancing selection is the major histocompatibility complex (MHC), which includes many examples of long-term maintenance of transspecies polymorphisms (trans-SNPs), i.e. ancient polymorphisms that survived in derived taxa (Takahata and Nei, 1990;Clark, 1997;Grimsley et al., 1998;Ségurel et al., 2013;Azevedo et al., 2015).The study of these trans-SNPs has revealed a common theme, where individuals heterozygous for genes with key roles in the immune system seem to be more effective in their defense against pathogens, while at the same time presenting only a moderate inflammatory response (Cagliani et al., 2008;Leffler et al., 2013;Azevedo et al., 2015;Teixeira et al., 2015).
Because of the expected loss of shared polymorphic sites over time due to genetic drift, polymorphisms shared between species that diverged a long time ago are rare under neutrality (Clark, 1997).Trans-SNPs common to species separated by a relatively deep evolutionary split, and without recent admixture, are therefore probably adaptive and maintained by balancing selection.Examples of such adaptive trans-SNPs were reported by Leffler et al. (2013) and Teixeira et al. (2015) in their comparison between humans and the Pan genus, which are thought to have diverged about 8 million years ago (Moorjani et al., 2016).In more recently-diverged species, the presence of trans-SNPs must be interpreted with greater caution.For instance, for humans, which have an estimated effective population size (Ne) of ~10,000 individuals, 1% of neutral trans-SNPs will be preserved in the genome even after 53,000 generations (~1,6 million years) (Clark, 1997).Assuming that the split between Homo sapiens and Homo neanderthalensis occurred around 400,000-275,000 years ago (Endicott et al., 2010;Prüfer et al., 2014), some trans-SNPs occurring in these species are expected to be neutral and occur due to stochastic events.Additionally, the retention of ancestral polymorphisms can be due to introgression, the exchange of genetic material between different species due to hybrid-ization (Fijarczyk and Babik, 2015).This has been described for the hybridization between archaic humans (including H. neanderthalensis and the closely related Denisovans) and some modern human populations (Green et al., 2010;Reich et al., 2010;Meyer et al., 2012).
Despite these difficulties, the investigation of trans-SNPs in the genus Homo is an exciting research goal.Hominin evolution is characterized by adaptive solutions rooted in behavioral and cognitive changes.For instance, creative thinking promotes change from prevailing modes of thought or expression, a change that can be associated with fitness gain (Nettle, 2006).In addition to the benefits of change over time, there are adaptive advantages to the parallel maintenance of different behavioral strategies within a species (Sih et al., 2004;Korte et al., 2005;Cagliani et al., 2009;Taub and Page, 2016).Assuming a genetic basis for these behavioral strategies, their parallel persistence can be seen as the result of balancing selection.In support of the idea for balancing selection, there have been several reports of polymorphisms in genes with known roles in modulating complex behavior in modern human and other mammals, which have likely persisted through balancing selection (Cagliani et al., 2009;Schaschl et al., 2015;Taub and Page, 2016).
Expanding on this idea, it is of great interest to investigate the role of balancing selection in the evolution of hominin, including human, behavior on a greater scale.If balancing selection has indeed had an important and longlasting impact on the evolution of behavior in hominins, it can be hypothesized that genes associated with behavior should carry an excess of trans-SNPs across hominin species.Based on this hypothesis, we investigated the role of balancing selection in the evolution of behavior in hominins by studying the pattern of trans-SNPs in genes relevant to these traits in Neanderthals and modern humans.
Importantly, as previously mentioned, many available methods used to detect balancing selection were originally designed to target balancing selection operating over more than 4Ne generations (Clark, 1997).Due to the relatively recent split between modern humans and Neanderthals, these methods are ineffective in detecting balancing selection when studying these two species.To overcome this limitation, we implemented an approach that enables the identification of an excess of trans-SNPs in groups of genes of interest in comparison to the exome background (i.e. the null distribution), which serves as a control for demographic effects, while we added control for gene length biases, GC content (and thus indirectly mutation rate), number of polymorphisms per gene and background selection.Using this approach, we were able to identify polymorphisms shared between modern humans and Neanderthals, many of which located in genes related to immunology and a few in genes playing a potential role in behavior, including genes that may contribute to personality traits and psychiatric disorders.

Defining gene sets
To find genes underlying immune and behavioral systems, we defined two target gene sets, which we named IMMS (genes related to the immune system) and BEHS (genes related to behavior).We populated these gene sets by searching the AmiGO database (http://amigo.geneontology.org/cgi-bin/amigo/browse.cgi) using GeneOntology terms directly related to immune system (for IMMS) and behavior (for BEHS; supplementary material Table S1; http://amigo.geneontology.org/cgi-bin/amigo/browse.cgi).We further added to the BEHS gene set, genes associated with autism spectrum disorder (Iossifov et al., 2014;Yuen et al., 2015), schizophrenia (Carrera et al., 2009;Li et al., 2015;Srinivasan et al., 2015), major depression (Cagliani et al., 2009), and finally the OMIM database (https://www.omim.org)was also consulted for psychopathology associated genes (i.e.schizophrenia, major depressive disorder, autism spectrum disorder, asperger syndrome, attention-deficit disorder, antisocial behavior, and obsessive-compulsive disorder) expanding our total dataset.We then excluded genes that were not available for Neanderthal exome analysis, making for a final count of 1,780 genes in IMMS and 278 in BEHS, with a total of 17,246 analyzed genes including target (IMMS or BEHS, accordingly) and control genes.

Genetic datasets
The high quality exomes of three Neanderthals were retrieved from the Max Planck Institute database (http://cdna.eva.mpg.de/neandertal/exomes/;Castellano et al., 2014).Modern human exome data were obtained from phase 3 of the 1000 Genomes Project (The 1000 Genomes Project Consortium, 2015).To avoid any confounding effects due to interbreeding among archaic humans and modern Eurasians (as reported by Green et al., 2010;Reich et al., 2010;Condemi et al., 2013) in our analysis of balancing selection, we used only Yoruba genomes, as they are assumed to have no admixture history with Neanderthals or Denisovans.We included only autosomal single nucleotide biallelic loci, therefore excluding insertions and deletions (InDels), polymorphic sites with more than two alleles, variants on the sex chromosomes and mitochondrial variants.Only loci found to be heterozygous in the Neanderthal exomes were considered for shared and non-shared polymorphisms; this observed heterozygosity at the individual level was assumed to reflect population-wide polymorphism.Regarding the Yoruba (YRI) sample of 108 individuals, we only considered polymorphisms that both ancestral and derived alleles were segregating in this population.Importantly, as these Neanderthal samples from Croatia and Spain dated to more than 40 Kya (Castellano et al., 2014), we do not expect that any Neanderthal polymorphism is a result of modern human introgression.
Due to the spontaneous deamination of 5-methylcytosine, methylated CpG sites are more prone to mutation than other sites, which raises the probability of allelic identity by state rather than by descent (Azevedo et al., 2015).Because of this, we performed all analyses both including and excluding SNPs located within putatively methylated CpG sites (similarly to Leffler et al., 2013).
The SIFT4G Software (Ng and Henikoff, 2003;Vaser et al., 2016) was used to classify SNPs into synonymous vs. nonsynonymous substitutions and to perform a phenotype prediction for the disruption effect of mutations, allied with known references in literature.Polymorphisms within untranslated regions (UTRs) were excluded from all further analyses.

Evaluating the excess of trans-SNPs in gene sets
We searched in each gene for polymorphisms shared between Neanderthal and Yoruba exomes (i.e.trans-SNPs).We then compared the number of trans-SNPs in each one of the two target gene sets (IMMS and IBHS) to that of the 10,000 gene sets made according to random permutations of all remaining human genes for which Neanderthal exome sequences were available (total of 17,246 genes).In doing so, we always removed genes already accounted for in the target gene set accordingly.For instance, IMMS has 1,780 genes in its dataset, therefore 10,000 random gene sets were built using 15,466 genes as control.Each of the 10,000 random gene sets consisted of as many genes as the target gene set it was simulating, namely 1,780 genes for the comparison to IMMS and 278 to BEHS.The rationale behind the construction of random gene sets, was to generate a null genomic distribution for trans-SNPs, against which each target gene set was then tested.Statistical significance was determined by assessing the deviation in the number of trans-SNPs in the target gene sets in comparison to the background genomic distribution of trans-SNPs generated with the random gene sets.The bash script was used to run this analysis (https://github.com/cegamorim/excess_transSNPs).
Because all loci in each genome were subject to the same demographic history, this approach implicitly controls for demography.However, it does not automatically control for possible effects of background selection, varying mutation rates across sites, and gene size.These effects are known to affect genetic diversity and may therefore bias our results.To control for these effects, genes in the control sets were matched, on a gene-by-gene basis, to those in the IMMS and BEHS target gene sets for background selection, gene length and GC content as follows: Background selection was measured with the B statistic developed by Mcvicker et al. (2009), which was computed based on the expected reduction in nucleotide diversity at a neutral site due to purifying selection at other sites, as a function of re-Balanced polymorphisms between Neanderthals and Modern Humans combination rates, selected site locations, deleterious mutation rate, and the distribution of selection strengths, and indicates the expected fraction of neutral diversity that is present at a given site.A value close to 0 represents near complete removal of diversity as a result of selection, while a value close to 1 indicates that selection has had little effect on diversity (McVicker et al., 2009).To be matched with a target gene, the value of B for a control gene needed to be within 0.1 units of the value of B for the target gene.Gene size was measured as total exonic length, in accordance with the USCS build hg19 refGene table (https://genome.ucsc.edu/).To be matched, the length of target and control genes needed to be within 400 bp of each other.GC content was calculated considering the gene coordinates described in the refGene table of UCSC build hg19.In a first step, we used BEDTools (Quinlan and Hall 2010) to extract the coding exon sequence based on these coordinates, and then used in-house scripts to determine the GC percentage.To be matched, the GC percentages of target and control genes needed to be within 5% of each other.The criterion for thresholds applied was chosen after many trials where at least one control gene in the exome was found for at least 98% of the target genes in IMMS and BEHS datasets.Those target genes that could not be matched to at least one control gene in the exome were excluded from all further analyses in both the target and control gene sets (Table S2).
All data were handled with vcftools 0.1.13(Danecek et al., 2011) and bcftools, as well as using in-house Python and bash scripts.Plots and analyses were implemented in the R environment (www.r-project.org).

Population genetics analyses
In addition to the analysis of trans-SNPs, we considered the classical population genetics statistics Tajima's D and Fst as potential markers of balancing selection, indicated by an excess of polymorphisms with low population differentiation.While shared polymorphisms can detect long-term balancing selection, Tajima's D highlights regions with an excess of (not necessarily shared) polymorphisms, due to balancing selection or population size change.On the other hand, Fst estimates genetic differentiation among populations.Both, Tajima's D and Fst are of interest, since together they indicate potential targets for balancing selection over an intermediate time span, as evaluated using the interval between 0.4 Ne and 4 Ne (Fijarczyk and Babik, 2015).Tajima's D for the YRI population and intercontinental Fst scores (Global Fst) for Yoruba versus Europeans (CEU) and Asians (ASN) were obtained from the 1000 Genomes Selection Browser 1.0 (http://hsb.upf.edu/)(Pybus et al., 2014).It should be stressed that, despite the use of just H. sapiens populations, the polymorphisms selected were those shared with Neanderthals.Negative Fst values were interpreted as 0. Furthermore, since intermediate allele frequencies are a hallmark of balancing selection (Nielsen et al., 2009), we retrieved the average heterozygosity and standard error for all trans-SNPs from the dbSNP database (https://www.ncbi.nlm.nih.gov/projects/SNP/).A threshold for Fst of 0.04, as used by Cagliani et al. (2013) to represent low values for Fst among human populations, and an average heterozygosity greater than 0.400 (as employed by Pakstis et al., 2010) were considered to flag up SNPs potentially affected by balancing selection.
All relevant SNPs identified in our analyses were queried for known associations with psychiatric disease using the GWAS catalog implemented in the UCSC Table Browser and the available literature.

Results
We retrieved 22,832 heterozygous sites from the Neanderthal exome (Castellano et al., 2014), which 4,117 were trans-SNPs (Table S3) found within 2,519 genes.Here trans-SNPs are defined as heterozygous sites in Neanderthals that were also polymorphic in modern humans, represented by Yoruba people, according to the 1000 Genomes phase 3 data (The 1000 Genomes Project Consortium, 2015).We note that when loci located within CpG sites are included the numbers of such polymorphisms and genes almost doubles (Table S4), which is expected due to the high mutability of such sites (Kong et al., 2012), some of which will potentially not be trans-SNPs but only identical by state.
We sought to evaluate if gene sets related to the immune system (IMMS; 1,754 genes) and behavior (BEHS; 271 genes) were enriched for these trans-SNPs in comparison to the genome as a whole.To do so, we built and permutated 10,000 random sets of genes to be compared to these, by matching each gene in these two sets of genes to others in the genome, controlling for background selection, gene size, and mutation rate (see Material and Methods), since these factors are known to affect the number of polymorphisms in each gene.Implicitly we were also controlling for demography, since we were comparing target genes to others in the same genome, which were therefore subject to the same demographic history.Moreover, because there is no evidence up to date of interbreeding between archaic and modern humans before H. sapiens migrated out of Africa, we used Yoruba samples to control for the effect of archaic introgression (Green et al., 2010;Reich et al., 2010;Meyer et al., 2012).Below we describe our findings using this approach for each gene set individually (IMMS and BEHS).

Signals of balancing selection in the Immune System
We identified 547 trans-SNPs in IMMS genes.This number was significantly higher than the null distribution for trans-SNP observed in the 10,000 random permutations of genes from the control set (p-value= 0.016); Table 1, Figure 1 and Figure S1; see Methods).This pattern remains even if we exclude from the analysis the CpG transitions, which are known to present a higher mutation rate (Table S3).This pattern is borderline significant when considering only shared non-synonymous substitutions (297 polymorphisms; p-value = 0.05).These results suggest that in the genus Homo, genes underlying immune system function are more likely than non-immune-genes to maintain longterm shared polymorphisms, possibly through balancing selection.Additionally, we found that Neanderthals IMMS genes harbor more heterozygous sites (2,024 polymorphisms; shared or not with modern humans) than the random sets generated by permutation (p-value < 0.001; Figure 1, column "Polymorphisms"; Table S5).Many of these loci have SNP ID numbers (rs) and have thus been found to be at least biallelic in a modern human population.Due to the known hybridization between Neanderthals and some non-African H. sapiens populations, it is difficult to determine whether they represent polymorphisms shared since their split from the common ancestral; nevertheless, these findings reinforce that genes of the immune system maintain a high level of heterozygosity in the genus Homo.
Six of the 547 trans-SNPs shared between Neanderthals and Yoruba that we identified in IMMS (rs2240464, rs56318802, rs5899, rs118014438, rs377657111, and rs14178; Table 2) are located within genes that were previously associated with schizophrenia by the Schizophrenia Working Group of the Psychiatric Genomics Consortium (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014).It is noteworthy that another two of the trans-SNPs in IMMS are found in genes that have been associated with this disorder according to another study: rs28919579 is located in the CD4 gene, a locus that has been linked to schizophrenia through an imbalance of CD4 cell subtypes, and rs374886374 is located within the C4A gene, a potential MHC locus that has recently been associated with schizophrenia (Sekar et al., 2016).C4 is a fundamental element of the classical complement cascade pathway, which rapidly recognizes and eliminates pathogens and cellular debris.In the brain, C4A is expressed in neurons and promotes synaptic pruning, which is impaired in schizophrenic patients (Sekar et al., 2016).The possibility that these trans-SNPs have influenced the behavioral plasticity of Neanderthal and modern humans is specula-    tive, potentially a lot so, but these findings suggest that at least part of the common genetic repertoire that links the IMMS with schizophrenia in modern humans has been maintained polymorphic for thousands of years.We note, however, that immune system genes are known to be highly pleiotropic (Cotsapas et al., 2011;Andreassen et al., 2015;Wang et al., 2015), and this picture may also be true for other traits besides those related to psychiatric disorders.Still, we identified another 60 trans-SNPs located within these 108 loci associated with schizophrenia (Table 2), but that were neither part of our set for IMMS nor BEHS, suggesting that other systems may have trans-SNPs in pleiotropy with schizophrenia.

Potential signals of balancing selection in Behavior System
No excess of trans-SNPs (shared polymorphisms) was found in BEHS genes (51 SNPs in 41 genes; Table 3; Figure 1 and Figure S1).Likewise, and in contrast to the immune system, the number of heterozygous loci found in Neanderthals, shared and non-shared with Yoruba, corresponds to that expected by chance (p-value = 0.1; Figure 1, column "Polymorphisms").Additionally, Tajima's D was not significant for positive values (Table 3).In this regard, significant positive values for Tajima's D would reflect more pairwise differences than segregating sites due to the increased diversity of the region surrounding the selected site, indicating old balancing selection (Nordborg et al., 1996).On the other hand, we found low Fst and high average heterozygosity, both important indicators of potential signals of balancing selection, in five trans-SNPs (rs11176013, rs12628, rs310617, rs438042, and rs362331 (Table 3) shared between Neanderthals and Yoruba.All of these SNPs are associated with psychiatric and neurodevelopmental disorders, including schizophrenia.For instance, rs438042, located near the intron/exon boundary of THBS4 exon 3, is associated with Alzheimer Disease and might be important for splicing, since THBS4 acts in inflammatory responses and synaptogenesis (Cagliani et al., 2013;Cocchi et al., 2015).Another trans-SNP identified in the BEHS gene set, rs310617, located in the EEF1A gene, has been found to be heterozygous in the Denisova specimen.Although the specific role of rs310617 is not known, other substitutions in this gene have been associated with severe intellectual disability and epileptic encephalopathy (Nakajima et al., 2015;Inui et al., 2016).

Discussion
The evolutionary history of hominins is characterized by several notable and peculiar features.Of particular importance to the successful evolutionary trajectory of the genus Homo was the emergence of a large brain capable of sustaining complex and plastic adaptive behaviors, as well as cognitive skills.In recent years, a growing body of research has countered the notion that the Neanderthals were devoid of symbolism and presented lower cognitive abilities and less sophisticated behavior than the early H. sapiens of the Paleolithic (Akazawa et al., 1993;Zilhão and Trinkaus, 2002;Zilhão et al., 2010;Pike et al., 2012;Rendu et al., 2013).This is in line with the findings of recent genetic and paleoneurological research.For instance, research by Mounier et al. (2016) and Ponce- de-León et al. (2016) revealed strong similarities between modern humans and Neanderthals in both endocranial anatomy and general brain development, while a study of 162 genes related to cognition by Paixão-Côrtes et al. (2013) identified a genetic repertoire shared between extinct archaic humans and modern humans.Assuming that the use of cognitive skills and complex behavior as an adaptive strategy represent a central element of the human evolutionary trajectory (Cagliani et al., 2009;Schaschl et al., 2015;Taub and Page, 2016), we sought to evaluate whether natural selection, particularly in the form of balancing selection, played any role in the evolution of genes potentially related to human behavior.Examples of a role for balancing selection in behavioral plasticity have been reported from primates (Babb et al., 2010;Claw et al., 2010;Dobson and Brent, 2013;Goto et al., 2016;Taub and Page, 2016), rodents (Lonn et al., 2017) and even arthropods (Fitzpatrick et al., 2007).To detect balanced polymorphisms in genes related to behavior, we implemented an approach based on the search for an excess of shared polymorphisms (trans-SNPs) between archaic and modern humans, in comparison to the genome as a whole.
Our analyses revealed no excess of trans-SNPs in genes known to underlie behavioral traits (captured in our BEHS gene set, see Material and Methods).These findings are consistent with the current knowledge about the genetics underlying H. sapiens brain function.Genes expressed in the brain have a large number of functions and the interactions between them are complex (giving rise to basal and specific behavioral phenotypes).Therefore, these genes are subject to functional constraint (Nielsen et al., 2009).Recently, Aggarwala and Voight (2016) developed the concept of genic tolerance to assess the probability of nucleotide substitution in the human genome, based on factors such as population history and selection, among others.In their analysis, genes playing a role in neurodevelopmental and psychiatric disorders were found to have a strong genic intolerance to nucleotide substitution.In the context of functional complexity and constraint, relatively few, key genetic changes can lead to larger effects on certain phenotypes, in response to a specific selective pressure, while at the same time maintaining original functions.Among the results of our analyses, five trans-SNPs located within our set of BEHS (rs11176013, rs12628, rs310617, rs438042, and rs362331; Table 3) presented high average heterozygosity and low Fst values, suggesting a homogeneous distribution of both alleles between populations.These results imply that balancing selection did not have a 74 Viscardi et al.   significant role in the evolution of genes implicated in human behavior as a whole, but may have been important for the evolution of particular genes within this set.Alternatively, our approach may not have had enough statistical power to detect the effect of balancing selection on the evolution of human behavior.That could be the case if frequency dependent selection, rather than overdominance, was the main mode of balancing selection, since our test is best suited to detect overdominance.Alternatively, the signals for balancing selection in the BEHS set may significantly pre-or postdate the evolutionary split between Neanderthals and modern humans.
Another possibility that deserves to be discussed in the light of our results comes from recent findings suggesting an association between behavioral traits and genes previously implicated in the immune response (Power et al., 2015;Sekar et al., 2016).Through our inter-Homo exonic trans-SNP approach, we found that genes underlying immune function (found in the IMMS gene set) contain more ancestral polymorphisms than expected by chance in both Neanderthal and modern humans.These genes have played an important immunological defense role during speciation and migration of the genus Homo in a probable similar context of their hominin common ancestral.Our results support the idea that the variability of immune genes is both a target and an outcome of balancing selection (Grimsley et al., 1998;Ségurel et al., 2013).Interestingly, and perhaps surprisingly, several studies have revealed a connection between the genetics of the immune system and human behavior.For instance, some of the strongest genetic associations found for schizophrenia at the population level involve variation in the immune system loci (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014; Sekar et al., 2016).It has been suggested that some proteins of the immune system work to promote synaptic pruning (Sekar et al., 2016), which is impaired in schizophrenic patients (Lips et al., 2012).Other mechanisms involved both in the etiology of schizophrenia and in the immune system have also been suggested, including deviant immune responses (Pandarakalam 2013).Some of the IMMS trans-SNPs identified here have previously been found to be associated with schizophrenia in modern humans (rs2240464, rs56318802, rs5899, rs118014438, rs377657111, rs14178, Table 2).Furthermore, the trans-SNP rs374886374 is located in the gene C4A, at the MHC locus, which has recently also been associated with schizophrenia in a landmark study using several Psychiatric Genomics Consortium cohorts (Sekar et al., 2016;nearly 65,000 individuals;p-value < 10-8).
Schizophrenia is known to have a high heritability of around 60-80%, and interestingly, it is frequently contextualized in hypotheses that attempt to explain the evolution of modern human complex behavior (Srinivasan et al., 2015).One hypothesis aiming to reconcile the relatively high prevalence (~0.5 to 1%) of this disorder across human populations with its negative effect on fitness is that balancing selection is maintaining several alleles at loci contributing to creative thinking, a trait that increases fitness.Under unfavorable circumstances, however, the same alleles are thought to increase vulnerability to psychiatric disorders, including schizophrenia (Power et al., 2015;Srinivasan et al., 2015).Our findings contribute to this hypothesis and suggest that some components of the immune genetic repertoire that were maintained polymorphic in both archaic and modern humans could have indirectly influenced the evolution of human behavior.This would represent an extraordinary case of evolutionary co-option, in which IMMS genes under balancing selection harbor ancestral adaptive polymorphisms related to the behavioral plasticity of the genus Homo.Allied to our conclusion, recent studies have contributed to unveil the physiological process of autoimmunity in cognition, being proposed as the driving force of cognitive evolution in genus Homo (Nataf, 2017).
Finally, a range of other molecular and biological processes certainly play an important role in the evolution of the behavioral plasticity characteristic of Homo species, such as gene regulation and epigenetic mechanisms.Moreover, beyond the role of the heterozygote advantage in maintaining these polymorphisms, other forms of natural selection (frequency-dependent, directional, etc.) at multiple levels (i.e., individual, kin, and/or group levels; Polimeni and Reiss, 2003;Wilson and Hölldobler, 2005;Zhang and Perc, 2016), together with the unequivocal role of culture (Mesoudi, 2016), have shaped and, in the case of our species, continue to shape, human evolution.A full exploration of these topics is well beyond the scope of the present study, which intends only to explore and discuss some genetic and evolutionary pieces of this complex puzzle.Future studies may help to build a more complete picture of the evolution of hominin behavior.

Figure 1 -
Figure 1 -Density distribution of the average number of polymorphisms per gene observed for random sets genes (blue shade) in the Neanderthal samples matched to those included in the immune system (IMMS) and behavioral system (BEHS) target gene sets (red bars).

Table 1 -
Balanced polymorphisms between Neanderthals and Modern Humans 71 Percentiles of the distribution of the mean values of polymorphism per gene in 10 000 random combinations simulating each of the proposed target gene systems (IMMS and BEHS) 1 .
1 Values close to the mean number of polymorphism per gene for each of the target gene sets (IMMS and BEHS) are in italic and underlined, while significant values are in bold.Mean values for the IMMS gene set: Total Neanderthal SNPs = 0.894647, non-synonymous trans-SNPs = 0.115604, trans-SNPs = 0.220957.Mean value for the BeHS gene set: Total Neanderthal SNPs = 1.1875, non-synonymous trans-SNPs = 0.07353, trans-SNPs = 0.169118.

Table 2 -
Trans-SNPs within loci associated with schizophrenia according to the Schizophrenia Working Group of the Psychiatric Genomics Consortium (2014).Loci are sorted by chromosome position.Genes in bold are those belonging to the IMMS target gene set.Genes shaded in gray are those considered as Credible causal schizophrenia SNPs (sets of SNPs that are 99% likely to contain the causal variants) according to the Schizophrenia Working Group of the Psychiatric Genomics Consortium (2014).

Table 3 -
Trans-species polymorphisms identified in genes of the BEHS gene set.