Signatures of Archaic Adaptive Introgression in Present-Day Human Populations Article Open Access

22 Pages·2017·2.33 MB·English

Signatures of Archaic Adaptive Introgression in Present-Day
Human Populations
Fernando Racimo,‡,1 Davide Marnetto,2 and Emilia Huerta-S(cid:2)anchez*,3
1DepartmentofIntegrativeBiology,UniversityofCaliforniaBerkeley,Berkeley,CA
2DepartmentofMolecularBiotechnologyandHealthSciences,UniversityofTorino,Turin,Italy
3SchoolofNaturalSciences,UniversityofCaliforniaMerced,Merced,CA
‡Presentaddress:NewYorkGenomeCenter,NewYork,NY
*Correspondingauthor:E-mail:[email protected]
Associateeditor:JohnNovembre
Abstract
D
ComparisonsofDNAfromarchaicandmodernhumansshowthatthesegroupsinterbred,andinsomecasesreceivedan ow
n
evolutionary advantage from doing so. This process—adaptive introgression—may lead to a faster rate of adaptation lo
a
thanispredictedfrommodelswithmutationandselectionalone.Withinthelastcoupleofyears,aseriesofstudieshave d
e
d
identified regions of the genome that are likely examples of adaptive introgression. In many cases, once a region was fro
ascertained as being introgressed, commonly used statistics based on both haplotype as well as allele frequency infor- m
mationwereemployedtotestforpositiveselection.Introgressionbyitself,however,changesboththehaplotypestruc- http
tTuhreereafonrde,tphaettdeirsntrsibguentieornateodf balyleilnetrforgeqreusesniocnieas,lotnheums acoynlefoaudntdoinfaglsteraindfietiroennaclesteosftpsofsoirtivdeesteelcetcintigonp.oHsietrievewseeelexcptlioorne. s://ac
a
modelsinvolvingbothintrogressionandpositiveselectiontoinvestigatethebehaviorofvariousstatisticsunderadaptive de
m
introgression. In particular, we find that the number and allelic frequencies of sites that are uniquely shared between ic
archaic humans and specific present-day populations are particularly useful for detecting adaptive introgression. We .ou
p
then examine the 1000 Genomes dataset to characterize the landscape of uniquely shared archaic alleles in human .c
o
populations.Finally,weidentifyregionsthatwerelikelysubjecttoadaptiveintrogressionanddiscusssomeofthemost m
/m
promisingcandidategeneslocatedintheseregions. b
e
/a
Keywords:neanderthal,denisova,adaptiveintrogression,ancientDNA. rtic
le
-a
b
s
Introduction 2011;Xuetal.2011;Jeongetal.2014;Hackingeretal.2016). trac
Theselected Tibetan haplotype is likely tohavebeen intro- t/3
There is now a large body of evidence supporting the idea ducedinthehumangenepoolbyDenisovansorapopulation 4/2
that certain modern human populations admixed with ar- /2
A closelyrelatedtothem(Huerta-S(cid:2)anchezetal.2014;Huerta- 96
r chaic groups of humans after expanding out of Africa. In SanchezandCasey2015). /26
t particular, non-African populationshave1–2% Neanderthal Inthisstudy,wefirstusesimulationstoassessthepowerto 33
i 3
c ancestry (Green et al. 2010; Pru¨fer et al. 2014), and detectAIusingdifferentexploratorysummarystatisticsthat 71
le MelanesiansandEastAsianshave3%and0.2%ancestry,re- donotrequiretheintrogressedfragmentstobeidentifieda by
spectively, from Denisovans (Reich et al. 2010; Meyer et al. g
priori.Someoftheseareinspiredbythesignaturesobserved u
e
2012;Pru¨feretal.2014). s
in EPAS1, which contains an elevated number of sites with t o
Recently,ithasbecomepossibletoidentifythefragments n
allelesuniquelysharedbetweentheDenisovangenomeand 0
ofthehumangenomethatwereintrogressedandsurvivein 8
Tibetans.Wethenapplythesestatisticstorealhumangeno- A
present-dayindividuals(Pru¨feretal.2014;Sankararamanetal. p
2014; Vernot and Akey 2014; Sankararaman et al. 2016; micdatafromphase3ofthe1000GenomesProject(Auton ril 2
et al. 2015), to detect AI in human populations, and find 0
Vernotetal.2016).Researchershavealsodetectedwhichof candidategenes.Whilethesestatisticsaresensitivetoadap- 19
theseintrogressedregionsarepresentathighfrequenciesin
tive introgression, they may also be sensitive to other phe-
certain present-day non-African populations. Someof these
nomena that generate genomic patterns similar to those
regions are likely to have undergone positive selection in
generated by AI, like ancestral population structure and in-
thosepopulationsaftertheywereintrogressed,aphenome-
complete lineage sorting. These processes, however, should
non known as adaptive introgression (AI). One particularly
not generate long regions of the genome where haplotypes
striking example of AI is the gene EPAS1 (Hu et al. 2003)
fromthesourceandtherecipientpopulationarehighlysim-
which confers a selective advantage in Tibetans by making
ilar.Asadditionalconfirmationthatthecandidateswefound
themlesspronetohypoxiaathighaltitudes(Bealletal.2010;
withourstatisticsaregeneratedbyAI,weexploredthehap-
Bighametal.2010;Yietal.2010;Pengetal.2011;Wangetal.
lotype structureofsomeofthemost promising candidates,
(cid:2)TheAuthor2016.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionNon-CommercialLicense
(http://creativecommons.org/licenses/by-nc/4.0/),whichpermitsnon-commercialre-use,distribution,andreproductioninany Open Access
medium,providedtheoriginalworkisproperlycited.Forcommercialre-use,[email protected]
296 Mol.Biol.Evol.34(2):296–317 doi:10.1093/molbev/msw216 AdvanceAccesspublicationOctober18,2016
MBE
ArchaicAdaptiveIntrogressioninPresent-DayHumanPopulations . doi:10.1093/molbev/msw216
Table1. Summarystatisticsmentionedinthemaintext.
Statistic Explanation References
D D-statistic:measuresexcessallelesharingbetweenatestpopulationandanoutgroupusinga (Greenetal.2010;
sisterpopulationthatismorecloselyrelatedtothetestthantheoutgroup Durandetal.2011)
f SimilartotheD-statistic,butservestobettercontrolforlocalvariationindiversitypatternsifone (Martinetal.2015)
D
isinterestedinfindinglociwithexcessancestryfromanadmixingpopulation.
R Averageratioofthesequencedivergencebetweenanindividualfromthesourcepopulationand Thisstudy
D
anindividualfromtheadmixedpopulation,andthesequencedivergencebetweenanindi-
vidualfromthesourcepopulationandanindividualfromthenon-admixedpopulation.Thisis
computedbytakingtheaverageoverallpairsofadmixedandnon-admixedindividuals.
UA;B;Cðw;x;yÞ NumberofsitesinwhichanyalleleisatafrequencylowerthanwinpanelA,higherthanxinpanel Thisstudy
BandequaltoyinpanelC.
UA;B;C;Dðw;x;y;zÞ NumberofsitesinwhichanyalleleisatafrequencylowerthanwinpanelA,higherthanxinpanel Thisstudy
B,equaltoyinpanelCandequaltozinpanelD.
D
Q95A;B;Cðw;yÞ 95%quantileofthedistributionofderivedallelefrequenciesinpanelB,forsiteswherethe Thisstudy ow
derivedalleleisatafrequencylowerthanwinpanelAandequaltoyinpanelC. n
lo
Q95A;B;C;Dðw;y;zÞ 95%quantileofthedistributionofderivedallelefrequenciesinpanelB,forsiteswherethe Thisstudy ad
derivedalleleisatafrequencylowerthanwinpanelA,equaltoyinpanelCandequaltozin e
d
panelD. fro
p Expectedheterozygosity,measuredastheaverageof2pð1(cid:2)pÞoverallsitesinawindow,wherep (Crowetal.1970) m
isthefrequencyofanarbitrarilychosenallele. h
D0½intro(cid:3) AmfreeqausuernecyofolfinhkaapgloetdyipseeqXuYil,ipbXriuismth.eCofrmeqpuuetnedcyaosfDa=llDelmeaxX,wphYeirsetDhe¼frepqXYue(cid:2)ncpyXpofY,apllXeYleisYt,haend (Lewontin1964) ttps://a
DmaxisthemaximumtheoreticalvaluethatDcantake.D0½intro(cid:3)iscomputedonlyusing cad
frequenciesfromtheintrogressedpanel.Whenwecomputethisinawindow,thevalueis e
m
obtainedbytakingtheaverageoverallpairsofSNPs. ic
D0½comb(cid:3) D0computedusinghaplotypeandallelefrequenciesfromthecombinationoftheintrogressed (Lewontin1964) .o
u
andnon-introgressedpanels. p
.c
r2½intro(cid:3) Ameasureoflinkagedisequilibrium.ComputedasD2=ðp ð1(cid:2)p Þp ð1(cid:2)p ÞÞ.r2½intro(cid:3)is (HillandRobertson1968) o
X X Y Y m
computedonlyusingfrequenciesfromtheintrogressedpanel.Whenwecomputethisina /m
window,thevalueisobtainedbytakingtheaverageoverallpairsofSNPs. b
e
r2½comb(cid:3) r2computedusinghaplotypeandallelefrequenciesfromthecombinationoftheintrogressed (HillandRobertson1968) /a
andnon-introgressedpanels. rtic
le
-a
b
s
tra
andusedaprobabilisticmethod(Seguin-Orlandoetal.2014) Second,underAIwewouldexpectalargenumberofsites ct/3
thatinfersintrogressedsegmentsalongthegenomebylook- containing archaic alleles at high frequency in the admixed 4/2
ing at the spatial arrangement of SNPs that are consistent population,butabsentoratlowfrequencyinanonadmixed /2
9
6
withintrogression.Thisallowsustoverifythatthecandidate population. Therefore, wedefine thestatistic UA;B;Cðw;x;yÞ /2
6
regionscontainintrogressedhaplotypesathighfrequencies:a tobeequaltothenumberofsiteswithinagenomicwindow 3
3
3
hallmarkofAI. where a sample C from an archaic source population has a 7
1
particularalleleatfrequencyy,andthatalleleisatafrequency b
y
smallerthanwinapanelAofanonadmixedpopulationbut g
Results u
e
largerthanxinapanelBofanadmixedpopulation(fig.1). s
StatisticsforDetectingAI Throughout the text, we denote panels A, B and C as the t on
0
Webeganbyevaluatingtheperformanceofvariousstatistics “outgroup”,“target”,and“bait”panels,respectively.Ifwehave 8
A
atdetectingAI.Inadditiontotestingstatisticsthathaveal- samples from two different archaic populations (e.g., a p
readybeenpreviouslydefinedintheliterature(D,f ,D0,r2,p), NeanderthalgenomeandaDenisovagenome),wecandefine ril 2
D 0
wedefinethreenewtypesofstatisticsthatwefindarepartic- UA;B;C;Dðw;x;y;zÞ asthenumberofsiteswherethearchaic 19
ularlypowerfulatdetectingAI(table1).Webrieflydescribe sampleChasaparticularalleleatfrequencyyandthearchaic
thenewstatisticshere,butmoreextensivedescriptionsofall sampleDhasthatalleleatfrequencyz.Additionally,atthose
testedstatisticscanbefoundinthe“Methods”sectionbelow. sites,thesamealleleshouldbeatafrequencysmallerthanw
First,inaregionunderAI,onewouldexpectthesequence inanoutgrouppanelAandlargerthanxinatargetpanelB
divergencebetweenanindividualfromthesourcepopulation (supplementaryfig.S1,SupplementaryMaterialonline).
andanadmixedindividualtobesmallerthanthesequence Finally,ifwedonotwanttosetahardcutoffforwhatwe
divergencebetweenanindividualfromthesourcepopulation consider “high-frequency” archaic alleles, we can just com-
andanonadmixedindividual.Thus,wedefineR astheav- pute a summary statistic of the site-frequency spectrum in
D
erageratioofthesequencedivergencebetweenapanelfrom thetargetpanel,conditionalonthearchaicallelebeingatlow
thesourcepopulationandapanelfromtheadmixedpopu- frequencyintheoutgroup.Thisstatisticshouldbehighwhen
lation, over the sequence divergence between the source aregioncontainsmanyallelesatespeciallyhighfrequenciesin
panel and a nonadmixed panel, computed in a window of thetarget.WethereforedefineQ95 ðw;yÞtobeequalto
A;B;C
thegenome. the95thpercentileofderivedfrequenciesinatargetpanelB
297
MBE
Racimoetal. . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
tra
c
FIG.1. SchematicillustrationofthewaytheUA;B;CandQ95A;B;Cstatisticsarecalculated. t/3
4
/2
/2
9
ofallSNPsthathaveaderivedallelefrequencyyinthebait settingtheminimumfrequencyofthearchaicalleleinthetest 6/2
6
(archaic) panel C, but where the derived allele is at a fre- populationatx¼0%meansthatanysitewithsomearchaic 3
3
quencysmallerthanwinanoutgrouppanelAfromanonad- alleleinthetestpanelwillbecounted,regardlessoftheallele 37
1
mixedpopulation(fig.1). frequency,solongasthearchaicalleleisatlowfrequencyin b
y
the outgroup panel. At high admixture rates, low- and g
u
SimulationsunderAI medium-frequency archaic alleles would naturally occur un- es
Weusesimulationstoassesstheperformanceofthestatistics derneutrality,sotheywouldnotbeinformativeaboutAI. t on
mentionedaboveatdetectingAI.SupplementaryfiguresS3– WealsoevaluatedtheeffectivenessofLD-basedstatistics 08
A
S5, Supplementary Material online show the distribution of at detecting AI (supplementary fig. S6, Supplementary p
statisticsthatrelyonpatternsofsharedalleleconfigurations Material online). We tested the performance of D0 and r2 ril 2
0
between source and introgressed populations (p, D, fD, by either computing each of these in the admixed panel 19
U ; Q95 ,andR ),fordifferentchoicesoftheselection only (D0½intro(cid:3); r2½intro(cid:3)) or in the combination of the ad-
A;B;C A;B;C D
coefficient s, and under 2%, 10%, and 25% admixture rates, mixed and nonadmixed panels (D0½comb(cid:3); r2½comb(cid:3)).
respectively. For Q95 ðw;100%Þ and U ðw;x;100%Þ, WhereasD0½intro(cid:3); D0½comb(cid:3) andr2½comb(cid:3) aremodestly in-
A;B;C A;B;C
we tested different choices of the outgroup cutoff w (1%, creasedbyAI,thisisnotthecasewithr2½intro(cid:3)understrong
10%)andthetargetcutoffx(0%,20%,50%,and80%). selectionandadmixtureregimes.Thisisbecauser2willtend
Some statistics, like Q95 ð1%;100%Þ and f show todecreaseiftheminorallelefrequencyisverysmall,which
A;B;C D
strongseparationbetweentheselectionregimes.Forexample, willoccurifthisfrequencyisonlymeasuredinthepopulation
withanadmixturerateof2%,Q95 ð1%;100%Þhas100% undergoingAI.Ingeneral,thesestatisticsarenotaspowerful
A;B;C
sensitivityataspecificityof99%,forboths¼0.1ands¼0.01. fordetectingAIasalleleconfigurationstatisticslikeUorQ95.
SomeparameterizationsoftheUstatisticarenotaseffective, To jointly explore the power and specificity of all these
however. For example, U ð1%;0%;100%Þ shows some statistics, we generated receiving operating characteristic
A;B;C
power when the admixture rate is low (2%), but almost no (ROC)curvesundervariousselectionandadmixtureregimes
powerwhentheadmixturerateishigh(25%).Thisisbecause (fig. 3 and supplementary fig. S7, Supplementary Material
298
MBE
ArchaicAdaptiveIntrogressioninPresent-DayHumanPopulations . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
FIG.2. Demographicmodelsdescribedinthemaintext. trac
t/3
4
/2
online).Ingeneral,Q95 ð1%;100%Þ; Q95 ð10%;100%Þ, increases with increasing selection intensity and admixture /2
A;B;C A;B;C 9
6
and fD are very powerful statistics for detecting AI under rates,pincreaseswithincreasingadmixturerates,butdecreases /26
strong(s¼0.1)andintermediate(s¼0.01)selectionpres- with increasing selection intensity. Thus, under AI the two 3
3
3
sures.ThenumberofuniquelysharedsitesUA;B;Cðw;x;yÞis forcescanceleachotherout,andweobtainasimilarvalueof 71
alsopowerful,solongasthepopulationinthetargetpanel pasunderneutrality.Furthermore,thejointdistributionsofQ b
y
(B)islarge.Additionally,fordifferentchoicesofx,usingw 95A;B;Cð1%;100%Þ andfDorRDshowparticularlygoodsepa- gu
e
¼1%yieldsamorepowerfulstatisticthanusingw¼10%. rationamongthedifferentAIscenarios. s
t o
WealsotestedAIscenarioswithweakselection(s¼0.001), Anotherjointdistributionthatisespeciallygoodatsepa- n
0
in which all statistics performed rather poorly, with Q95 ratingdifferentAIregimesisthecombinationofQ95 ðw; 8
A;B;C A
andf performingcomparablybetterthantherest(supple 100%ÞandU ðw;x;100%Þ.Infigure4,weshowthisjoint p
mentDaryfig.S8,SupplementaryMaterialonline).However, distribution, fAo;Br;Cdifferent choices of w (1% and 10%) and ril 2
0
undertheseconditions,itisveryunlikelythattheselected x(20%and50%).Here,withincreasingintensityofselection 19
allelewillreachappreciablefrequencies(supplementaryfig. andadmixture,thenumberofuniquelysharedsitesandthe
S2,SupplementaryMaterialonline),sothelackofsensitiv- quantilestatistic increase,butthequantilestatistic tendsto
ityofallstatisticshereislargelyareflectionofthefactthat only reach high values when selection is strong, even if ad-
in most simulations the selected allele is not successful, mixtureratesarelow.
especially when the probability of admixture is low.
Conditioning on survival of the allele should therefore in- AlternativeDemographicScenarios
creasesensitivity. Weevaluatedtheperformanceofourstatisticsundervarious
We were also interested in the joint distribution of pairs alternative demographic scenarios. First, we simulated a 5X
of these statistics. Supplementary figure S9, Supplementary bottleneckoccurringinpopulationB1600generationsbefore
Material online shows the joint distribution of Q95 ð1%; theadmixtureevent,andlasting200generations,toobserve
A;B;C
100%Þinthey-axisandfourotherstatistics(R ,p,D,andf )in itseffectsonthepowerofthestatisticsfordetectingAI(fig.
D D
thex-axis,underdifferentadmixtureandselectionregimes.One 2B).Thoughweobserveareductioninpower—mostevident
can observe, for example, that whereas Q95 ð1%;100%Þ intheheterozygositystatistics—noneofthestatisticsarevery
A;B;C
299
MBE
Racimoetal. . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
tra
c
t/3
4
/2
/2
9
6
/2
6
3
3
3
7
1
b
y
g
u
e
s
t o
n
0
8
A
p
FIG.3.Receiveroperatingcharacteristiccurvesforascenarioofadaptiveintrogression(s¼0.1)comparedwithascenarioofneutrality(s¼0), ril 2
0
using1,000simulationsforeachcase.PopulationsAandBsplitfromeachother4,000generationsago,andtheirancestralpopulationsplitfrom 1
9
populationC16,000generationsago.Populationsizeswereconstantandsetat2N¼20;000.Theadmixtureeventoccurred1,600generationsago
frompopulationCtopopulationB,atrate2%(toppanels)or25%(bottompanels).Therightpanelsarezoomed-inversionsoftheleftpanels.
strongly affected by this event (supplementary fig. S10, distribution under neutrality, which explains the low power
SupplementaryMaterialonline).Wealsosimulatedabottle- (supplementary figs. S12 and S13, Supplementary Material
neckofequalsizebutoccurringaftertheadmixtureevent— online). Additionally, U (but not Q95) seems to display
A;B;C
starting 1,400 generations ago, and lasting 200 generations more elevated values under neutrality in the bottleneck
(fig. 2C). In this case, the sensitivity of all the statistics is model than in the constant population size model.
stronglyreducedwhentheadmixturerateislow(supplemen However, the relative performance of each statistic with re-
tary fig. S11, Supplementary Material online). For example, specttoalltheothersdoesnotappeartochangesubstantially
when looking at the raw values of the U and Q95 (supplementaryfig.S11,SupplementaryMaterialonline).
A;B;C A;B;C
statistics,weobservethatforlowadmixtureratesthedistri- We next explored a model where the introgressed hap-
bution under selection has a larger overlap with the lotype was not immediately adaptive in the Eurasian
300
MBE
ArchaicAdaptiveIntrogressioninPresent-DayHumanPopulations . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
tra
c
t/3
4
/2
/2
9
6
/2
6
3
3
3
7
1
b
y
FIG.4. JointdistributionofQ95A;B;Cðw;yÞandUA;B;Cðw;x;yÞfordifferentchoicesofw(1%,10%)andx(20%,50%).Wesetyto100%inallcases.100 gue
individualsweresampledfrompanelA,100frompanelB,and2frompanelC.Thedemographicparameterswerethesameasinfigure3. s
t o
n
0
population,butinsteadunderwentanintermediateperiod the most sensitive ones (supplementary fig. S15, 8
A
ofneutraldrift,beforeitbecameadvantageous(fig.2D).In SupplementaryMaterialonline).Thisisexpected,giventhat p
such a situation, our power to detect AI is reduced, for all classicalselectivesweepsareknowntostronglyaffectpatterns ril 2
0
statistics (supplementary fig. S14, Supplementary Material ofheterozygosityandlinkagedisequilibriumintheneighbor- 19
online).Thisisparticularlyanissuewhentheadmixturerate hoodoftheselectedsite(Barton1998;KimandStephan2002;
islow,asinthosecasesthestartingfrequencyoftheselected Kim and Nielsen 2004). Since all other statistics have very
alleleintheEurasianpopulationislow,soitismorelikelyto poor sensitivity todetect SDN, weexpect tobeabletodis-
drift to extinction during the neutral period, before it can tinguishsignaturesgeneratedfromSDNandAI.Onecaveat
becomeadvantageous. tothisisthescenarioinwhichadenovoselectedmutation
Wealsoevaluatedtheperformanceofourstatisticsunder occurs on an introgressed haplotype immediately after an
selectivescenariosthatdidnotinvolveadaptiveintrogression, introgression event—before the haplotype has a chance to
tocheckwhichofthemweresensitivetothesemodelsand expandandrecombineinthepopulation—inwhichcaseour
which were not. Under a model of selection from de novo statisticswillnotbeabletodistinguishSDNfromAI.
mutation(SDN,fig.2E)—inwhichasinglemutationappears Wealsosimulatedamodelofselectionfromstandingvar-
inthereceivingpopulationafterthesplittimebetweenitand iation(fig.2F),byrandomlyselecting20%ofhaplotypeswithin
the nonadmixed population—the heterozygosity (p) and theintrogressedpopulationtobeadvantageous,afterthesplit
linkage disequilbrium statistics (r2½intro(cid:3) and D0½intro(cid:3)) are time between it and the nonintrogressed population. In this
301
MBE
Racimoetal. . doi:10.1093/molbev/msw216
case,allstatisticsperformpoorly,especiallywhenadmixtureis shared archaic alleles than Europeans and South Asians, for
low.Interestingly,whenadmixtureishigh(supplementaryfig. both x¼20% and x¼50% (fig. 5B–C). Population-specific
S16,SupplementaryMaterialonline),Q95 ð1%;100%Þand D-statistics (using YRI as the nonadmixed population) also
A;B;C
U ð1%;0%;100%Þarethebestperformingstatistics.Thisis follow this trend (supplementary fig. S20, Supplementary
A;B;C
likelybecausesomeofthehaplotypesthatarerandomlycho- Material online) and we observe this pattern when looking
sentobeselectedalsohappentobeancestrallypolymorphic only at the X chromosome as well (supplementary fig. S21,
andpresentinthearchaichumans. Supplementary Material online). These results hold in com-
Whenwesetancestralstructuretobeournullmodel,we parisonswithbotharchaichumangenomes,butweobservea
observedifferentbehaviorsdependingonthestrengthofthe stronger signal when looking at Neanderthal-specific shared
migrationrates.Whenthemigrationratesarestrong(supple alleles. To correct for the fact that some panels have more
mentary fig. S17, Supplementary Material online), we have segregatingsitesthanothers(andmaythereforehavemore
excellentpowertodetectAIwithseveralstatistics,including archaic-like segregating sites), we also scaled the number of
Q95A;B;Cð1%;100%Þ, D, fD, RD, and UA;B;Cð1%;50%;100%Þ. uniquelysharedsitesbythetotalnumberofsegregatingsites D
o
Whentheratesareofmediumstrength(supplementaryfig. per population panel (fig. 5D–F), and we see in general the w
n
S18, Supplementary Material online), the power is slightly samepatterns,withtheexceptionofaPeruvianpanel,which lo
a
reduced, but the same statistics are the ones that perform we discuss further below. We also observe similar patterns de
d
best.Whenthemigrationratesareweak—meaningancestral when calculating Q95A;B;Nea;Denð1%;y;zÞ genome-wide (sup fro
structure is very strong—Q95 ð1%;100%Þ loses power, plementaryfig.S22,SupplementaryMaterialonline).Theel- m
A;B;C h
and the best-performing statistics are RD, D, and fD (supple evation in UA;B;Nea;Den and Q95A;B;Nea;Den in East Asians may ttp
mentary fig. S19, Supplementary Material online). We note, possibly result from higher levels of archaic ancestry in East s://a
though, that the genome-wide D observed under this last Asians than in Europeans (Wall et al. 2013), which some c
a
d
ancestralstructuremodel(D¼0.24)ismuchmoreextreme studies argue could be due to additional admixture events e
m
thanthegenome-wideDobservedempiricallybetweenany occurring in East Asians (Kim and Lohmueller 2015; Vernot ic
.o
Eurasian population and Neanderthals or Denisovans, sug- andAkey2015;). u
p
gestingthatiftherewasancestralstructurebetweenarchaic Surprisingly,thePeruvians(PEL)harborthelargestamount .c
o
m
andmodernhumans,itwaslikelynotofthismagnitude. ofhighfrequencymutationsofarchaicoriginthananyother /m
singlepopulation,especiallywhenusingNeanderthalsasbait b
e
GlobalFeaturesofUniquelySharedArchaicAlleles (figs. 5B–C and supplementary fig. S21, Supplementary /a
Before identifying candidate genes for adaptive introgres- Material online). It is unclear whether this signal is due to rticle
sion,weinvestigatedthefrequencyandnumberofuniquely increaseddriftorselectioninthispopulation.Skoglundand -a
b
s
shared alleles at the genome-wide level. Specifically, we Jakobsson(2011)argueviasimulationsthatifoneanalyzesa tra
wanted to know whether human populations varied in population with high amounts of recent genetic drift and ct/3
the number of sites with uniquely shared archaic alleles, excludes SNPs where the minor allele is at low frequency, 4/2
andwhether they alsovaried in thefrequency distribution some statistics that are meant to detect archaic ancestry— /2
9
6
of these alleles. Therefore, we computed UA;B;Nea;Den likeD—maybeartificiallyinflated.Ourfilteringprocedureto /2
6
(1%,x,y,z)andQ95 (1%,y,z)fordifferentchoicesof selectuniquelysharedarchaicallelesnecessarilyexcludessites 3
A;B;Nea;Den 3
3
x,y,andz.Weuseddifferentcutoffsforthefrequencyofthe wherethearchaicalleleisatlowfrequencyinthetargetpanel, 7
1
archaic allele (x) in the target population B: 0%, 20%, and andthePELpanelcomesfromapopulationwithahistoryof b
y
50%. To look for alleles uniquely shared with the Altai loweffectivepopulationsizes(highdrift)relativetoothernon- g
u
e
Neanderthal genome only, we set y¼100% and z¼0%. Africans(Autonetal.2015),whichcouldexplainthispattern. s
t o
To look for alleles uniquely shared with the Denisovan ge- Thiscouldalsoexplainwhytheeffectisnotseenwhenx¼0% n
0
nomeonly,wesety¼0%andz¼100%.Finally,tolookfor (fig. 5A), or when computing D-statistics (supplementary fig. 8
A
uniquely shared alleles matching both of the archaic ge- S20, Supplementary Material online), both of which include p
nomes,wesety¼100%andz¼100%. sites with low-frequency alleles in their computation. ril 2
0
We used each of the non-African panels in the 1000 Additionally,scalingtheuniquelysharedsitesbythetotalnum- 19
Genomes Project phase 3 data (Auton et al. 2015) as the ber of segregating sites per population panel mitigates (but
“target” panel (B), and chose the outgroup panel (A) to be does not completely erase) this pattern. After scaling, PEL
thecombinationofallAfricanpopulations(YRI,LWK,GWD, shows levels of archaic allele sharing within the range of the
MSL, and ESN), excluding admixed African-Americans. We EastAsianpopulationsatx¼20%(fig.5E),butisstillthepanel
note that this is a conservative reference panel, as some of withthelargestnumberofarchaicsitesatx¼50%(fig.5F).
the African panels—like LWK—are from populations with a Furthermore,weplottedthevaluesofU ð1%;x;
AFR;X;Nea;Den
substantial amount of Eurasian ancestry (Auton et al. 2015), y;zÞandQ95 ð1%;y;zÞjointlyforeachpopulation
AFR;X;Nea;Den
likely preventing the detection of introgressed segments at X,underdifferentfrequencycutoffsx.Whenx¼0%,thereis
someloci. a generally inversely proportional relationship between the
Whensettingx¼0%(i.e.,notimposingafrequencycutoff two scores (supplementary fig.S23, Supplementary Material
in the target panel B), South Asians as a target population online),butthisbecomesadirectlyproportionalrelationship
showthelargestnumberofarchaicalleles(fig.5A).However, when x¼20% (fig. 6) or x¼50% (supplementary fig. S27,
EastAsianshavealargernumberofhigh-frequencyuniquely SupplementaryMaterialonline).Here,wealsoclearlyobserve
302
MBE
ArchaicAdaptiveIntrogressioninPresent-DayHumanPopulations . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
tra
c
t/3
4
/2
/2
9
6
/2
6
3
3
3
7
1
b
y
g
FIG.5.WecomputedthenumberofuniquelysharedsitesintheautosomesandtheXchromosomebetweenparticulararchaichumansand ue
differentchoicesofpresent-daynon-AfricanpanelsX(x-axis)fromphase3ofthe1000GenomesProject.Weusedasharedfrequencycutoffof0% st o
(A), 20% (B), and 50% (C). Nea-only¼ U ð1%;20%;100%;0%Þ. Den-only¼ U ð1%;20%;0%;100%Þ. Nea-all¼ n
Afr;X;Nea;Den Afr;X;Nea;Den 0
U ð1%;20%;100%Þ.Den-all¼ U ð1%;20%;100%Þ.Both¼ U ð1%;20%;100%;100%Þ.Finally,wescaledeachofthesta- 8
Afr;X;Nea Afr;X;Den Afr;X;Nea;Den A
tisticsfrompanelsAtoCbythenumberofsegregatingsitesineach1000Genomespopulationpanel,yieldingpanelsD–F. p
ril 2
0
1
9
thatPELisanextremepanelwithrespecttoboththenumber analysesofthe1000GenomessamplessuggestthatPeruvians
andfrequencyofarchaicsharedderivedalleles,andthatEast havethelargestamountofNativeAmericanancestry(Auton
Asianand Americanpopulationshavemorehigh-frequency et al. 2015) and show a bottleneck with a lack of recent
archaicsharedallelesthanEuropeans. population growth, which could explain this pattern. We
Wecheckedviasimulationsiftheobservedexcessofhigh also observe an increase in the variance of the distribution
frequencyarchaicderivedmutationsinAmericansandespe- ofUandQ95inthepresenceofabottleneck,especiallywhen
ciallyPeruvianscouldbecausedbygeneticdrift,asaconse- long and severe (supplementary figs. S25 and S26,
quence of the bottleneck that occurred in the ancestors of SupplementaryMaterialonline).
NativeAmericansastheycrossedBeringia.Weobservethatif
the introgressed population B undergoes a bottleneck, this CandidateRegionsforAdaptiveIntrogression
can lead to a larger number of UA;B;Cðw;x;y;zÞ for large Toidentifyadaptivelyintrogressedregionsofthegenome,we
values of x (supplementary figs. S12, S13, and S24, computedU ðw;x;y;zÞandQ95 ðw;y;zÞin40kb
A;B;C;D A;B;C;D
SupplementaryMaterialonline).Indeed,populationstructure nonoverlapping windows along the genome, using the
303
MBE
Racimoetal. . doi:10.1093/molbev/msw216
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
-a
b
s
tra
c
t/3
4
/2
/2
9
6
/2
6
3
FIG.6. Foreachpopulationpanelfromthe1000GenomesProject,wejointlyplottedtheUandQ95statisticswithanarchaicfrequencycutoffof 33
>20% within each population. Nea-only¼ UAfr;X;Nea;Denð1%;20%;100%;0%Þ and Q95Afr;X;Nea;Denð1%;100%;0%Þ. Den-only¼ 71
UUAAffrr;;XX;;NDeean;ðD1en%ð;12%0;%2;01%0;00%%Þ;a1n0d0%QÞ95anAfdr;XQ;D9en5ðA1fr%;X;;N1ea0;0D%enÞð.1B%o;t0h%¼;10U0A%fr;ÞX.;NNeae;Dae-naðll1¼%;U2A0f%r;X;;N1e0a0ð%1%;1;0200%%Þ;1a0n0d%QÞ9a5nAdfr;QX;9Ne5aA;Dfre;Xn;ðN1e%að;11%0;01%0;01%0Þ0.%DÞe.n-all¼ by gue
s
t o
low-coverage sequencing data from phase 3 of the 1000 all the CCDS-verified genes located inside each window n
0
GenomesProject(Autonetal.2015).Weusedthiswindow (Pruittetal.2009). 8
A
sPirzu¨efebrecaeutseatlh.em(20ea1n4)lenwgaths of44in,0t7ro8gbrpesse(dSuhpapplleomtyepnetsairny woFuoldregxupideacntcuendaesrtnoeuhtorwalithyi,gwheauvseadlutehoefsiUmualnadtioQn9s5f,rowme pril 2
0
1
Information13).Ourmotivationwastofindregionsunder figure2toobtain95%empiricalquantilesofthedistribution 9
AIinaparticularpanelB,usingpanelAasanonintrogressed ofthesescoresunderneutrality.SupplementarytablesS1and
outgroup (generally Africans, unless otherwise stated). We S2, Supplementary Material online show the 95% quantiles
used the high-coverage Altai Neanderthal genome (Pru¨fer forthesetwostatisticsundervariousmodelsofadaptivein-
etal.2014)asbaitpanelCandthehigh-coverageDenisova trogression and ancestral structure, for different choices of
genome (Meyer et al. 2012) as bait panel D. We deployed parametervalues(see“Methods” section).Whenexamining
these statistics in three ways: (a) to look for Neanderthal- ourcandidatesforAIbelow,wefocusedonwindowswhose
specific AI, we set y¼100% and z¼0%; (b) to look for values for U ðw;x;y;zÞ and Q95 ðw;y;zÞ
A;B;Nea;Den A;B;Nea;Den
Denisova-specific AI, we set y¼0% and z¼100%; (c) to were both in the 99.9% quantile of their respective
look for AI matching both of the archaic genomes, we set genome-widedistributions.Weverifiedviasimulationsthat,
y¼100% and z¼100% (supplementary fig. S1 and table underasimplemodelofneutraladmixtureatagenome-wide
S3,SupplementaryMaterialonline).Totrytodeterminethe rate of 2%, the estimated probability of obtaining values as
adaptivepressurebehindtheputativeAIevent,weobtained high as these (or the false positive rate, FPR) was between
304
MBE
ArchaicAdaptiveIntrogressioninPresent-DayHumanPopulations . doi:10.1093/molbev/msw216
10.6% and 0%, depending on the target population chosen. frequency, which were absent or almost absent in
ThehighestratescorrespondtotheAfrican-American(ASW) Europeans,EastAsiansandAfricans.Inthiscase,wefocused
admixed panel, as this panel contains high proportions of onx¼20%becausewefoundthatx¼50%leftuswithno
ancestry from the outgroup panel (unadmixed Africans) candidate regions. Among the candidate regions sharing a
and are therefore not well-suited for our statistics. large number of high-frequency Neanderthal alleles in South
ExcludingASW,thehighestestimatedFPRwas5.5%. Asians,wefindgenesASTN2,SFMBT1,MUSTN1andMAML2
We also calculated D and f along the same windows (supplementary fig. S30, Supplementary Material online).
D
(using Africans as the nonadmixed population), and saw ASTN2isinvolvedinneuronalmigration(Wilsonetal.2010)
goodagreementwiththenewstatisticspresentedhere(sup and is associated with schizophrenia (Vrijenhoek et al. 2008;
plementarytableS3,SupplementaryMaterialonline).Finally, Wangetal.2010).SFMBT1isinvolvedinmyogenesis(Linetal.
we further validated the regions most likely to have been 2013)andisassociatedwithhydrocephalus(Katoetal.2011).
adaptively introgressed by searching for archaic tracts of in- MUSTN1playsaroleintheregenerationofthemuscoskeletal
trogressionwithinthemthatwereathighfrequency,usinga system(Krauseetal.2013).Finally,MAML2codesforasignal- D
o
HiddenMarkovModel(seebelow). ingprotein(Linetal.2002;Wuetal.2005),andisassociated w
n
with cutaneouscarcinoma(Winnes etal. 2007) and lacrimal lo
a
glandcancer(VonHolsteinetal.2012). de
d
ContinentalPopulations fro
Whenfocusingonadaptiveintrogressionincontinentalpop- m
Eurasia h
ulations, we first looked for uniquely shared archaic alleles We then looked for AI in all Eurasians (EUA¼EURþ ttps
specific to Europeans that were absent or almost absent SASþEAS, ignoring American populations) using Africans ://a
(< 1% frequency) in Africans and East Asians. In addition, c
as the nonadmixed population (AFR, ignoring admixed ad
we also looked for uniquely shared archaic alleles in East e
African-Americans). Figure 8 shows the extreme outlier re- m
Asians,whichwereabsentoralmostabsentinAfricansand ic
Europeans. In this continental survey, we ignored Latin gionsthatareinthe99.9%quantilesforbothUEUA;AFR;Nea;Den .ou
ð1%;20%;y;zÞandQ95 ð1%;y;zÞ,whereassup p
American populations as they have high amounts of EUA;AFR;Nea;Den .c
plementaryfigureS31,SupplementaryMaterialonlineshows om
European and African ancestry, which could confound our theentiredistribution.Wefocusedonx¼20%becausewe /m
analyses. Figure 7 shows the number of sites with uniquely b
shared alleles for increasing frequency cutoffs in the intro- foundthatx¼50%leftuswithalmostnocandidateregions. e/a
g(Nreesasneddeprtohpaul-lasptieocnif,ica,nDdefnoirsodviaff-esrpeenctifitcypoerscoofmamrcohnaictoablleoltehs Iunnitqhuiselycasshea,rtehdearrecghioainc waliltehlesbyisftahrethoenelacrgoensttainnuinmgbgeerneosf rticle-a
OAS1 andOAS3,involvedininnateimmunity (Knappetal. bs
archaic humans). In other words, we calculated 2003;Hamanoetal.2005;Fedetzetal.2006;Limetal.2009). tra
UAFR;EUR;Nea;Denð1%;x;y;zÞ and UAFR;EAS;Nea;Denð1%;x;y;zÞ This region was previously identified as a candidate for AI ct/3
fordifferentvaluesofx(0%,20%,50%,and80%)anddifferent from Neanderthals in non-Africans (Mendez et al. 2013). 4/2
choicesofyandz,dependingonwhichtypeofarchaicalleles /2
Anotherregionthatwerecoverandwaspreviouslyidentified 96
wewerelookingfor.Weobservethattheregionsintheex- /2
as a candidate for AI is the one containing genes TLR1 and 6
treme of the distributions for x¼50% corresponded very 3
TLR6(Dannemannetal.2016;Deschampsetal.2016).These 33
welltogenesthathadbeenpreviouslyfoundtobecandidates 7
genes are also involved in innate immunity and have been 1
foradaptiveintrogressionfromarchaichumansinthesepop- b
ulations, using more complex probabilistic methods shown to be under positive selection in some non-African y g
populations(Akiraetal.2006;Barreiroetal.2009). ue
(Sankararamanetal. 2014; Vernot and Akey 2014) orgene- s
Interestingly,wefindthataverystrongcandidateregionin t o
centric approaches (Ding et al. 2013). These include BNC2 n
Eurasia contains genes TBX15 and WARS2. This region has 0
(involved in skin pigmentation [Vanhoutteghem and Djian 8
been associated with a variety of traits, including adipose A
2006;Jacobsetal.2013]),POU2F3(involvedinskinkeratino- p
cytedifferentiation[Cabraletal.2003;Takemotoetal.2010]), tissuedifferentiation(Gburciketal.2012),bodyfatdistribu- ril 2
tion (Heid et al. 2010; Liu et al., 2013, 2014; Shungin et al. 0
HYAL2(involvedintheresponsetoUVradiationonhuman 19
2015),hairpigmentation(Candilleetal.2004)facialmorphol-
keratinocytes[Ha(cid:3)sov(cid:2)aetal.2011]),SIPA1L2(involvedinneu-
ogy(Lauschetal.2008;Pallaresetal.2015),earmorphology
ronal signaling [Spilker and Kreutz 2010]), and CHMP1A (a
(Curry1959),stature(Lauschetal.2008),andskeletaldevel-
regulator of cerebellar development [Mochida et al. 2012]).
opment (Singh et al. 2005; Lausch et al. 2008). It was previ-
Tobemorerigorousinoursearchforadaptiveintrogression,
ously identified as being under positive selection in
welookedatthejointdistributionoftheUstatisticandthe
Greenlanders(Fumagallietal.2015),anditshowsparticularly
Q95 statistic for the same choices of w, y, and z, and then
striking signatures of adaptive introgression, so we devotea
selectedtheregionsthatwereinthe99.9% quantilesofthe
separatestudytoitsanalysis(Racimoetal.2016).
distributionsofbothstatistics(fig.8,supplementaryfigs.S28
and S29, Supplementary Material online). We find that the
strongest candidates here are BNC2, POU2F3, SIPA1L2, and Population-SpecificSignalsofAdaptiveIntrogression
theHYAL2region. To identify population-specific signals of AI, we looked for
WealsoscannedforregionsofthegenomewhereSouth archaic alleles at high frequency in a particular non-African
Asians (SAS) had uniquely shared archaic alleles at high panelX,whichwerealsoatlessthan1%frequencyinallother
305