inappropriate immune activity can increase susceptibility to infectious , inflammatory , and autoimmune diseases , the clinical manifestations of which vary considerably between individuals and populations ( brinkworth and barreiro , 2014 , casanova et  al . ,
the contribution of host genetic factors in explaining such heterogeneity is increasingly documented by genome - wide association studies ( gwass ) , which have identified variants , often located in non - coding regions , associated with disease risk ( parkes et  al . , 2013 , schaub et  al . , 2012 ) .
genetic variants exerting regulatory effects on gene expression , known as expression quantitative trait loci ( eqtls ) , have proven to be of significant biomedical interest ( montgomery and dermitzakis , 2011 ) , as they help to establish links between intermediate phenotypes and organismal traits , such as immunity to infection ( fairfax and knight , 2014 ) . while eqtl studies have mostly focused on steady - state expression measurements ( lappalainen et  al .
, 2010 , spielman et  al . , 2007 , stranger et  al . , 2012 ) , recent work has characterized response eqtls in human cells exposed to various immune or infectious challenges ( barreiro et  al . , 2012 ,
alkan et  al . , 2015 , fairfax et  al . , 2014 ,
. however , the extent and genetic determinants of inter - population transcriptional differences upon immune stimulation remain largely unexplored , yet this is critical to increase knowledge on the varying susceptibility to immune disorders observed at the population level . understanding how natural selection has shaped genome variability represents a powerful approach to identify genes that have played a major role in host survival , complementing immunological as well as clinical and epidemiological genetic studies ( casanova et  al .
indeed , microorganisms are among the strongest selective pressures encountered by humans , and multiple host genes and  variants associated with immune functions and diseases are reported to evolve adaptively ( fumagalli and sironi , 2014 , karlsson et  al . , 2014 , quintana - murci and clark , 2013 ) .
furthermore , there is growing evidence that regulatory variants play a major role in population adaptation and contribute to the diversity of complex phenotypes ( fraser , 2013 , pickrell , 2014 , schaub et  al . , 2012 ) . besides the occurrence of new advantageous mutations ,
functional variants can be introduced in human populations through interbreeding with now - extinct lineages ( vattathil and akey , 2015 ) .
recent data showed that 1%6% of the genome of modern eurasians derives from neandertals or denisovans ( prfer et  al .
furthermore , humans appear to  have acquired genetic diversity from archaic hominins at several immune genes , such as hla , tlr1 , or the oas cluster ( abi - rached et  al .
, 2016 , mendez et  al . , 2013 ) . despite these findings , the  impact of selection and archaic admixture on driving population differences in immune responses remains to be investigated .
here , we determined the degree of immune response variation , and of its genetic and evolutionary sources , within and between human populations .
we used rna sequencing ( rna - seq ) to characterize the responses of primary monocytes , from individuals of european and african descent , to various toll - like receptor ( tlr ) ligands and influenza a virus , and we mapped eqtls .
we found that marked differences in immune responses exist between populations due to the contribution of cis- and trans - acting regulatory variants .
we also show that natural selection has contributed to the observed population differences of immune responses and establish that admixture with neandertals introduced regulatory variants affecting responsiveness to immune stimuli into european genomes .
population variation in immune responses was characterized in primary monocytes , as a model of an innate immune cell type , from 200 healthy individuals of self - reported african and european ancestry ( 100 individuals of each population ) ( see star methods ; figure  s1 ) .
cells were exposed , for 6  hr , to ligands activating tlr4 ( bacterial lipopolysaccharide [ lps ] ) and tlr1/2 ( pam3csk4 , a synthetic triacylated lipopeptide ) , responsible principally for sensing bacterial components , tlr7/8 ( r848 , an imidazoquinoline compound ) , responsible predominantly for sensing viral nucleic acids , and to a human seasonal influenza a virus ( iav ) .
rna - seq data were collected from unstimulated and stimulated monocytes , yielding a final dataset of 970 transcriptional profiles with 34 million single - end reads per sample .
high - density genome - wide genotyping and whole - exome sequencing data were generated for all individuals and used to characterize their genetic ancestry , map eqtls , explore patterns of allele - specific expression ( ase ) , and perform population and evolutionary genetic analyses . to assess the transcriptional responses of cd14 monocytes , we first processed and normalized the rna - seq data and accounted for potential batch effects and confounding factors that could impact gene expression variation ( figures s2 and s3 ; see star methods ) .
we obtained a final set of 12,578 expressed genes , 6,752 of which were differentially expressed upon stimulation in at least one condition ( |log2(fold change [ fc])| > 1 ) at a false discovery rate ( fdr )  < 0.05 ( table s1a ) .
using adjusted expression data , principal component ( pc ) analysis revealed that pc1 and pc2 , accounting for 65% of total variation , corresponded primarily to iav infection and tlr activation , respectively ( figure  1a ) .
we next used weighted correlation network analysis ( wgcna ) ( langfelder and horvath , 2008 ) to define modules of genes that  showed similar behaviors ( up-/downregulation ) upon immune stimulation and identified ten modules , each comprising 2574,070 genes ( figures 1b and s4 ) .
the gene modules upregulated upon stimulation ( modules 14 ) were enriched in  diverse gene ontology ( go ) functions related to host defense , including an inflammatory response restricted to tlr activation and a global antiviral response exacerbated upon iav infection ( figure  1c ; table s1b ) .
the gene modules downregulated ( modules 59 ) , or containing similarly expressed genes across conditions ( module 10 ) , were enriched in functions such as rna processing and translational termination ( table s1b ) . using the transcription factor affinity prediction ( trap ) method ( thomas - chollier et  al . , 2011 )
, we found that the annotated promoter regions of the genes within each module were enriched in binding motifs for specific transcription factors , such as nuclear factor b ( nf-b ) , irf1 , and gata2 ( table s1c ) .
these analyses show that cellular context is the major determinant of transcriptional variability and provide a genome - wide view of common and specific responses of cd14 monocytes to activation by tlr ligands and iav infection .
the estimated genetic ancestry of individuals accurately reflected self - reported ethnicity , indicating negligible levels of admixture between the two groups ( figures s5a and s5b ) .
we searched for genes that show population differences in expression ( popdegs ; fdr  < 0.05 ) and further considered the magnitude of such differences by setting different thresholds of fold change between populations ( fcpop ) .
we identified 5,501 popdegs with a |log2(fcpop)| > 0.2 in at least one condition , a figure that dropped to 821 and 70 when increasing the magnitude of fold change ( |log2(fcpop)| > 0.5 and > 1 , respectively ; table s1d ) . among genes displaying the largest population differences ( table 1 )
, we observed the scavenger receptor marco , involved in early inflammatory responses to influenza ( areschoug and gordon , 2009 ) ; the chemokine receptor cx3cr1 , mediating skin wound healing ( ishida et  al . , 2008 ) ; and , more generally , several interferon - stimulated genes .
we next searched for genes presenting population differences in their response to treatment , relative to non - stimulated cells ( popdrgs ) .
we found 3,841 popdrgs ( fdr  < 0.05 , 70% of popdegs ) , the majority of which were treatment specific ( 2,687 popdrgs ; table s1e ) .
popdrgs displaying stronger responses in africans were enriched in go functions from metabolic processes to defense responses , while popdrgs responding more strongly in europeans were essentially restricted to defense functions in the tlr conditions and enriched in translational processes upon iav infection ( table s1f ) .
popdrgs showing the greatest population differences ( |log2(fcpop)| > 1 ) were enriched in cytokines and chemokines ( fisher s exact test , odds ratio [ or ]  = 36.7 , p  < 10 ) , including il12b and csf3 , responding more strongly to pam3csk4 in africans , and ccl8 , ccl13 , ccl15 , ccl23 and cxcl10 , being more responsive to lps in europeans ( table 1 ) .
these results indicate that while population transcriptional differences  of moderate effect are widespread , strong differences predominantly affect antiviral and inflammation - related genes that differ  markedly in responsiveness between africans and europeans .
we next mapped eqtls by testing for associations between 10,278,745 snps ( the set of genotyped and imputed snps presenting a minor allele frequency [ maf ] > 0.05 ) and gene expression phenotypes .
we first mapped local , likely cis - acting eqtls within 1 mb of each gene in africans and europeans separately .
we used an additive linear model ( shabalin , 2012 ) that included the first two pcs of the genetic data ( figures s5c and s5d ) to account for possible population substructure . considering only eqtls having an effect size of |eqtl|
> 0.2 at a fdr of 5% , we found 2,665 genes with an eqtl in at least one condition ( figure  s6a ; table s2a ) .
of these , 917 genes presented a response eqtl ( reqtl ) , an eqtl with a significantly larger effect size after treatment than at the basal state ( |eqtl| > 0 and p  < 10 , figure  2a ) .
, 2014 , lee et  al . , 2014 ) ; most reqtls were treatment specific ( 62% , 570 genes ) , indicating strong context specificity of the genetic regulation of immune responses . to investigate the functional features of ( r)eqtls
, we used the predicted regulatory elements of cd14 monocytes ( zerbino et  al . , 2015 ) and identified a strong enrichment in such elements , particularly in promoter sequences ( or > 10.4 , p  < 10 ; figure  s6b ) .
furthermore , we observed strong enrichments of basal eqtls and reqtls in binding sites for several transcription factors ( tfs ) , including kdm5a and thap1 at the basal state , tbp and stat3 after tlr activation , and stat2 , hmgn3 , and irf1 following r848 and iav treatments ( figure  2b ) , highlighting mediators of cellular responses to immune activation .
we found that popdrgs were enriched in reqtls in all conditions ( or > 2.6 , p  < 10 ) , an enrichment that increased with the magnitude of the population fold change in gene expression ( |log2fcpop| ; figure  2c ) .
this suggests that differences in transcriptional responses between populations are , at least partially , under genetic control . to test this hypothesis
, we evaluated the fraction of population transcriptional differences that can be explained by genetics ( see star methods ) .
we found that ( r)eqtls account on average for 50% of such expression differences and for up to 70% when focusing on ( r)eqtls of strong effect size ( i.e. , fifth quintile ; figure  2d ) .
furthermore , reqtls associated with popdrgs showed a stronger degree of population differentiation ( mean difference in derived allele frequency |daf|  = 0.24 for popdrgs versus 0.16 for non - popdrgs , p  < 2.2   10 ) , suggesting that differences in transcriptional responses are mainly accounted for by population variation in allele frequency of reqtls .
an example is provided by ncf2 , which is downregulated specifically in africans upon tlr activation , due to the higher daf of reqtl rs2274065 , with respect to europeans ( dafafb  = 0.50 versus dafeub  = 0.07 ) ( figure  2e ) .
we next searched for population - specific ( r)eqtls , i.e. , snps present at similar population frequencies ( maf > 0.05 ) but having a regulatory effect in one population only .
we found 16 eqtls presenting significant differences in effect size between populations ( pinteraction  < 0.001 ) , 5 of which were reqtls ( table s2b ) .
for example , rs1051712 was associated with decreased st3gal6 expression upon r848 stimulation in europeans only ( figure  2f ) .
our analyses suggest that while population - specific gene regulation can occur , population differences in immune responses are mostly the result of regulatory variants presenting different allele frequencies between africans and europeans . to provide a more accurate evaluation of cis effects affecting immune response variation , we mapped allele - specific expression qtls ( aseqtls ) ( figures s6c and s6d ) .
aseqtl mapping is constrained by not only the availability of heterozygotes and read depth but also by effect size , which strongly impacts the power of detection ( figure  s6e ) . to ensure sufficient power , we focused on the 233 genes with large - effect eqtls ( |eqtl| > 0.5 ) that could be tested and found 200 with an aseqtl ( 86% ) , including 160 assessed with high confidence ( paseqtl  < 10 ) ( figure  s6f ; table s2c ) . similarly , among the 42 reqtl genes that could be tested
, we detected 33 ( 78% ) with a stimulus - induced allele - specific response qtl ( asrqtl ) , including 20 assessed with high confidence ( figures 3a and s6 g ; table s2d ) . among these
, we found the tlr - induced ncf2 and pcid2 and the iav - induced arl5b ( figure  3b ) , which regulates the rig - i - like receptor mda5 ( kitai et  al . , 2015 ) .
we next assessed the contribution of common regulatory variants ( maf > 0.05 ) to ase at the individual level . out of 5,889 genes
for which ase could be tested , we identified 1,942 genes with at least one ase event ( |log2(nalternative / nreference)| > 0.2 , fdr  < 0.05 ) , yielding an average of 188 ase events per individual ( table s2e ) .
of these , 275 genes presented evidence of allele - specific responses ( i.e. , significant differences in ase before and after stimulation ) , suggesting g  e interactions ( table s2f ) . focusing on the 160 aseqtls detected at the population level , we consistently observed stronger allelic imbalance in heterozygous individuals , 70% of whom displayed ase in both the presence and absence of stimulation ( p  < 2.2   10 ; figure  3c ) .
our results indicate that , upon immune stimulation , a large fraction of ase events can be accounted for by common regulatory variants , as shown for steady - state expression ( battle et  al . , 2014 ,
finally , we evaluated whether rare coding variants , presenting a frequency 1% and characterized through whole - exome sequencing , impact ase upon immune stimulation . a significant increase in ase
was observed in individuals carrying rare missense variants in stimulated conditions ( or  = 1.34 , p  < 5.0   10 ; figure  3d ) .
notably , nonsense variants contributed to the strong increase in ase in both basal and stimulated states ( or  = 6.8 , p  < 5.9   10 and or  = 10.6 , p 
this is consistent with a role of rare coding variants in the generation of allelic imbalance in monocytes , particularly premature stop variants , possibly through nonsense - mediated decay , as reported for other cell types and tissues ( kukurba et  al . , 2014 , lappalainen et  al . , 2013 ,
2012 ) . besides the contribution of common regulatory variants and rare coding mutations , our results identified a fraction of ase events that are not explained by nearby eqtls ( i.e. , 17% of homozygotes display ase ) .
this suggests the occurrence of secondary mechanisms regulating ase , including undetected eqtls of small effect size or epigenetic effects . to detect master regulators underlying population differences in immune responses , we mapped trans - eqtls , i.e. , snps regulating gene networks over long distances .
our genome - wide mapping across stimulations , correcting for multiple testing , resolved a total of 42 trans - eqtls regulating 165 genes at an fdr of 5% ( p  < 2.7   10 ; figure  4a ; table s3a ) .
of these , 62% ( 103 genes ) were trans - regulated in one condition only , highlighting the high degree of context specificity .
we assessed the contribution of trans regulation to population differences in immune responses , and found that trans - regulated genes , upon tlr4 and tlr1/2 treatments , were strongly enriched in popdrgs ( or > 10.3 , p  < 1.6   10 ; figure  4b ) . to decrease the multiple testing burden of detecting trans - associations
, we further interrogated the 42 trans - eqtls on a single snp basis ( see star methods ) .
this enabled the detection of 794 trans - regulated genes ( p  < 4.4   10 , bonferroni - corrected p  < 0.05 ) , the large majority of which ( 98% ) were associated to a single trans - eqtl .
furthermore , we observed that only two loci , ifnb1 and tlr1 , account for 88% of these associations ( figure  4a ; table s3a ) .
the ifnb1 locus , previously reported upon lps treatment for 24  hr in europeans ( fairfax et  al . , 2014 ) , was the strongest trans - regulatory hotspot .
we found that this locus controlled , in both populations , a tlr4- and tlr1/2-mediated antiviral gene network ( table s3b ) , corresponding mostly to genes belonging to module 2 ( 96% of overlap ) .
genes in this network were enriched in popdrgs ( or  > 9.2 , p  < 10 ) , owing to population differences in ifnb1 response .
local ifnb1 regulatory variants had similar population frequencies ( maximum |daf|  = 0.1 ) and explained only up to 9% of the differences in ifnb1 response .
thus , the population differences observed for ifnb1 trans - regulated genes are not due to variation in the cis regulation of ifnb1 itself but instead are due to yet - unidentified genetic and non - genetic factors .
we identified a pam3csk4-induced gene network that is trans - regulated by the tlr1 missense variant rs5743618 ( i602s ) .
this european - specific trans - eqtl ( dafeub  = 0.71 , dafafb  = 0.01 ) was also associated with the expression of one of the largest networks ( 432 genes , bonferroni - corrected p  < 0.05 ; figures 4a , 4c , and 4d ) .
genes downregulated by the rs5743618-derived variant were enriched in genes belonging to module 1 ( 67% of genes ) and module 4 ( 18% of genes ) ( or > 11.5 , p  <  2.5   10 ; figure  1b ) .
these genes were preferentially involved in responses to bacterial infection ( or  = 6.3 , p  = 6.5   10 ; table s3c ) and included regulators of inflammation , such as ccl5 , il10 , il12b , and ptgs2 ( figure  4d ) . among upregulated genes , signaling - related functions were overrepresented and included bid , ikbke , and pak1 , involved in tnfr1 signaling ( figure  4d ) .
remarkably , tlr1 trans - associated genes displayed strong enrichment in popdrgs ( or  = 8.6 , p  <  10 ) , and , contrary to ifnb1 , such population transcriptional differences were largely explained by genetics ( i.e. , the rs5743618 variant ; figure  4e ) .
we then investigated the effects of this trans - regulatory variant on the inflammatory response , as we previously showed its functional impact on nf-b activity ( barreiro et  al . , 2009 ) , and assessed the correlation of rs5743618 genotypes with the expression of the 81 inflammatory response genes of module 1 ( figures 1b and 1c ; table s1b ) .
the derived c allele ( 602s ) was associated with a significant overall decrease in the expression of inflammatory response genes ( p  = 1.2   10 ; figure  4f ) .
these results reveal major population differences in tlr1/2-mediated responses , which are largely explained by a european - specific tlr1 trans - regulatory hotspot that contributes significantly to differences in the strength of the inflammatory response between africans and europeans .
we next assessed how natural selection , as opposed to genetic drift , has contributed to differences in immune responses between populations .
we computed two metrics  fst , based on the degree of population differentiation ( holsinger and weir , 2009 ) , and ihs , based on haplotype homozygosity ( voight et  al .
, 2006)to detect signals of old and recent events of positive selection , respectively . after matching for maf and linkage disequilibrium ( ld ) patterns
, we found that basal eqtls and reqtls were enriched in stronger values of fst ( p  < 0.005 for eqtls and p  < 1   10 for reqtls , respectively ) and ihs
( p  <  0.002 and p  < 1   10 , respectively ) , relative to genome - wide expectations , in africans and europeans ( figure  5a ) .
significant enrichments in selection signals were also obtained using a composite selection score ( css ) combining fst and ihs , which detects signals of recent , strong positive selection , and the xp - clr method , which uses allele frequency differentiation at linked loci to detect selective sweeps ( chen et  al . , 2010 ) .
among reqtls , the strongest enrichments were observed for the iav condition in both africans ( ihs p  = 0.04 , xp - clr p  < 10 ) and europeans ( fst p  < 10 , css p  = 0.002 ) ( table s4a ) .
this supports a history of positive selection targeting immune - responsive regulatory variants , particularly those involved in responses to viral infection . to highlight specific ( r)eqtl candidates that may have participated in population adaptation at different timescales , we considered loci presenting extreme values of fst or ihs at the genome - wide level ( > 99 percentile ; tables s4b and s4c ) . among these
, we conservatively retained ( r)eqtls that were located in genomic regions presenting a significant enrichment in selection signals and deviated from neutral expectations based on validated demographic scenarios ( see star methods ) . among local eqtls , the strongest signal detected by fst involved an eqtl associated with stronger expression of the methyltransferase gene pcmtd1 in europeans ( fst  = 0.8 , ihs  = 3 ) , while the strongest signal of ihs involved a reqtl associated with a reduced expression of ccr1 following tlr1/2 activation in africans ( ihs  = 4 ; fst  = 0.4 ) ( figure  5b ) .
with respect to trans - eqtls , the master regulatory snp rs5743618 at tlr1 also presented a strong signal of local adaptation in europeans ( fst  = 0.7 , ihs  = 1.5 , pempirical - fst  = 0.002 , psim - fst  = 0.007 ) .
together , our results provide genome - wide support for the important role of regulatory variants affecting basal gene expression and responses to immune stimuli in driving human adaptation .
this , together with the enrichments in genes showing differential expression between populations ( or  = 2.1 , p  = 1.1   10 ) among ( r)eqtls with selection signatures ( tables s4b and s4c ) , emphasizes the contribution of natural selection to the differences in immune responses observed between human populations .
we investigated the impact of admixture between neandertals and the ancestors of europeans on genome - wide expression profiles ( see star methods ) .
we first defined a set of 197,959 variants as of putative neandertal ancestry ( archaic snps [ asnps ] ) if the neandertal allele was present in europeans and absent in africans and located in genomic regions with a high probability of neandertal ancestry ( sankararaman et  al . , 2014 ) .
we identified a total of 52 loci harboring at least one asnp overlapping a local eqtl ( archaic eqtl ) .
interestingly , relative to genome - wide expectations , an enrichment in asnps was observed for basal eqtls ( p  < 0.003 ) and reqtls in r848 and iav conditions ( p  <  0.014 and p  < 10 , respectively ; figure  6a ) . to identify archaic eqtls with high - confidence , we next focused on those located in haplotypes longer than expected under a scenario of incomplete lineage sorting ( figure  s7a ) . among the 19 eqtls presenting strong evidence of neandertal origin ( table s5 ) , 9 corresponded to r848- and iav - induced reqtls , implicating genes encoding ras gtpases such as rab3ip and rapgef3 .
some of these ( r)eqtls carry archaic alleles that are at appreciable frequencies in europeans , suggestive of adaptive introgression . to test this hypothesis , we reasoned that an archaic allele that introgressed into europeans and east asians and was advantageous in one population only should present today unusually high levels of genetic differentiation , relative to genome - wide expectations ( vernot and akey , 2014 ) . when comparing fst between europeans and east asians at archaic ( r)eqtls against the genome - wide distribution of asnps ( table s5 ) , we identified a haplotype that regulates the response of pnma1 to r848 and iav as a significant genomic outlier ( fst  =  0.28 ; pemp  = 0.01 ; figures 6b and s7b  s7e ) .
this archaic haplotype is present at very high frequency in europeans ( 33.5% ) , while it is absent in east asians ( figure  s7b ) .
using simulations that make conservative assumptions about the past frequency spectrum of archaic alleles ( see star methods ) , we found that the high frequency of the pnma1 haplotype in europeans is not compatible with neutral expectations ( psim  < 0.05 ; table s5 ) , providing support to the adaptive nature of this introgression event . collectively , these results indicate that regulatory variants affecting steady - state gene expression and transcriptional responsiveness to immune challenges , particularly those that are viral related , were preferentially introduced into european genomes via admixture with neandertals , of which some may have conferred a selective advantage to modern populations .
recent studies have offered proof of concept that eqtl mapping detects key genetic variants relevant to immunity and infection ( fairfax and knight , 2014 ) . here , using rna - seq data , we characterized , at an unprecedented level of resolution , the transcriptional response of primary monocytes to inflammatory and infectious cues .
we defined the respective contributions of natural selection and archaic admixture to differences in immune response regulation between populations . in doing so
, we identify regulatory variants and molecular phenotypes that have been important to human survival and that are of biomedical interest for the understanding of genetic susceptibility to immune - related diseases .
our analyses uncovered extensive variation , globally of moderate effect , in transcriptional responses to immune challenges between individuals of african and european descent , with the strongest differences being observed for genes with antiviral and inflammatory - related functions .
these genes are enriched in associations with cis- and trans - eqtls , and regulatory variants presenting different allele frequencies between populations account for a large fraction of the population differences in immune responses observed .
highlighting one pertinent example , we identify a reqtl ( rs2274065 ) , whose cis - action was supported by our analyses of ase , leading to tlr - mediated ncf2 downregulation in africans , where this variant is present at high frequency ( 50% ) . that this mutation has been associated with systemic lupus erythematosus ( jacob et  al . , 2007 )
suggests that lower levels of ncf2 expression may contribute to the higher prevalence and severity of this disease in africans ( fernndez et  al .
this example illustrates the value of mapping response eqtls across populations to uncover mechanisms that might explain ethnic disparities in the clinical manifestation of immune disorders .
this study also establishes that natural selection has contributed to the differences in immune responses observed between populations by providing genome - wide support that regulatory variants associated with different responsiveness to immune challenges have been targeted by positive selection . in doing so
for example , selection appears to have increased the frequency of the african - specific reqtl rs7426702 ( 39% ) , leading to stronger ccr1 downregulation following tlr1/2 activation .
interestingly , the inhibition of ccr1 limits leukocyte recruitment and prevents inflammatory responses in experimental settings ( gladue et  al . , 2006 ) .
our results thus suggest that ccr1 downregulation has conferred a selective advantage in africans , likely to favor diminished inflammation .
further support for this concept is provided by the strong selection signature detected for the european trans - eqtl at tlr1 , spanning a region shown to have evolved adaptively ( barreiro et  al . , 2009 ,
the tlr1 variant is a strong trans - regulatory hotspot associated with a gene network presenting marked population differences in the response to immune activation .
we also found that the advantageous rs5743618 allele , which impairs nf-b activity ( barreiro et  al . , 2009 ) , is associated with a global decreased expression of inflammatory response genes , consistent with an attenuated tlr1-mediated signaling beneficial to europeans .
together , our findings highlight the evolutionary tradeoff between activating efficient responses to sense microorganisms , both pathogenic and commensal , while avoiding aberrant , deleterious inflammation . genetic variation transmitted through admixture with neandertals
can also represent a source of functional , potentially advantageous variants ( vattathil and akey , 2015 ) .
relative to genome - wide expectations , we show that genetic segments introgressed from neandertals have preferentially introduced regulatory variants into european genomes , affecting steady - state expression and responses to tlr7/8 stimulation and iav .
furthermore , we report several loci presenting strong evidence of archaic ancestry that exert a regulatory effect in cis . among these
, we find the iav - induced reqtl of pnma1 , which encodes a protein that physically interacts with the iav protein pb2 and stimulates interferon production ( shapira et  al . , 2009 ) . that the pnma1 haplotype presents a frequency in europeans that is not compatible with neutral evolution , together with its strong levels of population differentiation between modern europeans and east asians , supports its contribution to european adaptation and provides a case of adaptive introgression .
the functional roles of the introgressed regulatory variants require further investigation , but our results clearly establish that archaic admixture , whether adaptive or not , has increased the diversity of the immune repertoire of contemporary europeans . collectively , our analyses provide a comprehensive view of the impact of population genetic differences on transcriptional responses to innate immunity activation and highlight evolutionarily important determinants of host immune responsiveness .
the regulatory variants identified here constitute a useful resource for evaluating the role of such variants in the molecular and cellular mechanisms underlying host immunity to infection and susceptibility to disease , both at the individual and population levels .
reagent or resourcesourceidentifierantibodiesmagnetic cd14 microbeads , humanmiltenyi bioteccat#130 - 050 - 201cd14-apc , humanmiltenyi bioteccat#130 - 091 - 243chemicals , peptides , and recombinant proteinsficoll - paquege healthcarecat#17 - 1440 - 03fetal calf serumpaa laboratoriescat#a15 - 502dimethyl sulfoxidesigma aldrichcat#d2438penicillin / streptomycinlife technologiescat#15140 - 122propidium iodidemiltenyi bioteccat#130 - 093 - 233lpsinvivogencat#tlrl-3pelpspam3csk4invivogencat#tlrl - pmsr848invivogencat#tlrl - r848 - 5critical commercial assaysnucleospin mirna kitmacherey nagelcat#740971.250rna 6000 nano kitagilent technologiescat#5067 - 1511quant - it picogreen dsdna assay kitlife technologiescat#p7589truseq rna sample prep kit v2illuminacat # rs-122 - 2001truseq sr cluster kit v3-hsilluminacat # gd-401 - 3001truseq sbs kit v3-hsilluminacat # fc-401 - 3001humanomni5-quad beadchipsilluminacat#wg-311 - 5001nextera rapid capture expanded exome kitilluminacat#fc-140 - 1006deposited datagenotyping , exome and rna sequencing dataeuropean genome - phenome archive ( ega)egas00001001895experimental models : organisms / strainshuman primary monocytesthis papern / ainfluenza a virus , strain a / ussr/90/1977(pothlichet et  al .
2012)http://cole - trapnell - lab.github.io / cufflinks / cuffdiff / combat , sva r package(johnson et  al . , 2007)https://www.bioconductor.org / wgcna(langfelder and horvath , 2008)https://www.bioconductor.org / pastaa(roider et  al . ,
2009)http://trap.molgen.mpg.de / cgi - bin / pastaa.cgitrap(thomas - chollier et  al . , 2011)http://trap.molgen.mpg.de / cgi - bin / download.cgifdrtool , r package(strimmer , 2008)http://cran.r - project.org / goseq , r package(ashburner et  al . , 2000)https://www.bioconductor.org / plink v1.9(chang et  al . ,
, 2010)http://people.virginia.edu/wc9c / king / index.htmladmixture(alexander et  al . , 2009)https://www.genetics.ucla.edu / software / admixture / bwa v.0.7.7(li and durbin , 2009)http://bio - bwa.sourceforge.net / picard tools v.1.94n
/ picardgatk v.3.2.2(depristo et  al . , 2011)https://software.broadinstitute.org / gatk / shapeit2(delaneau et  al . , 2013)http://www.shapeit.frimpute v.2(howie et  al . , 2009)http://mathgen.stats.ox.ac.uk / impute / impute_v2.1.0.htmleigenstrat(patterson et  al . , 2006)http://genepath.med.harvard.edu/reich / eigenstrat.htmmatrixeqtl , r package(shabalin , 2012)http://www.bios.unc.edu / research / genomic_software / matrix_eqtl / wasp(van de geijn et  al . , 2015)https://github.com / bmvdgeijn / waspsamtools mpileup(li , 2011)http://samtools.sourceforge.net / variant effect predictor(mclaren et  al . ,
2010)http://www.ensembl.org / info / docs / tools / vep / index.html further information and requests for reagents may be directed to , and will be fulfilled by the corresponding author lluis quintana - murci ( quintana@pasteur.fr ) .
we recruited 100 healthy , male donors of self - reported european descent ( eub ) and 100 of self - reported african descent ( afb ) , all living in belgium , at the center for vaccinology ( cevac ) of ghent university hospital ( ghent , belgium ) .
samples were collected after written informed consent had been obtained , and the study was approved by the local ethics committee ( ethics committee of the ghent university ) , the ethics board of institut pasteur ( evoimmunopop-281297 ) and the relevant french authorities ( cpp , ccitrs and cnil ) .
inclusion was restricted to donors between 19 and 50 years of age , nominally healthy at the time of sample collection .
a  case report form was obtained for all donors , including information on vital sign measurements , medication , medical history and travel .
no overrepresentation of any particular disease was observed relative to official report statistics published by the world health organization or in epidemiological studies .
serological testing was performed for all donors at the cevac , and those with serological signs of past or ongoing infection with human immunodeficiency virus ( hiv ) , hepatitis b virus ( hbv ) or hepatitis c virus ( hcv ) were excluded .
for each participant , we collected 300  ml of whole blood into anticoagulant edta - blood collection tubes and peripheral blood mononuclear cells ( pbmcs ) were isolated on ficoll - paque density gradients .
pbmcs were frozen in 90% fetal calf serum ( fcs ) and 10% dimethyl sulfoxide , at a density of 50   10 pbmcs / ml and transported in dry shipper from cevac to institut pasteur .
300   10 pbmcs were thawed , washed twice and resuspended in pre - warmed rpmi-1640 glutamax medium , supplemented with 10% fcs and penicillin / streptomycin ( complete medium ) .
monocytes were then positively selected with magnetic cd14 microbeads , according to the manufacturer s instructions .
the number of monocytes was determined with a kova glasstic slide 10 with a grid in the presence of trypan blue . for each donor ,
30   10 monocytes were split between five 25  cm non - treated flasks ( i.e. one flask per condition and five conditions per donor ) , each containing 6   10 monocytes in 9  ml of complete medium .
monocytes were allowed to rest for one hour at 37c under 5% co2 before stimulation .
purity and cell death of the isolated monocytes were assessed for all donors on a fraction of 10 cd14 monocytes stained , according to the manufacturer s instructions , with fluorescent apc - conjugated anti - cd14 antibodies and propidium iodide , respectively .
samples were then analyzed on a macsquant analyzer 10 benchtop flow cytometer ( miltenyi biotec ) and using flowjo vx.0.6 software .
the mean values obtained for all samples were 96.8% for monocyte purity and 2.1% for initial cell death rates .
monocytes were exposed to five different conditions for 6  hr , in order to capture transcriptional signatures from both an early response and the beginning of a late response , i.e. , an  intermediate response  ( huang et  al . , 2001 ) .
the choice of this time point was also based on a pilot study on the kinetics of gene expression of several key inflammatory and antiviral response genes ( il1a , il23a , il6 , il8 , tnf , irf1 and stat2 ) upon immune activation .
our results showed that 6  hr of stimulation was the best time point to capture simultaneously expression signals from early , intermediate and late response genes , with respect to other time points at 2 , 4 , 8 and 24 
one monocyte flask was left untreated as a baseline control , while the others were each exposed to one of four different immune stimuli .
these stimuli included synthetic ligands specifically activating three toll - like receptor ( tlr ) signaling pathways : 1  ng / ml ultrapure lps from e.  coli , 0.2  g / ml synthetic triacylated lipoprotein pam3csk4 , and 0.3  g / ml imidazoquinoline compound r848 .
monocytes were also infected with strain a / ussr/90/1977(h1n1 ) of the human seasonal influenza a virus ( iav ) at a moi  = 1 , and iav particles were produced as previously described ( pothlichet et  al . , 2013 ) .
after stimulation , cells were collected by centrifugation , lysed in a guanidinium thiocyanate solution provided in the nucleospin mirna kit , according to the manufacturer s instructions , and stored at 80c until rna extraction .
cellular assays were performed per batch of 30 samples from 6 individuals , including 3 africans and 3 europeans , across all 5 conditions .
total rna was extracted with the nucleospin mirna kit from macherey nagel , including the enzymatic digestion of genomic dna .
extractions were performed in batches of 30 samples ( i.e. 5 conditions for 3 africans and 3 europeans ) , and rna quality and quantity were assessed with a nanodrop spectrometer and the agilent bioanalyzer rna 6000 nano kit .
we generated a final set of 978 samples from the 200 donors fulfilling the quality and quantity criteria ( rin > 7 , quantity > 2.5  g ) for high - throughput rna - sequencing , including 200 , 188 , 197 , 193 and 200 samples for the non - simulated , lps , pam3csk4 , r848 and iav conditions , respectively .
rna was obtained from 978 of the 1000 samples , and was sequenced on an illumina hiseq2000 .
samples were then randomized before library preparation in order to obtain a balanced number of samples across ethnicity and cellular conditions per sequencing batch / lane / machine / index .
standard reagents were used for transcriptome sequencing : truseq rna sample prep kit v2 for mrna library construction , truseq sr cluster kit v3-hs for cluster generation and truseq sbs kit v3-hs for sequencing .
we pooled six samples per lane to generate outputs of around 30 million 101-bp single - end reads per sample ( ranging from 27.7 to 94.8 million reads , mean 34.4 ) ( figure  s2a ) .
genomic dna was extracted from the cd14-negative cell fraction ( i.e. non - monocyte cells ) by a standard phenol / chloroform protocol followed by ethanol precipitation .
the dna was quantified by nanodrop spectrometry and with the quant - it picogreen dsdna assay kit .
the 200 subjects studied were genotyped for a total of 4,301,332 snps on the illumina humanomni5-quad beadchips .
whole - exome sequencing was carried out for the same individuals with the nextera rapid capture expanded exome kit , on the illumina hiseq 2000 platform , with 100-bp paired - end reads .
this kit delivers 62 mb of genomic content per individual , including exons , untranslated regions ( utr ) , and micrornas .
reads were assessed for multiple quality metrics , including number of reads , nucleotide distribution and sequencing quality , and the last base of all reads was trimmed due to a fall in sequencing quality .
rna reads were then mapped onto the human grch37 genome with tophat ( kim et  al . , 2013 ) , resulting in the successful mapping of 89.9% of reads per sample on average ( minimum 67.3% ; maximum 93.7% ) .
we used the rseqc package to assess the alignment of reads with various genomic features , gc content , and gene body coverage ( figures s2b  s2e ) .
samples with uneven gene body coverage were found to be more likely to be outliers .
we used gene body coverage regularity as an indicator of library quality , removing eight samples due to irregular gene body coverage .
the remaining 970 samples were used for subsequent analyses and consisted of 200 non - stimulated ( eub : 100 , afb : 100 ) , 184 lps ( eub : 96 , afb : 88 ) , 196 pam3csk4 ( eub : 100 , afb : 96 ) , 191 r848 ( eub : 98 , afb : 93 ) , and 199 iav samples ( eub : 99 , afb : 100 ) .
cufflinks / cuffdiff ( v2.0.2 ) ( trapnell et  al . , 2012 ) was used to quantify expression levels in fpkm ( fragments per kilobase of transcript per million mapped reads ) for each annotated transcript of the genome in ensembl ( v.70 ) , and fpkm values for which cufflinks returned fail status ( < 0.5% of quantified transcripts ) were set to missing values .
gene expression data were filtered to remove genes with low levels of expression ( mean fpkm  < 1 in all conditions ) and their quality was checked by principal component analysis ( pca ) .
pca captured differences between conditions and populations on the first two axes , but we tested for additional causes of technical variability , by fitting , for each gene , a mixed model of gene expression as a function of condition , population , and technical covariates , including total rna concentration , rin , percentage of high - quality bases ( q30 ) , mean gc content , library concentration , 5/3 coverage bias ( measured as the mean difference in coverage between the 5 and 3 ends of the gene ) as continuous covariates , and date of experiment , library preparation batch , sequencing batch , sequencer used , sequencing index , and sequencing lane as putative batch effects .
putative batch effects were modeled as random effects to prevent the loss of degrees of freedom , whereas all other covariates ( condition , population and continuous covariates ) were included as fixed effects , giving the following model for gene i and sample j : log2(1+fpkmij)=+cond+pop+cccovariatec+kk , jbatch+ijwhere  is the intercept , cond and pop are the fixed effect of the condition and population on sample j , c are the fixed effect of continuous covariates on sample j , the k , jbatch are the random effects of batch covariate k , on sample j , and ij are the residuals . the proportions of genes affected by each factor are reported in figure  s2f for various levels of explained variance .
we observed that gc content , 5/3 bias , date of the experiment and library batch were among the strongest confounding factors , and accordingly corrected the data for these factors before analysis , following the pipeline detailed in figure  s3 .
first , we adjusted the data for gc content and 5/3 bias using linear models .
then , we imputed missing values by k - nearest neighbor imputation , and adjusted for experiment date and library batch by sequentially running combat ( johnson et  al . , 2007 ) for each batch effect , with condition and population as covariates .
batch - corrected gene expression levels , in fpkm , were then recalculated from the adjusted transcript level estimates .
refitting our linear mixed model confirmed that correction was satisfactory for most of the technical covariates ( figure  s2f ) .
the reproducibility of our rna - seq experiments was assessed by performing technical and biological replicates on seven independent donors ( 4 afb and 3 eub ) across the five experimental conditions .
we showed that ( i ) the coefficients of variation of genes within technical replicates were consistently , and significantly , smaller in magnitude and less variable than those within biological replicates ( wilcoxon rank - sum test , p  < 10 ; figure  s2 g ) , and ( ii ) technical replicates exhibit higher correlation coefficients ( r ) between samples with respect to the distribution of r values calculated from pairwise comparisons between biological replicates ( figure  s2h ) .
modules of genes presenting correlated expression patterns , extracted from log - transformed fpkm data , were defined by weighted correlation network analysis ( wgcna ) ( langfelder and horvath , 2008 ) . in our setting of immune response activation , this analysis detects modules of correlated genes that can reflect either shared regulation by common transcription factors , or regulation by independent transcription factors with similar patterns of activation upon stimulation .
tukey s biweight correlation was used as a measure of gene relatedness to reduce the influence of outliers , and correlations were measured across all 970 samples .
the scale - free topology of the networks was assessed for various values of the  shrinkage parameter , according to wgcna user manual , and the default value of   = 6 appeared to give a satisfactory fit to scale - free topology .
signed clustering of genes ( grouping only positively correlated genes ) was used to simplify the interpretation of the extracted modules .
we also found that varying the level of shrinkage (   = 5 or 6 ) or the depth of the clustering ( deepsplit parameter set to 3 or 4 ) had only a mild impact on the number of clusters or the enrichments obtained , confirming the robustness of these analyses .
for each module , we used pastaa ( roider et  al . , 2009 ) to identify transcription factor binding site motifs overrepresented within the annotated proximal promoters of the genes within each module .
we first defined the proximal promoter region for each gene as the region extending 200  bp on either side of the transcription start site ( tss ) of the most abundant transcript on the basis of cufflinks fpkm estimates .
we then used the transcription factor affinity prediction ( trap ) method ( thomas - chollier et  al . , 2011 ) to measure the binding affinities of each transcription factor present in the jaspar core vertebrate database ( mathelier et  al . , 2014 ) with the proximal promoters of the 12,578 expressed genes , and these affinities were then used as the input for pastaa enrichment analysis .
we reported only enrichments significant at a false discovery rate ( fdr ) of 0.05 with a fold - change ( i.e. observed / expected ) greater than 1.2 . for each module
, we represent the transcription factor binding sites with the highest value for the lower limit of the odds ratio confidence interval .
differential expression was assessed directly from log - transformed fpkm , using t tests for each condition .
fdr was then calculated jointly for all conditions , with the r package fdrtool , and genes differentially expressed between populations ( popdegs ) were defined as genes presenting an absolute log2 fold change between populations  |log2fcpop| 
we then calculated the fold - change in expression after stimulation relative to the basal state , and used t tests to determine whether there was a differential response .
population differential response genes ( popdrgs ) were then defined as popdegs for which there was a differential response between populations under stimulated conditions , at fdr  < 0.05 ( i.e. the transcriptional response to treatment , relative to the basal state , differed between populations ) , resulting in a larger difference in expression after stimulation.|log2fcpopstimulated|>|log2fcpopbasal| all gene ontology ( go ) enrichment analyses were performed with goseq package ( ashburner et  al . , 2000 ) , using the default settings , with the 12,578 expressed genes as the background set . only enrichments significant at fdr of 0.05 and with a fold - change ( i.e. observed / expected ) greater than 1.2 are reported .
using plink v1.9 ( chang et  al . , 2015 ) , we removed snps that : ( i ) were typed with probes mapping to several genomic locations ( n  =  12,440 ) , ( ii ) presented a poor genotype clustering ( gentrain score  < 0.35 ; n  = 809 ) , ( iii ) had the same chromosomal position as another snp in dbsnp b138 ( n  = 6,968 ) , ( iv ) were not reported in dbsnp b138 ( n  = 5,311 ) , ( v ) presented a call rate  < 95% ( n  = 79,310 ) , ( vi ) were monomorphic in our sample ( n  = 652,385 ) , ( vii ) were located on the sex chromosomes ( n  = 50,994 ) , and ( viii ) presented a hardy - weinberg p  < 10 in afb or eub populations ( n  = 4,007 ) .
the snp call rate for the 200 individuals was 99.8% on average , ranging from 98.89% to 99.99% .
no evidence was found for 2-degree cryptic relatedness ( kinship coefficient > 0.07 ) in king ( manichaikul et  al . , 2010 ) , or for sex mismatch , for any of the individuals .
two afb individuals presented an excess of heterozygosity ( <  3sd from the population average ) , as a result of their moderate levels of non - african ancestry , as estimated using admixture .
read - pairs were first mapped onto the human grch37 genome with bwa v.0.7.7 ( li and durbin , 2009 ) , and reads duplicating the start position of another read were marked as duplicates with picard tools v.1.94
we used gatk v.3.2.2 ( depristo et  al . , 2011 ) for base quality score recalibration (  baserecalibrator  ) , insertion / deletion realignment (  indelrealigner  ) , and snp and insertion / deletion discovery for each sample (  haplotype caller  ) .
individual variant files were combined with  genotypegvcfs  and filtered with  variantqualityscorerecalibration . 
individual coverage was 52.32   on average , ranging from 33.84 to 100.59   , and individual breadth of coverage at 5   was 92.42% , ranging from 83.5% to 95.0% .
we removed those of the 540,990 snps obtained that : ( i )  were triallelic ( n  = 11,925 ) , ( ii ) presented a call rate  < 95% ( n  = 44,716 ) , ( iii ) were located on the sex chromosomes ( n  = 8,369 ) , and ( iv ) presented a hardy - weinberg p  < 10 in afb or eub populations ( n  = 4,510 ) .
the application of these quality - control filters resulted in the retention of 471,740 snps . before merging the omni5 and exome datasets , we first checked genotype concordance for 169,406 snps common to the two platforms .
we flipped alleles for 8,025 snps with incompatible allelic states , and removed 119 snps with alleles that remained incompatible after allele flipping .
the concordance rates for each of the 200 individuals exceeded 99% , confirming an absence of errors during dna sample processing .
of the 8,155 snps with discordance rates > 1% , 296 were due to c / g or a / t snps , and high genotype concordance between the two dna typing technologies was restored by allele flipping .
the entire omni5 and exome datasets ( 3,489,108 and 471,740 snps , respectively ) were then merged , yielding a final concordance rate of 99.93% , for a total of 3,782,260 snps . before imputation
, we phased the data with shapeit2 ( delaneau et  al . , 2013 ) , using 500 conditioning haplotypes , 50 mcmc iterations , 10 burn - in and 10 pruning iterations .
snps and allelic states were then aligned with the 1,000 genomes project imputation reference panel ( phase 1 v3.2010/11/23 ) .
we removed 8,705 snps with identical positions in our data and in the reference panel but incompatible alleles , even after allele flipping , and 4,137 snps with c / g or a / t alleles .
genotype imputation was performed with impute v.2 ( howie et  al . , 2009 ) , considering 1-mb windows and a buffer region of 1 mb . of the 38,098,530 snps obtained after imputation , we removed snps that : ( i ) presented an information metric below 0.8 ( n  =  18,085,215 ) , ( ii ) had a duplicate ( n  = 59,914 ) , ( iii ) presented a call rate  < 90% ( n  = 329,910 ) , and ( iv ) were monomorphic ( n  =  4,053 ) .
, we estimated correlation coefficients r between true genotypes ( i.e. , obtained by illumina array genotyping or exome sequencing ) and imputed genotypes for the same snps ( i.e. , obtained by artificially removing genotyped snps from the data before imputation and then imputing them ) . in very good agreement with recent studies ( auton et  al . ,
2015 ) , the average correlation coefficient was 95.6% across all genotyped snps with information metric > 0.8 ( 93.6% for snps with maf  < 0.10 and 97.7% for snps with maf > 0.10 ) .
this shows that our stringent quality filters ensure that only accurately imputed snps are analyzed .
two methods were used to infer the genetic structure of our population set of 100 african - descent and 100 european - descent belgians ( afb and eub , respectively ) .
because both methods assume linkage equilibrium among snps , we pruned the datasets for snps in linkage disequilibrium ( ld ) , using plink v1.9 ( chang et  al . , 2015 ) . specifically , we removed snps in 50-snp windows that present ld r > 0.5 (  indep - pairwise 1000 10 0.5  option ) .
the first model - based approach , admixture ( alexander et  al . , 2009 ) , estimates the proportions of each individual s genome originating from k ancestral populations , k being specified a priori .
this analysis was performed on 229,320 independent snps and 789 individuals from 22 populations , including eub and afb , together with a selection of representative populations from sub - saharan africa , north africa , the near east and europe ( altshuler et  al . , 2010 , behar et  al . , 2010 , patin et  al . , 2014 ) .
we made k vary from 2 to 10 . to obtain the most supported results and test for their stability
, all admixture analyses were run five times with different random seeds , for each k value .
we kept results providing the lowest cross - validation error ( cv ) among iterations .
the second model - free approach is the principal component analysis ( pca ) implemented in eigenstrat ( patterson et  al . , 2006 ) .
we used this approach to describe the local genetic sub - structure of afb and eub separately .
the analysis for afb was performed on 341,593 independent snps and 511 individuals from 7 western and central african populations , while the analysis for eub was performed on 182,572 independent snps and 220 individuals from 13 european populations ( altshuler et  al . , 2010 , behar et  al . , 2010 ,
, only variants with a minor allele frequency ( maf )  0.05 in the population studied were retained in the analysis , resulting in a set of 10,278,745 snps ( i.e. , corresponding to the merged genotyping and imputed snp dataset ; 8,913,090 snps in africans and 6,178,808 snps in europeans , ) .
pc1 and pc2 of the genotype matrix were included in the model to account for possible population stratification .
the inclusion of additional pc in the model ( up to pc6 ) was tested and showed highly consistent results ( i.e. , correlation of -log10 p - values of eqtl > 0.95 ) .
for each gene , snps were considered  local ,  likely cis - acting , if they were located less than 1 mb away from the gene transcription start or end site .
eqtl mapping was performed separately for each population and condition , and false - positives due to outliers were prevented by discarding , from the analysis , eqtl associations that did not pass a p - value threshold of 10 for local eqtls , and 10 for trans - eqtls in kruskal - wallis rank tests . for both cis- and trans - eqtls ,
fdr was computed by mapping eqtls on 100 datasets with genotypes permuted within each population .
we then kept , after each permutation , the most significant p - value per gene , across all conditions and populations .
finally , we computed the false discovery rate associated with each p - value threshold in cis or in trans , and subsequently selected the p - value  threshold that provided a 5% fdr , leading to p  = 7.67   10 and p  = 2.7   10 for cis- and trans - eqtls , respectively . for local eqtls , we report only the snp at which the strongest association was observed ( i.e. , eqtl peak - snp ) .
when multiple snps in perfect ld fell within the peak , only one snp is reported .
eqtls for which the eqtl peak - snp had an allelic effect size ( |eqtl| ) below 0.2 were discarded from further analysis .
we next mapped fold - changes between the basal and stimulated states using matrixeqtl , and defined response eqtls ( reqtls ) as stimulated eqtls associated to a significant difference in response to stimulation ( p  < 10 , |eqtlstim|>|eqtlbasal| ) . for trans - eqtls , we reported , within each 1 mb window , the snp for which we observed both the largest number of trans - associated genes and strongest p - value of association . furthermore , for each trans - eqtl that passed genome - wide significance at p  = 2.7   10 ( fdr of 5% ) , we performed a snp - based analysis to identify genes regulated in trans by the eqtl at a bonferroni p  < 0.05 , correcting for the 12,578 genes tested within the condition where the eqtl was found . to estimate the fraction of population differences in gene expression that can be attributed to genetic variants , we used a two - step strategy .
first , we consider the set of all snps in ld ( r > 0.5 ) with the eqtl peak - snp , in the population where the eqtl was discovered and fine map the eqtl signal by fitting across populations the following linear model : expressionj=+.snpj+.popj+jwhere snpj is the genotype of the individual j for the variant under study , popj is a binary variable indicating the population origin ( 0 for europeans and 1 for africans ) , and j is a random , normally distributed residual . in this model
,  is the intercept ,  reflects the effect of the derived allele of the snp on gene expression , and  estimates the fold change in expression between populations observed for individuals with identical genotype ( i.e. gene expression differences that are not explained by genetics ) .
we next focused on the snp showing the strongest association p - value with gene expression across populations , and estimated the difference in population expression that is attributable to the snp as : fcsnp = fcpopwith  representing  set to ensure that the ratio of fcsnp / fcpop is between 0 and 1 , i.e.   = 0 , if the sign of  differs from that of fcpop ;   = fcpop , if || > |fcpop| ; and   =  otherwise .
the percentage of population differences in expression that is attributable to genetics is then given by the ratio fcsnp / fcpop .
we aimed at distinguishing population specific eqtls ( i.e. , snps present at similar frequencies in both populations but having an effect on gene expression in one population only ) from eqtls detected in one population only due to population differences in allelic frequencies .
to do so , we first focused on the 1,109 genes associated with an eqtl ( including 363 genes associated with a reqtl ) where all snps in ld ( r > 0.5 ) with the eqtl peak - snp were present at frequency > 5% in both populations .
we then tested these eqtls for replication at a relaxed threshold of p  < 0.05 across all snps at the locus , to decrease the false negative rate , and focused on the 127 genes for which the eqtl was not replicated ( including 28 genes with a reqtl ) . finally , we considered as population - specific , eqtls whose effect size was significantly different between populations . to do so ,
we fit , for each snp at the locus ( r > 0.5 with the eqtl peak - snp ) , the following linear model : expressionj=+.snpj+.popj+.snpjpopj+jwhere snpj is the genotype of the individual j for the variant under study , popj is a binary variable indicating the population origin ( 0 for europeans and 1 for africans ) , and j is a random , normally distributed residual . in this model ,  reflects the effect of the derived allele of the snp on gene expression ,  estimates the fold change in expression between populations observed for individuals with identical genotype , and  captures the differences in eqtl effect size between populations .
such a model allows to test for a difference in eqtl effect size between populations by testing the null hypothesis ,   = 0 ( interaction test ) . to be conservative and to account for the uncertainty in detecting the causal variant at the eqtl
, we considered an eqtl as population specific if all snps in ld ( r > 0.5 ) with the eqtl peak - snp presented a significant interaction p-value.pinteraction(locus)=maxsnplocuspinteraction(snp ) .
we then considered eqtls ( or reqtls ) as being population specific when the interaction p - value at the locus was lower than 10 ( corresponding to fdr  < 0.01 ) , leading to a final set of 16 population - specific eqtls ( including 5 reqtls ) .
regulatory features were extracted from ensembl regulatory build v80 ( zerbino et  al . , 2015 ) , which contains regulatory element predictions based on open chromatin regions and histone marks from encode and the roadmap epigenomics datasets ( ernst and kellis , 2015 , kundaje et  al . , 2015 ) .
snps overlapping a regulatory element were then classified into four categories : promoter , promoter flanking , enhancer , and ctcf binding sites .
similarly , encode uniformly processed transcription factor binding site ( tfbs ) clusters ( v3 ) ( ernst and kellis , 2015 ) were downloaded from ucsc , and their overlap with the physical position of all snps was determined .
we then used fisher s exact test to assess the eqtl enrichment of specific tfbs or regulatory elements , considering the peak - snp at each locus , or a randomly selected snp if multiple snps in perfect ld were found .
all snps with a maf  0.05 located less than 1 mb away from an expressed gene were used to constitute the background set . in each condition
( or combination of conditions ) , only the tfbs with the highest values for the lower limit of the odds ratio confidence intervals are reported . for the quantification of allele - specific imbalance , we focused on exonic snps genotyped as heterozygous in our exome data , excluding snps with discordant genotypes in the omni5 data .
we used bwa mem ( v.0.7.7 ) ( li and durbin , 2009 ) to remap rna - seq reads onto the hg19 genome for all 970 samples , and extracted all reads aligned with a genetic variant .
we reduced mapping bias , by using wasp ( van de geijn et  al . , 2015 ) to exclude reads overlapping with known variants ( based on dbsnp138 ) likely to alter the read mapping location .
briefly , for each read overlapping one or more dbsnp variants , wasp creates alternative reads consisting of all possible combinations of reads given these snps .
it then remaps the alternative reads to the genome , and keeps the original read only if all alternative versions of the read map to the same position .
finally , samtools mpileup ( li , 2011 ) , with option -d 10000 , was used to count the number of reads mapping to each allele at heterozygous loci . the allelic ratio ( ar )
was defined for each site as the proportion of minor alleles among all reads , and the allelic imbalance ( ai ) was defined as the absolute deviation from a balanced ratio of 0.5 ( i.e. ai  = |ar-0.5| ) .
we mapped allele - specific expression qtls ( aseqtls ) , by estimating the allelic ratio on the subset of eqtl - genes with sufficient expression coverage at heterozygous exonic snps ( n  10 reads ) in at least five individuals of each eqtl genotype ( heterozygous / homozygous ) .
we extracted the phase information between the strongest local eqtl and the exonic snp , and tested the correlation between the ar and the phased eqtl genotypes ( coded 0 for homozygous , and  1 for heterozygotes with variants in phase or in the opposite phase ) , in a gene- , condition- and population - specific manner .
similarly , allele - specific response qtls ( asrqtls ) were mapped by assessing the correlation between the phased reqtl genotypes and the change in ar at the exonic site after stimulation .
the power to detect aseqtls was computed for various eqtl effect sizes || , number of observations
n and number of reads per exonic snp n. we assumed the same number of observations for heterozygous and homozygous genotypes at the eqtl , and equal coverage across all exonic snps .
power was then computed for a standard t - test assuming a mean allelic ratio nalternative / nreference of 0.5 in homozygous individuals and 2/(1 + 2 ) in heterozygous individuals .
residual variance was set to 0.25/n to match that of a binomial distribution with parameters ( 0.5 , n ) . to ensure sufficient power when exploring ase within single individuals , we considered a higher coverage of heterozygous exonic snps ( n  30 reads ) , and used a binomial test to evaluate allelic imbalance .
we also excluded sites at which one allele accounted for  less than 2% of the reads or less than 3 reads in total , as such sites might be subject to genotyping errors or systematic mapping  biases .
the fdr was first calculated across all snps , individuals and conditions , using fdrtool , and ase was defined as the combination of significance at fdr  = 0.05 and an absolute log2 fold change of expression between alleles of more than 0.2 ( |log2(nalt / nref)| > 0.2 ) . for each significant ase event in stimulated conditions
, we checked for differences in allelic imbalance relative to the non - stimulated condition , and defined allele - specific response as the subset of ase displaying significantly higher allelic imbalance ( p  < 10 , fisher s exact test ) after stimulation with respect to the basal state .
finally , we used simulations to evaluate fdr among the set of genes with at least one ase / asr event .
we generated 1,000 null datasets , by randomly reassigning reads to the alternative and reference allele with equal probability , and estimated the number of genes with at least one significant ase or asr event at each p - value threshold .
we then computed fdr as the ratio of the average number of genes with ase in our resampling to the observed number of genes with ase at the same p - value threshold . for each exonic snp for which we quantified ase , we used variant effect predictor ( vep ) ( mclaren et  al . , 2010 ) with ensembl v.70 transcript annotation to identify the set of transcripts overlapping the variant , and cufflinks fpkm to identify the most strongly expressed overlapping transcript in the individual / condition concerned
. vep annotation was then used to classify variants , according to the most abundant transcript , as synonymous ( synonymous_variant/ non_coding exon_variant ) , missense ( missense_variant ) or nonsense ( stop_gained , stop_lost , initiator_codon_variant ) .
enrichment in rare coding variants was then assessed using fisher s exact test comparing each category with synonymous variants .
we used two metrics , fst and ihs , which detect signals of population - specific positive selection , i.e. , mutations that provided a selective advantage to a specific human population .
fst measures population differentiation by comparing the variance of allele frequencies within and between populations ( holsinger and weir , 2009 ) , as local positive selection tends to increase allele frequency differences between populations .
as fst is a population pairwise comparison , we derived a directional fst , equal in absolute value to the pairwise fst but with a positive sign if the derived allele was more frequent in the population studied , and a negative sign otherwise .
this enables to distinguish selection events that likely occurred in africans from those that likely occurred in europeans .
the integrated haplotype score ( ihs ) measures the degree of extended haplotype homozygosity of the putatively selected allele over that of the putatively neutral allele ( voight et  al . , 2006 ) , as the long - range associations of the selected mutation with neighboring snps are not disrupted by recombination .
furthermore , we used a composite selection score ( css ) allowing to capture signals of recent , strong selective events , by combining fst and ihs .
the css was designed to identify variants with both a higher derived allele frequency in one population ( positive value of directional fst ) , and a longer haplotype length around the derived allele of the variant in that population ( characterized by a negative ihs value ) .
it was computed for all variants with derived allele frequency 0.2   daf  0.95 from genome - wide ranks of both directional fst 
 ( rfst ) and ihs ( rihs ) , attributing the highest rank to positive values of fst and negative values of ihs , respectively .
we defined the css as following : css = rank(rfst.rihs)nobswith nobs being the total number of variants with 0.2   daf  0.95 in the population studied .
css ranges from 0 to 1 , and increases with the strength of positive selection targeting the derived allele .
finally , we used the cross - population composite likelihood ratio score , xp - clr , a region - based metric detecting extended regions where the allele frequencies of multiple contiguous markers are distorted from the prediction under neutrality ( chen et  al . , 2010 ) .
xp - clr detects classical selective sweeps as well as selection events on pre - existing alleles ( standing variation ) .
xp - clr was scored every 2,000  bp , using windows of 0.2 cm and downsampling to 200 snps per window . to map selection signals at haplotypes containing eqtls , we determined , for each statistic ( ihs , fst , or css ) , the strongest signal of selection on derived alleles of all snps in high ld ( r > 0.8 ) with the eqtl peak - snp . to assess significance
, we then compared , for each population and condition , the mean of these values across all eqtls / reqtls , with the expected distribution obtained from resampling 10,000 sets of random snps matched for maf ( using bins of maf of 0.05 ) and the number of snps in ld ( r > 0.8 , using bins of 0 - 2 , 3 - 5 , 6 - 10 , 11 - 20 , 21 - 50 , and > 50 snps in ld ) .
similarly , for xp - clr , we compared the mean of xp - clr scores at eqtls / reqtls ( considering the region that contains the eqtl peak - snp ) , to the expected distribution obtained from resampling 10,000 sets of random snps matched for maf and ld patterns . to identify candidate eqtls under selection , we used an outlier approach where we computed the top 1% values of fst and ihs at the genome - wide level , focusing on signals consistent with selection on derived alleles , within each population separately . to support the adaptive nature of candidate eqtls ,
we computed neutral p values for each statistic using simulations based on validated demographic models of africans and europeans ( grossman et  al . , 2013 ) .
furthermore , we tested for local enrichment of outliers ( top  1% signals ) within a 100kb - window around each eqtl ( 50 kb on each side ) , similarly to previous work ( grossman et  al . , 2013 ,
the proportion of outliers of fst or ihs ( 1% threshold ) was computed from snps with daf  0.2 in a 100 kb window around each putatively selected locus .
significance was assessed from a beta binomial distribution fitted , in each population separately , to the observed genome - wide distribution of the proportion of outliers , to account for variations in the number of snps at each locus .
we determined the level of neandertal ancestry of the detected eqtls , by first defining an  archaic eqtl  as an eqtl for which regulatory variants were introduced into european genomes by introgression from archaic hominins .
we identified such eqtls using the complete genome sequence of neandertal from altai ( prfer et  al . , 2014 ) .
briefly , the 1000 genomes phase 3 variants ( auton et  al . , 2015 ) were considered as of putative archaic origin ( archaic snps , or asnps ) if the neandertal allele was present in at least one non - african individual and absent from the yoruba population . according to this definition
we rendered the analysis more conservative , by further restricting the definition of asnps to those in regions of the modern human genome for which neandertal ancestry has been predicted with a high degree of confidence ( marginal probability of neandertal ancestry  0.9 and a genetic length  0.02 cm ) ( sankararaman et  al . , 2014 ) .
this resulted in a final set of 197,959 asnps , of which 77,823 presented a maf > 0.05 .
more than 96% of these asnps had an archaic allele frequency below 1% in our african samples ( who are slightly differentiated from the yoruba of 1000 genomes phase 3 ) , consistent with a strong enrichment in true archaic variants . to account for ld between asnps and characterize haplotypes that were inherited from neandertal , we used plink ( chang et  al . , 2015 ) to extract a set of 924,362 genome - wide snps tagging all european variants at an r > 0.8 . among these
, 9,677 tagged all asnps in europeans and are referred to here as  archaic tagging snps . 
they were not necessarily asnps themselves , reflecting the fact that haplotypes inherited from neandertals can harbor a mixture of different variants ( i.e. variants that appeared in the neandertal lineage , and ancient variants pre - existing in both lineages before admixture , but for which one allele is carried almost exclusively by neandertal haplotypes in modern europeans ) .
we explored the effect of introgression from neandertals on the immune repertoire of europeans , by counting , for each condition , the number of eqtls overlapped by at least one archaic tagging snp ( or for which the archaic tagging snp overlapped the reqtl in  stimulated conditions ) , referred to here as archaic eqtls .
we then compared the number of archaic eqtls detected with the number of snps expected to overlap the eqtl , when resampling snps tagging non - archaic haplotypes , at random from genic regions ( < 1 mb from a gene ) .
we resampled 1,000 sets of 9,677 tagsnps with the same allele frequency spectrum as the 9,677 archaic tagging snps found in europeans , and determined their overlap with ( r)eqtls , to assess the significance of our observations .
we  finally report only the archaic ( r)eqtls for which at least 2 asnps were found to be in high ld with the eqtl peak - snp , and ( ii ) the haplotype containing the largest number of archaic alleles within the eqtl was sufficiently long for the formal exclusion of incomplete lineage sorting .
the presence of asnps in present - day humans can be explained either by introgression or by incomplete lineage sorting ( ils ) .
ils occurs when an ancestral variant predating the split between humans and neandertals is retained in both lineages , but lost from a specific human population ( i.e. the african population ; figure  s7a ) .
given the time since introgression ( 47,000 - 65,000 years ago ) , the haplotypes containing alleles resulting from ils would be expected to be shorter than those containing an asnp introgressed from neandertals .
we distinguished between these two scenarios by first defining the core archaic haplotype for each eqtl as the haplotype within the eqtl carrying the longest stretch of archaic alleles , and then determining whether its size exceeded the expected length of haplotypes assuming an ils model .
we used the approach described by ( huerta - snchez et  al . , 2014 ) and the most conservative parameters for the age of altai neandertal and denisovan bones reported by ( dannemann et  al . , 2016 ) .
we used the mean recombination rate calculated for a region composed of the core archaic haplotype in a region of 1 mb surrounding the haplotypes ( 500 kb on either side of the eqtl ) in the 1000 genomes ceu individuals ( phase 1 ) .
p values were adjusted for multiple testing with the benjamini - hochberg procedure . to test if archaic eqtls result from adaptive introgression , we used both empirical and simulation - based approaches .
the empirical approach tests if archaic alleles at eqtls are more frequent in europeans than expected , by comparing their levels of genetic differentiation between european and east asian populations with respect to the genome - wide distribution of asnps , similarly to a recent study ( vernot and akey , 2014 ) .
the rationale is that an archaic allele that introgressed into europeans and east asians 40,000 - 50,000 years ago and was advantageous in one population only should present today unusually high levels of genetic differentiation , relative to genome - wide expectations .
we thus computed the genome - wide distribution of fst between european and east asian populations at archaic snps , using the 1000 genomes project phase 3 ( auton et  al . , 2015 ) , and estimated the empirical p value for candidate archaic eqtl snps by dividing their rank by the total number of archaic snps .
we next tested if the high frequency of archaic alleles at candidate ( r)eqtls is compatible with a neutral model of evolution , using simulations .
importantly , while a detailed demographic model of neandertals is not specifically required for such simulations , we need an estimated site frequency spectrum ( sfs ) of neandertal alleles in europeans at the time of their introgression , which is unknown .
we thus used three different approximated sfs that rely on different assumptions , detailed below . in each case , the simulation - based p value for candidate archaic eqtl snps was obtained by comparing observed frequencies to the neutral simulated sfs in the current generation .
although each of the simulations used present different limitations related to the sfs of archaic alleles at the time of introgression , they can collectively provide information about the expected fate of introgressed alleles under simple scenarios .
the first approximated sfs of neandertal alleles in present - day european populations , called here  sankararaman sfs ,  was retrieved from a previous study ( sankararaman et  al . ,
this sfs was obtained assuming that archaic alleles in europeans 40,000 ya ( i ) evolve under neutrality , and ( ii ) could not have a frequency larger than 4% . for the second approximated sfs , called here  fixed - in - neandertal sfs ,
 we relaxed the assumption that archaic ancestral frequencies were lower than 4% and used 100,000 forward simulations based on the wright - fisher model to simulate frequency changes of archaic alleles in europeans since neandertal introgression , 1,440 - 2,200 generations ago ( vernot and akey , 2015 ) .
we used the best - fit demographic model of tennessen and colleagues ( tennessen et  al . , 2012 ) to model changes in the effective population size of europeans
( i.e. , from an ancestral ne of 1,032 , two successive exponential growths with rate 0.31% and 1.95% occur 920 and 205 generations ago , respectively ) .
we assumed here that ( i ) archaic alleles evolve under neutrality , and ( ii ) neandertal alleles that segregate today in europeans were most likely fixed in neandertals at the time they were introgressed .
this second assumption is conservative , as we neglect all rare neandertal alleles that had a higher probability to be lost by genetic drift .
we modeled the sfs of archaic alleles in europeans 40,000 ya by a gaussian distribution with average 5% ( i.e. , the estimated neandertal ancestry in european ancient dnas from this period ; ( fu et  al . , 2016 ) ) and 1% standard deviation .
for the third approximated sfs , called here  ancient dna - based sfs ,  we sought to circumvent the uncertainty inherent to the estimation of past neandertal allele frequencies , by retrieving them from maximum likelihood estimates in ancient dnas of european hunter - gatherers , early farmers and steppe herders ( mathieson et  al . , 2015 ) .
we computed the sfs of neandertal alleles in european populations 8,000 ya based on the 5,900 snps that were detected in eub as asnps and that were covered in the mathieson s study , to approximate the sfs of archaic alleles 320 generations ago . for convenience
, we fitted to this observed sfs a beta distribution (   = 1.21 ,   = 10.23 ) .
we then used 100,000 forward simulations under the same wright - fisher model with two exponential growths , to simulate the fate of neutral alleles during the last 320 generations .
this simulation analysis only tests if archaic alleles at candidate eqtls have been under positive selection in the last 8,000 years .
genome - wide snp genotyping , whole exome sequencing and rna - sequencing data generated in this study have been deposited in the european genome - phenome archive ( ega ) under accession code ega : egas00001001895 .
, m.r . , and l.q .- m . analyzed and interpreted results and wrote the paper with input from all authors
 summaryhumans differ in the outcome that follows exposure  to life - threatening pathogens , yet the extent of population differences in immune responses and their genetic and evolutionary determinants remain undefined . here , we characterized , using rna sequencing , the transcriptional response of primary monocytes from africans and europeans to bacterial and viral stimuli  ligands activating toll - like receptor pathways ( tlr1/2 , tlr4 , and tlr7/8 ) and influenza virus  and mapped expression quantitative trait loci ( eqtls ) . 
 we identify numerous cis - eqtls that contribute to the marked differences in immune responses detected within and between populations and a strong trans - eqtl hotspot at tlr1 that decreases expression of pro - inflammatory genes in europeans only . 
 we find that immune - responsive regulatory variants are enriched in population - specific signals of natural selection and show that admixture with neandertals introduced regulatory variants into european genomes , affecting preferentially responses to viral challenges . 
 together , our study uncovers evolutionarily important determinants of differences in host immune responsiveness between human populations . 
Introduction
Results
Discussion
STARMethods
Author Contributions