Simultaneous purifying selection on the ancestral MC1R allele and positive selection on the ...

12 Pages·2013·0.52 MB·English

Simultaneous Purifying Selection on the Ancestral MC1R
Allele and Positive Selection on the Melanoma-Risk
Allele V60L in South Europeans
Conrado Mart´ınez-Cadenas,y,1 Saioa Lo´pez,y,2 Gloria Ribas,3 Carlos Flores,4,5 Oscar Garc´ıa,6 Arrate Sevilla,2
Isabel Smith-Zubiaga,7 Maider Ibarrola-Villaba,3 Maria del Mar Pino-Yanes,4,5,8 Jesu´s Gardeazabal,9
Dolores Boyano,10 Alicia Garc´ıa de Galdeano,10 Neskuts Izagirre,2 Concepcio´n de la Ru´a,2 and
Santos Alonso*,2
1DepartmentofMedicine,UniversityofCastello´n,Castello´n,Spain
2DepartmentofGenetics,PhysicalAnthropologyandAnimalPhysiology,UniversityoftheBasqueCountry,Bizkaia,Spain
3INCLIVABiomedicalResearchInstitute,Valencia,Spain
4ResearchUnit,UniversityHospitalN.S.deCandelaria,Tenerife,Spain Do
w
5CIBERonRespiratoryDiseases,ISCIII,Madrid,Spain n
lo
6ForensicScienceUnit,ForensicGeneticsSection,BasqueCountryPolice,Erandio,Bizkaia,Spain ad
e
d
7DepartmentofZoologyandAnimalCellBiology,UniversityoftheBasqueCountry,Bizkaia,Spain fro
8DepartmentofMedicine,UniversityofCalifornia,SanFrancisco,CA,USA m
h
9DermatologyService,CrucesHospital,Bizkaia,Spain ttp
s
10DepartmentofCellBiologyandHistology,UniversityoftheBasqueCountry,Bizkaia,Spain ://a
c
yTheseauthorscontributedequallytothiswork. ad
e
m
*Correspondingauthor:E-mail:[email protected].
ic
Associateeditor:NaruyaSaitou .ou
p
.c
o
m
Abstract /m
b
e
In humans, the geographical apportionment of the coding diversity of the pigmentary locus melanocortin-1 receptor /a
(MC1R)is,unusually,higherinEurasiansthaninAfricans.Thisatypicalobservationhasbeeninterpretedastheresultof rtic
le
purifying selection due to functional constraint on MC1R in high UV-B radiation environments. By analyzing 3,142 -a
b
humanMC1RallelesfromdifferentregionsofSpaininthecontextofadditionalhaplotypicinformationfromthe1000 s
tra
Genomes(1000G)Projectdata,weshowthatpurifyingselectionisalsostronginsouthernEurope,butnotsoinnorthern c
Europe.Furthermore,weshowthatpurifyingandpositiveselectionactsimultaneouslyonMC1R.Thus,atleastinSpain, t/30
/1
regions at opposite ends of the incident UV-B radiation distribution show significantly different frequencies for the 2
/2
A melanoma-risk allele V60L (a mutation also associated to red hair and fair skin and even blonde hair), with higher 6
5
r frequencyofV60LatthoseregionsoflowerincidentUV-Bradiation.Besides,usingthe1000GsouthEuropeandata,we 4/1
t 0
showthattheV60Lhaplogroupisalsocharacterizedbyanextendedhaplotypehomozygosity(EHH)patternindicativeof 1
i 4
c positive selection. We, thus, provide evidence for an adaptive value of human skin depigmentation in Europe and 64
l 5
e illustratehowanadaptiveprocesscansimultaneouslyhelptomaintainadisease-riskallele.Inaddition,ourdatasupport b
y
the hypothesis proposed by Jablonski and Chaplin (Human skin pigmentation as an adaptation to UVB radiation. g
u
ProcNatlAcadSciUSA.2010;107:8962–8968),whichpositsthathabitationofmiddlelatitudesinvolvedtheevolution es
ofpartiallydepigmentedphenotypesthatarestillcapableofsuitabletanning. t on
1
Keywords:skinpigmentation,MC1R,adaptiveevolution,melanoma. 2 A
p
ril 2
Introduction 0
melanins). In the absence of MC1R signaling, pheomelano- 19
The human melanocortin-1 receptor (MC1R) is an integral genesis (synthesis of yellow/red melanins) is the default
membrane G-protein-coupled receptor (GPCR) of skin pathway (Ito 2003). Therefore, MC1R is a regulator of the
and hair follicle melanocytes. GPCRs communicate the cells amount and type of pigment production, and has thus
with their environment and upon receiving external stimuli been referred to as a major determinant of skin phototype
trigger an adaptive response. In this regard, the binding of (Garc´ıa-Borro´netal.2005).Inthisregard,somenonfunctional
a small peptide, the a-melanocyte-stimulating hormone MC1Rvariantsleadtophenotypescharacterizedbyredhair,
(a-MSH) toMC1Riscoupledtoa cAMPsignaling pathway fair skin, freckles, and poor tanning ability (the red hair and
that stimulates eumelanogenesis (synthesis of brown/black fairskinphenotype,RHC)inadominantmanner.Thesehave
(cid:2)TheAuthor2013.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.Allrightsreserved.Forpermissions,please
e-mail:[email protected]
2654 Mol.Biol.Evol.30(12):2654–2665 doi:10.1093/molbev/mst158 AdvanceAccesspublicationSeptember17,2013
MBE
VaryingSelectiononMC1RAlleles . doi:10.1093/molbev/mst158
been called R alleles, andinclude D84E,R151C, R160W, and V60L,consistentwithotherstudiesonEuropeanpopulations
D294H.Variantsthatdonotleadtoatotallossoffunction (Flanaganetal.2000;Fargnolietal.2006;Latreilleetal.2009).
andthathaveweakornoassociationtotheRHCphenotype To this sample, we added those sequences previously pub-
(NRHC phenotype) are called r alleles, and include, for lished by some of us (Fernandez et al. 2007; Ibarrola-Villava
instance,V60L,V92M,andR163Q(Valverdeetal.1995;Box et al. 2012), which included 354 individuals from Madrid
etal.1997;Sturmetal.2003;Duffyetal.2004;Beaumontetal. (Spain). The global sample showed a total of 42 SNPs, with
2007;Dessiniotietal.2011). atotalof42differenthaplotypesobserved(seetable1fora
As skin phototype is a risk factor for melanoma (Abdel- comparativediversityanalysis).
MalekandIto2013),considerableefforthasbeendedicated WefurtherresequencedMC1Rin127melanomapatients
to find association between variants at MC1R (and other and detected 15 SNPs. The most frequent variant was also
pigmentation genes) to melanoma risk (Valverde et al. V60L.Tothissample,weaddedthosemelanomasequences
1996; Ichii-Jones et al. 1998; Palmer et al. 2000; Kennedy previously published by some of us (Fernandez et al. 2007;
et al. 2001; Fernandez et al. 2007; Raimondi et al. 2008; Ibarrola-Villavaetal.2012),whichcorrespondedto595indi- D
o
Scherer et al. 2009; Williams et al. 2011; Davies et al. 2012). viduals.Then,theglobalmelanomasampleshowedatotalof wn
Less effort has been dedicated to understand pigmentation 37differenthaplotypesformedfrom28SNPs,7ofwhichdid loa
d
and melanoma risk from an evolutionary perspective. notshowupinthecontrolsamples(S41F,S41C,G89R,I180I, ed
Interestingly,MC1Rishighlydiverseinhumans:arecentlit- R213W,P268R,andN279K). fro
m
erature review cataloged up to 57 nonsynonymous and 25 In our healthy sample, the most likely root is the most h
synonymous variants in approximately 1kb of coding DNA common haplotype. Its likelihood was approximately 2 ttp
s
(Gerstenblithetal.2007).However,thegeographicalappor- orders of magnitude higher than the one corresponding to ://a
c
tionment of this coding diversity shows that there is an in- the nextmostlikely rootedtree. The mostcommon haplo- a
d
creaseddiversityinEurasianpopulationsincomparisonwith typeinoursampleisalsothemostcommonhaplotypeinthe em
African populations (Rana et al. 1999; Harding et al. 2000), 1000GAfricansample.MutationT314T,whichliesinaCpG ic.o
whichisanunusualobservation,giventhatAfricanstypically doublet,wasremovedbecauseitseemstohavearisenmore up
.c
exhibitgreatergenomicdiversitythanEurasians.Thisiseven than once (within different lineages). Its frequency in the o
m
more intriguing as the upstream region containing the pro- healthy samples is 0.07 (total 2N=3,142 alleles). Mutation /m
b
moter of MC1R does show the expected, opposite pattern, R151C was also discarded because it seems to have arisen e
/a
that is, a higher diversity in Africans compared with Asians several times (within different lineages) in the melanoma rtic
andEuropeans(Makovaetal.2001). samples(frequencyinthehealthysamples:0.02). le
-a
Thesepatterns of diversity in the coding region ofMC1R b
s
ledbothRanaetal.(1999)andHardingetal.(2000)toposit Diversity Analysis tra
c
that MC1R is under strong functional constraint in Africa. t/3
ThecodingregionofMC1Rshowsahighnucleotidediversity 0
However, their interpretation for the diversity patterns in /1
for such a short coding region: We observed an estimated 2
Eurasia differed. Thus, while Rana et al. (1999) proposed /2
nucleotide diversity (cid:2) of 0.08%, whereas in a small set of 6
localselectionformutationsthatresultinlighterpigmenta- 54
tion (diversifying selection) in Eurasia, Harding et al. (2000) European human exomes the average (cid:2) was 0.038% /10
(Tennessen et al. 2010). This increased nucleotide diversity 1
proposed instead that the level of MC1R polymorphism 46
seemstobetheresultofahighmutationratederivedfroma 4
simply reflects the expectations of a neutral evolution for 5
high CpG content. Figure 1 shows the CpG content of the b
thislocusinEurasians. y
codingregionofMC1R.Forthisregion,butnotfortheflank- g
In this work, we propose an explanation for the high ue
ingregions,theobservedversusexpectedCpGratioishigher s
diversity at MC1R and assess the role of selection in one t o
than 0.65,whichqualifiesthisregionasaCpG island(Takai n
large Spanish population sample in the European context. 1
andJones2002).ThishighCpGcontentisalsocorrelatedwith 2
We conclude that purifying and positive selection can act A
thehighest(cid:2),butthelowestTajima’sDvalues. p
simultaneously on the same locus. More interestingly, ril 2
MC1Rillustrateshowanadaptiveprocesscansimultaneously 0
1
helpmaintainadisease-riskallele. Evidence for Purifying Selection at MC1R 9
TherearemanymutationsatcodingMC1R,butthesedonot
Results
reachingeneralasubstantialfrequencyinoursample(sup-
AfterresequencingMC1Rin1,217healthySpanishindividuals, plementaryinformationS1,SupplementaryMaterialonline).
wedetected38single-nucleotidepolymorphism(SNPs),6of AnexceptioncouldbeV60L,anrallelewithafrequencyof
which have not been previously reported: V193M, A285T, approximately 0.15, which will be dealt as discussed later.
R229H,Q233X,S302N,andW317C.Wedidnotobservethe Neutrality tests indicate that the low frequency of new
R307G Neanderthal MC1R mutation described by Lalueza- mutantscanbeduetotheeffectofpurifyingselectionthat
Fox et al. (2007), which adds support to their findings. We tendstoremovethesenewmutationsfromthepopulation.
alsofoundonenewdeletion,520_522delGTC,andonenew Thus,table2showsthatTajima’sDandFu’sandLi’sDandF
insertion, 539insC. However, for the subsequent diversity testsare highly significant even under modelsthat incorpo-
analyses, indel polymorphisms were not considered due to ratearecentOut-of-Africademographichistory(Gutenkunst
theirdifferentnaturetoSNPs.Themostfrequentvariantwas etal.2009; Lavaletal. 2010),whereasthe combinedDHEW
2655
MBE
Mart´ınez-Cadenasetal. . doi:10.1093/molbev/mst158
Table 1. Diversity Parameters for the Coding Region of MC1R in Several Populations.
Samplea 2Nb Lengthc Sd No. of Hapse hf pg h^ h ki Taj’s D (P Value)j
w
This work 3,142 953 42 42 0.517 6.3 4.867 0.599 (cid:2)2.143 ((cid:3)0)
TSI 196 953 10 11 0.562 6.9 1.709 0.659 (cid:2)1.437 (0.048)
ALL-AFR 492 950 11 12 0.184 2 1.624 0.192 (cid:2)1.913 (0.001)
ASW 122 950 10 11 0.325 3.7 1.860 0.353 (cid:2)2.014 (0.001)
LWK 194 950 4 5 0.071 0.8 0.685 0.072 (cid:2)1.560 (0.010)
YRI 176 950 5 6 0.203 2.2 0.870 0.212 (cid:2)1.477 (0.048)
ALL-NEU 536 951 10 11 0.674 9 1.458 0.854 (cid:2)0.870 (0.210)
CEU 174 951 8 9 0.671 9 1.395 0.852 (cid:2)0.873 (0.215)
FIN 186 951 7 8 0.663 8.7 1.207 0.827 (cid:2)0.676 (0.282)
GBR 176 951 8 9 0.682 11 1.393 0.863 (cid:2)0.852 (0.224)
D
ALL-ASI 572 953 6 8 0.525 9.1 0.866 0.913 0.097 (0.614) o
w
n
CHB 194 953 4 5 0.555 10.2 0.685 0.967 0.736 (0.801) lo
a
CHS 200 953 4 5 0.538 10.6 0.681 1.009 0.856 (0.828) d
e
JPT 178 953 4 6 0.435 7.3 0.695 0.664 0.081 (0.541) d fro
aSamples labeled “This work” correspond to healthy individuals from Spain. The rest correspond to the 1000G Project populations: the south-European (TSI) population m
(ptohoelesdo:uAthL-LE-uArFoRp)e,atnheIBnSoprtohp-uEulartoiopneafnrompop1u0l0a0tiGonwsa(sinndoivtidcuoanllsyid:eCrEeUd,dFuINe,toGBitRs;ltohwestahmrepeleposiozele;d2:NA=LL2-8N)E;Uth),eaAnfdrictahnepAospiaunlatpioonpsul(aitniodnivsid(uinadllyiv:idAuSaWlly,:LCWHKB,,YCRHI,S,anJPdT,alalntdhreaell https
threepooled:ALL-ASI). ://a
b2N=Numberofhaplotypes. ca
cNucleotidesconsideredafterremovingcodingsynonymouspolymorphisms.InthecaseofsouthEuropeansamples,weexcludedT314T,causedbyapolymorphisminnucleotide de
m
942;inthecaseofthenorthEuropeansamples,weexcludedI264I,T314T,andS316S,causedbypolymorphismsatnucleotides318,531,and942,respectively.IntheAfrican ic
samples,weexcludedL106L,T177T,F300F,andT314T,causedbypolymorphismsatnucleotides318,531,900,and942,respectively.IntheAsiansamples,weexcludedA103A, .o
Y298Y,andT314Tcausedbypolymorphismsatnucleotides309,894,and942,respectively. up
dNumberofsegregatingsites. .co
eNumberofdifferenthaplotypesobservedinthesample. m
fHaplotypediversity. /mb
gNucleotidediversity((cid:4)10(cid:2)4). e/a
hTheta(persequence)estimatedfromthenumberofsegregatingsites. rtic
iAveragenumberofnucleotidedifferencesbetweenhaplotypes. le
jTajima’sDvalue(significanceinferredafter5.000standardcoalescentsimulationswithDnaSP5.1). -ab
s
tra
c
t/3
0
/1
2
/2
6
5
4
/1
0
1
4
6
4
5
b
y
g
u
e
s
t o
n
1
2
A
p
ril 2
0
1
9
FIG.1. SlidingwindowgraphofpolymorphismandCpGvaluesaroundMC1R.Thegraphshowsthevaluesobtainedfordifferentparametersinsliding
windows(ofsize500bp,shifted100bp)acrosstheregion.Ontop,thecodingregion(broadgraybox)andUTRs(graylines)ofMC1R.Belowthis,the
continuouslineindicatesthattheGCcontentacrosstheentireregionishigh(~65%).However,theobservedtoexpectedCpGratio(continuousblack
linewithregularlyintersperseddots)substantiallyincreasesinthecodingregion(~0.75).AhighcontentofCpGs,ifmethylated,couldresultinahigh
number of new variants (see Theta W, continuous line with regularly interspersed asterisks), which because ofpurifyingselectionresult in highly
negativeTajima’sDvalues(approximately(cid:2)2;discontinuousgrayline).
2656
MBE
VaryingSelectiononMC1RAlleles . doi:10.1093/molbev/mst158
test,whichisspecificofrecentpositiveselection(duetothe east Asian populations; table 1), which is characterized by
incorporationoftheHtest),isnot.Thisisinagreementwith the high frequency (~0.64) of R163Q, an r allele (Tajima’s
the observation that many of the polymorphisms at this D=(cid:2)0.393, P value=0.411, under standard coalescent
region are nonsynonymous and that the most common simulations).
alleleisalsotheancestralone.
Using data from the 1000 Genomes (1000G) Project, we
Evidence for Positive Selection at V60L in South
observethatthisscenarioofpurifyingselectioninourSpanish
Europeans
population could be shared with other southern European
The strong pattern of purifying selection maintaining the
countries(table1).Thus,forinstance,Tajima’sDtestapplied
totheTSIsampleshowsaDvalueof(cid:2)1.437(Pvalue=0.0438, ancestralMC1R allelecontrastswitha substantialfrequency
ofV60Linoursample.However,therelativelyhighfrequency
understandardcoalescentsimulations).Asexpected,asimilar
of V60L is not homogeneously distributed within Spain, as
patternofpurifyingselectioncanbeobservedforthepooled
1000G African sample, for whom Tajima’s D is (cid:2)1.913 thereseemstobesomespatialheterogeneityinthefrequency D
o
(P value=0.001, under standard coalescent simulations), of this allele (table 3; supplementary information S2, w
n
and also for the individual African populations (table 1). SupplementaryMaterialonline).Toexplorewhetherthishet- loa
However, the pooled northern European sample (or the erogeneity in the frequency of the V60L allele in Spanish de
d
north European populations individually; table 1) does not populations could be explained by a general underlying fro
sharethispattern(Tajima’sD=(cid:2)0.870,Pvalue=0.21,under demographic process in the region, we then assessed the m h
standard coalescent simulations) The south European– existence of population stratification by means of principal ttp
s
African purifying selection pattern is not observed in the component (PC) analyses based on data derived from 93 ://a
pooledeastAsian sampleeither (orin anyof the individual European ancestry informative markers (termed EuroAIMs), ca
d
that are well known to recover the largest features of the e
m
north-northwest (NNW) to south-southeast (SSE) axis of ic
.o
Table 2. Neutrality Tests for the Observed MC1R Diversity in the genetic differentiation in Europe (Price et al. 2008). In fact, up
Spanish Sample. analyzing 1,187 individuals from different European popula- .co
m
Taj’s Da F&W’s F&L’s F&L’s DHEWe tions including the Spanish, EIGENSOFT found significant /m
nHb Dc Fd onlythefirstPC(Tracy–Widomtest,Pvalue=2.06(cid:4)10(cid:2)34), be
Observed values (cid:2)2.143 0.161 (cid:2)7.300 (cid:2)6.000 N.A.f differentiating NNW-Europeans from SSE-Europeans, and /artic
P values, under models leaving the Spanish individuals at intermediate positions in le
Single, constant- (cid:3)0.000* 0.396 0.000* 0.000* (cid:3)0.112 the NNW–SSE axis of European differentiation. However, -ab
s
size popg usingthissetofEuroAIMs,EIGENSOFTwasunabletoidentify tra
c
Guntekunst et al. (cid:3)0.0002* (cid:3)0.597 (cid:3)0.0001* (cid:3)0.0003* (cid:3)0.325 a significant PC (Tracy–Widom test, P value=0.646 for the t/3
(2009) firstPC)whenrestrictingtheanalysistotheSpanishpopula- 0/1
Laval et al. (cid:3)0.001* (cid:3)0.942 (cid:3)0.0007* (cid:3)0.002* (cid:3)0.992 tions(includingalltheSpanishprovinces,theSpanishBasque, 2/2
(2010) 6
and the Canary Islanders). To check whether this heteroge- 5
4
aTajima’sDtest. neity could be reflecting the geographical UV-B irradiation /10
bFay’sandWu’snormalizedHtest. pattern, we selected two population groups, one including 14
cFuandLi’sDtest. 64
dFuandLi’sFtest. thoseprovincesforwhichUV-Birradiationvalues(fig.2)lie 5 b
eZengetal.’scombinedDHEWtest. inthelower5thpercentileintheIberianPeninsula(group1: y g
fNotapplicable. Asturias,Bizkaia,Gipuzkoa;allofthemnorthernpopulations), ue
gNorecombinationassumed. ananotheroneincludingthoseprovincesabovethe95thper- st o
*SignificantafterBonferronicorrection. n
centileofUV-BintheIberianPeninsula(theCanaryIslands, 1
2
A
p
ril 2
Table 3. Fisher Exact Test of V60L (rs1805005) Allele Frequency. 0
1
9
Group 1 vs. Group 2a 2N Alleles Frequency V60L Allele Fisher Exact Test (P Value)b
Group 1 Group 2 Group 1 Group 2
a) UV-B(cid:5)5th vs. UV-B(cid:6)95th 188 134 32 13 0.042
b) UV-B(cid:5)10th vs. UV-B(cid:6)90th 266 170 41 16 0.046
c) UV-B(cid:5)15th vs. UV-B(cid:6)85th 348 270 55 28 0.057
d) UV-B(cid:5)20th vs. UV-B(cid:6)80th 466 314 72 37 0.088
aPopulationgroupscompared:ina)group1includesthoseprovinceswhosedaily(annualmean)receivedUV-Bradiationislowerthanorequaltothe5thpercentileofthe
distributionofUV-BradiationintheIberianPeninsula.Group1:Asturias(1),Bizkaia(2),andGipuzkoa(3)(allofthemnorthernpopulations);group2includesthoseprovinces
whosedaily(annualmean)receivedUV-Bradiationisgreaterthanorequaltothe95thpercentileofthedistributionofUV-BradiationintheIberianPeninsula:theCanary
Islands(12),Ca´diz(13),andGranada(14).Similarlyforb),c),andd)butconsideringthepercentilesindicated.Inb)group1,asina)plusCantabria(4),Araba(5),andLugo(6);
group2asina)plusAlmer´ıa(15),Huelva(16),andCa´ceres(17).Inc)Group1,asinb)plusACorun˜a(7)andNavarra(8);Group2asinb)plusSevilla(18)andJae´n(19).Ind)
Group1asinc)plusPontevedra(9),Ourense(10),andLaRioja(11);Group2asinc)plusCiudadReal(17),Ma´laga(20),andMelilla(21).Onlyindividualswhosegrandparents
werealsobornintheseprovinceswereconsidered.Seefigure2tolocalizeprovincesorregionsnumericcodesonthemap.
bFisher’sexacttestone-tailedPvalues.
2657
MBE
Mart´ınez-Cadenasetal. . doi:10.1093/molbev/mst158
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
FIG.2. UV-BradiationmapofSpain.MapshowingtheannualaveragedailyUV-Bradiation(J/m2)perprovince(seeMaterialsandMethods).Note ://ac
a
thatforaestheticalpurposestheCanaryIslandsaredisplayedclosertotheIberianPeninsulathantheyactuallyare(i.e.,offthecostofthesouthern d
e
endofMorocco).Dotscorrespondtothecentroidsofthepolygonsrepresentingeachprovinceandnotthelocationofthecapitalofeachprovince. m
ic
MaprepresentationandcolorinterpolationhasbeencarriedoutwithQGISv1.8.0.Numericcodesforprovincesorregionsindicatedintables3and .o
u
4 are as follows: Provinces with the lowest UV values: Asturias (1), Bizkaia (2), Gipuzkoa (3), Cantabria (4), Araba (5), Lugo (6), A Corun˜a (7), p
.c
Navarra (8), Pontevedra (9), Ourense (10), and La Rioja (11); provinces with the highest UV values: the Canary Islands (12), Ca´diz (13), Granada o
m
(14),Almer´ıa(15),Huelva(16),Ca´ceres(17),Jae´n(18),Sevilla(19),Ma´laga(20),Melilla(21),andCiudadReal(22);andmelanomasamplesweretaken /m
fromCastello´n(23). be
/a
rtic
le
-a
Ca´diz, and Granada; all of them southern populations) plus informationS3,SupplementaryMaterialonlineforEHHanal- b
s
(group2),Fisher’sexacttestshowsthatthefrequencyofV60L ysisinsinglepopulations).Weperformedanadditionaltestof tra
c
ingroup1(0.17)issignificantlyhigherthanthatingroup2 selection based on a haplogroup-specific Tajima’s D test on t/3
0
(0.09) (one-tailed P value=0.042). This is confirmed if we a region of 10kb 50 of rs1805005 for the set of European /1
2
consider other percentiles (table 3) also. Nevertheless, the rs1805005-T (V60L) haplotypes (see Materials and /2
6
distribution of the V60L allele frequency is far from closely Methods). Tajima’s D for the rs1805005-T haplogroup in 54
following a strictly linear latitudinal pattern. Thus, when we this population was (cid:2)2.045, whereas for the alternative /10
1
groupsamplesaccordingtotheirUV-radiationpercentile(at rs1805005-G haplogroup Tajima’s D was 1.693. The P value 4
6
4
noninclusive 0.05 quantiles from the lower and upper ends, under the Gutenkunst et al. (2009) demographic model for 5
b
respectively),thereisactuallyanegativecorrelationbetween theThaplogroupis0.021(0.032undertheLavaletal.[2010] y
g
theV60LallelefrequencyandUVradiation(Pearson’scorre- model). This indicates that positive selection can be held ue
s
lationcoefficient=(cid:2)0.165),butthisisnotsignificant(ttestP responsible for the substantial frequency of V60L at least in t o
n
value=0.25).Thissuggeststhatotherevolutionaryprocesses, someareasofEurope. 1
2
including demographic history, may interfere. However, we A
The Collateral Damage of Positive Selection p
alsohavetotakeintoaccountthelowsamplesizesformany ril 2
of the bins, which result in imprecise estimates of allele fre- Itisinterestingtonotethatindividualswithmelanomashow 0
1
quencies (supplementary information S2, Supplementary a significantly higher frequency of V60L than their corre- 9
Materialonline). sponding controls (Fernandez et al. 2007; Ibarrola-Villava
Notwithstanding, the relative high frequency of V60L in et al. 2012; see also http://www.melgene.org/, last accessed
theSpanishpopulationcouldbeindicativeofsomeselective October3,2013).Wehaveimprovedthisdatasetwithnew
advantage. Therefore, we then investigated whether recent melanomacasesandcompareditwithourSpanishreference
positive selection could explain this pattern. For that, we setasthecontrolsample.Cochran–Armitagetestalsoshows
applied the extended haplotype homozygosity (EHH) test anassociationbetweenV60Landmelanomariskinanaddi-
tothe1000GsouthernEuropeansample,thebestproxyavail- tivemodel(table4)onlywhenconsideringascontrolstheset
able,astheTSIpopulationhassimilarlatitudetonorthSpain. ofpopulationsbelongingtotheprovinceswithhighestUV-B
Figure3showsthattheextentofhomozygosityfromV60Lis irradiation (above the 95th, 90th, 85th, or 80th percentiles).
infactlongerthanexpectedundertheneutraldemographic Thissuggeststhattheadaptiveevolutionofskindepigmen-
models considered. It is also longer than that observed for tationisconcurrentlyresponsibleforanincreasedhealthrisk
the 1000G north European sample set (see supplementary in,mosttypically,postreproductiveage.
2658
MBE
VaryingSelectiononMC1RAlleles . doi:10.1093/molbev/mst158
D
o
w
n
lo
a
d
e
d
fro
m
h
Fig.3. EHHanalysisofthe1000GEuropeansamples.ObservedandsimulatedEHHvaluesfortheextendedV60Lhaplotype.(A)simulationsusingthe ttp
s
demographicmodelforEuropeansofLavaletal.(2010)and(B)usingthedemographicmodelforEuropeansofGutenkunstetal.(2009).Shadesofblue ://a
correspondtoasmoothedcolordensityrepresentationofthescatterplotofallEHHsimulationsforeachregion,obtainedthroughakerneldensity c
a
d
estimate, usingthe smoothScatter function in R. Blue darkness isproportional to the density of pointsin that vicinity. The continuous blue line e
m
correspondstothe95thpercentileofthedistribution(seeMaterialsandMethods),smoothedusingtheloessfunctioninR.Orangedotsarethe ic
observedvaluesfromthe1000GdataforsouthEuropeans.Forcomparativepurposes,blackdotscorrespondingtotheEHHobservedprofileof1000G .ou
p
datafromnorthEuropeans(GBR,FIN,andCEU)arealsoshown.EvidenceforselectioninnorthEuropeansisnotsoevidentasinsouthEuropeans,and .c
o
inanycase,muchweaker. m
/m
b
e
/a
rtic
Table 4. Cochran–Armitage Association Test for V60L (rs1805005) and Melanoma Risk. le
-a
Sample N (inds) Genotypesa Exact HW Test (P Value) Cochran–Armitage bs
Trend Test (P Value)b tra
GG GT TT c
Melanomac 722 494 210 18 0.54 t/30
UV-B (upper 95%)d 72 52 8 1 0.66 0.0134 /12
/2
UV-B (upper 90%)e 85 62 11 1 0.87 0.0089 6
5
UV-B (upper 85%)f 139 101 21 2 0.70 0.0059 4/1
0
UV-B (upper 80%)g 157 115 28 3 0.60 0.0222 1
4
6
aTheTalleleisresponsiblefortheV60Laminoacidchange. 45
bTestedagainstthemelanomasampleunderanadditivemodel. by
cIndividualsfromourownsequencedmelanomasamplesfromCastello´n(22)plusthoseofFerna´ndezetal.(2007)andIbarrola-Villavaetal.(2012). g
u
dIncludesindividualsfromtheprovincesoftheCanaryIslands(12),Ca´diz(13),andGranada(14)andwhosegrandparentsarealsofromthesameprovince. es
eAsinthefootnote“d” plusAlmer´ıa(15),Huelva(16),andCa´ceres(17). t o
n
fAsinthefootnote“e”plusSevilla(18)andJae´n(19). 1
gAsinthefootnote“f”plusMa´laga(20),Melilla(21),andCiudadReal(22).Seefigure2forlocalizingtheprovincesorregionsonthemapbyusingthecorrespondingnumeric 2 A
codesbetweenparentheses. p
ril 2
0
1
On the other hand, there were no significant differences observations add support to the hypothesis that: 1) CpGs 9
whenusingascontrolseitherthesetofpopulationswiththe withingenes,contrarytoregulatoryCpGislands,tendtobe
lowestUV-BradiationorthewholesetofSpanishprovinces. highlymethylated(Jones1999);2)methylomeanalysisofH1
This is illustrative of how population structure in allele fre- humanembryonicstemcellsandIMR90fetallungfibroblasts
quencies can be an issue in genetic association tests, which shows that MC1R is highly methylated (Lister et al. 2009)
canleadtofalsenegativesifnotproperlyaccountedfor. (http://neomorph.salk.edu/human_methylome/browser.
html, last accessed October 3, 2013); 3) our own bisulfite
Discussion
sequencing data indicate that MC1R is methylated in
We have shown that MC1R shows high coding nucleotide sperm; 4) in mammalian genomes, CpG sites are hotspots
diversity in the (south European) population analyzed. for transition mutations because of methylation-mediated
Wehaveproposedthatmethylation-mediateddeamination deamination of 5-methylcytosine (Cooper and Krawczak
of5-methylcytosineisultimatelyresponsibleforthehighnu- 1993); 5) 88% (44 out of 50) of the SNPs observed in our
cleotide diversity observed at MC1R. The following sample (considering both controls and melanomas) are
2659
MBE
Mart´ınez-Cadenasetal. . doi:10.1093/molbev/mst158
transitions;and6)CpGsitesshowamutationrate10times On theotherhand,despite proving evidence thatdepig-
higherthannon-CpGs(NachmanandCrowell2000). mentationcanbeadaptiveinrelationtoscarcesunlight,we
This high variability is, however, shaped by contrasting however,cannotexplicitlysaywhichistheunderlyingphys-
evolutionary forces. It has been argued that environmental iological mechanism involved. The vitamin D hypothesis
pressure is responsible for maintaining the ancestral MC1R (Loomis 1967) would constitute a straightforward connec-
alleleathighfrequencyinAfrica,whereasinEurasiainstead, tion.Nevertheless,thereisstillagapthatneedstobefurther
thisfunctionalconstraintisabsent(Ranaetal.1999;Harding exploredinorderfortheearliermentionedhypothesistobe
etal.2000).However,giventheconspicuouspatternofpuri- proven.
fyingselectiondetectedinoursample,weshowthatstrong Finally,italsoremainstobeknownwhyotherR/rmuta-
functional constraint on the ancestral MC1R haplotype also tionsarenotsocommoninourpopulation.Onehypothesis
takes place at least in some areas of southern Europe. This isthatrmutationsresultinalighterphenotypebutinaway
observationmightsubsequentlyleadustothinkthatMC1R thatdoesnotcompromisetanningduringseasonalincreases
has played no role in the adaptive lightening of the human in UV-B, and thus enabling protection against skin cancer Do
w
skininEurasia.Actually,MC1Rvariabilityhasbeenfoundnot (Robins 1991) or against folate destruction (Jablonski and n
lo
tocontributetolightpigmentationinmicefromtheFlorida’s Chaplin2010).Inthisregard,V60Lhasbeenshowntoachieve a
d
Atlantic Coast (Steiner et al. 2009), although it had been significantlylowercAMPlevelsafteragoniststimulationthan ed
highlighted as an adaptive locus in beach mice from thewildtypeMC1Rprotein,whichindicatessubstantialim- fro
m
Florida’sGulfCoast(Hoekstraetal.2006).Similarly,innonhu- pairedfunction,similartothatoftheRvariants(Schio¨thetal. h
manprimates,MC1Rvariabilitydoesnotseemeithertoshow 1999;Beaumontetal.2007,2009).Thisresultissupportedby ttps
a simple relationship to coat color (Mundy and Kelly 2003; geneticdata,astheskincolorofindividualshomozygousfor ://a
c
Haitina et al. 2007). Besides, the presence of a functional, V60L is significantly fairer than wild-type individuals, to the ad
e
MSH-responsive MC1R in orangutan shows that alternative same extent, for instance, as homozygotes for the R allele m
mechanismsofcoatcolorgenerationcanexist(Haitinaetal. R160W (Beaumont et al. 2007). However, in heterodimers ic.o
u
2007). with the wild-type receptor, V60L does not seem to show p.c
However, our results indicate that two selective forces (butmostRmutantsdo)dominant-negativeeffectonwild- om
shape human MC1R diversity in Sothern Europe. Similarly type receptor’scell-surface expression.A dominant-negative /m
b
effect would result in a decreased ability of the wild-type e
toAfrica,purifyingselectionisstrongmaintainingtheances- /a
tral allele in Sothern Europe. But simultaneously, positive receptor (the ancestral allele) to elevate intracellular cAMP rtic
levels (Beaumont et al. 2007, 2009), and would add to the le
selection is responsible for an increased frequency of the -a
V60L allele, a low-penetrance “red hair and fair skin” allele decreased ability of the V60L allele itself. Thus, generally bs
speaking,V60Lwouldpreservetheabilitytoincreaseeume- tra
(alsoassociatedtoblondehair;Boxetal.1997).Thefactthat c
laninproduction(tanning)duringseasonalincreasesinUV-B t/3
V60L is under positive selection constitutes evidence that 0
betterthantheRalleles.Thiswouldbeusefulparticularlyin /1
MC1R is also playing an active adaptive role in human skin 2
depigmentation, at least for some areas of south Europe. southernEuropeanregions,wheresummerUV-Bradiationis /26
5
higher than in northern Europe. Thus, in conclusion, our 4
Actually, the world frequency distribution of V60L shows /1
resultswouldsupportthe hypothesisproposedbyJablonski 0
thatthisalleleiscircumscribedtoEuropeandtheneareast 1
4
and Chaplin (2010) that “habitation of middle latitudes 6
(Gerstenblithetal.2007;Galore-Haskeletal.2009;Eskandani 4
et al. 2010; Acar et al. 2012). However, we show that this between approximately 23(cid:7) and 46(cid:7) involved the evolution 5 b
y
increaseinfitnessforthepopulationhascollateraldamaging ofpartiallydepigmentedphenotypescapableoftanning.” gu
e
consequences for the individual’s health. In this regard, the Materials and Methods st o
samevariant(V60L)thatisunderpositiveselectionisalsoa n
1
melanoma-risk allele. As melanoma is typically a postrepro- Population Samples Sequenced 2 A
daruectpivoestdpiosenaisneg(Stha´eenpzaeytmael.n2t00d5u)e,wtoe coaunrcaodnaspidtaetriotnhattowae Wfroemadniafflyezreendt3r,e1g4i2onhsuomfSapnaiMnCth1Rroualglehledsir(e1c,t57D1NiAndsievqiduueanlcs)- pril 20
1
reduced-sunlight environment. Does this mean that mela- ingof thecodingregion.Ofthe1,571individuals,962came 9
noma,ormoregenerally,skincancerriskhasnofitnesscon- fromtheSpanishNationalDNABank(http://www.bancoadn.
sequencesitself?Wecannotanswerthatquestion withour org/en/home.htm,lastaccessedOctober3,2013)andrepre-
data, because although the ancestral MC1R allele is under sentdifferentregionsinSpain.Ofthese,670individualswere
strong purifying selection (not only in our sample but also known tohaveall grandparentsfrom asinglespecific prov-
inmoredemandingregionslikeAfrica),wewouldalsoneed inceofSpain.Inaddition,wesequencedafurthersampleof
toassessotherprotectiverolesthathavebeenproposedfor 172allelesfromtheBasqueCountry(86individuals),60alleles
darkskin,mostnotablyphotoprotectionofcirculatingfolate (30 individuals) from the Canary Islands (a subset of those
(JablonskiandChaplin 2010).However, the role of MC1R in described in Pino-Yanes et al. [2011], selected based on the
DNArepairandapoptosissuggeststhatMC1Risnotsimply lower ancestry to northern Africa as deduced from 93
encodingapassivefilter(melanin)(Garc´ıa-Borro´netal.2005). EuroAIMs), and 278 alleles (139 individuals) from the prov-
Inanycase,whateverthefitnessforcesatplay,thefinalfitness inceofCastello´n.SequenceshavebeendepositedinGenBank
balanceseemstofavordepigmentationinEurope. with accession numbers KC981760–KC984193. Finally, 708
2660
MBE
VaryingSelectiononMC1RAlleles . doi:10.1093/molbev/mst158
already sequenced MC1R alleles (354 individuals) from MC1R-R3 50-AAGGGTCCGCGCTTCAACACTTTCAGAG-30,
Madrid were recruited from previously published work MC1R-R4 50-CAACACTTTCAGAGATCATTTAGTCCATC-30,
(Fernandezetal.2007;Ibarrola-Villavaetal.2012). andMC1R-R550-CAGCACGGCCATGAGCACCAG-30.
Sequencingreactionswereperformedusing10ngofPCR
Data Compiled from Other Sources purified product (MicroSpin S-400 HR Columns, GE
Healthcare) and BigDye Terminator Cycle Sequencing kits
We obtained full MC1R sequence information as strings of
(Life Technologies). Sequencing reaction products were
SNP genotypes from the 1000G project data (http://www.
then washed with Performa DTR V3 96-well plates
1000genomes.org/, last accessed October 3, 2013) by means
(EdgeBio), and they were finally run on an ABI 3130
of SPSmart v.5.1.1 (http://spsmart.cesga.es/, last accessed
Genetic Analyzer (Life Technologies). Sequences were
October 3, 2013) (Amigo et al. 2008), which uses the
editedwithGenalysv2.8(Takahashietal.2003)andfurther
information 1000 Genomes Phase I data (data freeze as of
checkedmanuallyafterward.
May2011)forthefollowingpopulations:
D
o
European populations: 268 individuals from northern Bisulfite Methylation Analysis wn
Europe, comprising Utah residents with northern and Methylationwasassessedusingabisulfatetreatmentproto- load
western European ancestry from the CEPH collection e
colbasedonKanedaetal.(2004).Inbrief,3mgofDNAfrom d
(CEU, N=87), Finnish from Finland (FIN, N=93), and spermwasdigestedwithHinfIandincubatedwith6NNaOH fro
m
British from England and Scotland (GBR, N=88); and for15minat37(cid:7)C.Tothissolution,weadded120mLofthe h
112 individuals from southern Europe, comprising conversion solution: sodium bisulfite 4.04M (Sigma), ttps
Tuscans from Italy (TSI, N=98) and Iberian Spanish Hydroquinone 10mM (Sigma), and NaOH 6 N (pH 5), and ://a
(IBS,N=14). performed20cyclesofdenaturingfor30sat95(cid:7)Candincu- cad
African populations: 246 individuals of African ancestry bationfor15minat50(cid:7)C.DNAwasthenpurifiedusingthe em
from southwest USA (ASW, N=61), Luhya in Webuye, ic
QIAquick PCR purification kit (Qiagen), desulfonated by in- .o
NKe=ny8a8)(.LWK,N=97),andYorubainIbadan,Nigeria(YRI, cubationwithNaOH6Nfor5minatroomtemperature,and up.co
precipitatedwithethanol100%andammoniumacetate6M m
East Asian populations: 268 individuals from China (pH7)(Sigma). /m
(CHB), Han Chinese South (CHS), and Japanese in ConvertedDNAwasamplifiedbyPCRwiththefollowing be/a
Tokyo,Japan(JPT). primers:Fw50-TGTTGGTTCGGGTTTGTTA-30 andRev50-AC rtic
CTCCTTAAACGTCCTACG-30,whichweredesignedtoanneal le-a
b
Resequencing smpeetchifiyclaatliloyntointChpeGb)i.sPuClfiRtec-ocnondviteiortnesdwseerqeuaesnfcoello(wcosn:3simdeirninagt strac
All DNA samples were analyzed for the coding sequence of 96(cid:7)C; 45 cycles of 30s at 96(cid:7)C, 20s at 59(cid:7)C, and 1min at t/3
0
the MC1R gene by direct automated DNA sequencing. 70(cid:7)C; and a final step of 10min at 70(cid:7)C. Sequencing reac- /12
Polymerase chain reaction (PCR) primer sequences used to /2
tions were performed using the BigDye Terminator kit 3.1 6
framea1,150-bpfragmentcontainingtheentiresequenceof 5
theMC1Rgenesingleexonwereasfollows:50-CAGCACCAT (pLriifmeeTrescthhnaotlwogeirees)uasenddfrournthoenPaCnR.ABI 310 using the same 4/101
GAACTAAGCAGGACACCT-30 (sense)and50-AAGGGT 4
6
CCG CGC TTC AAC ACT TTC AGA G-30 (antisense). PCR 45
Probabilistic and Experimental Phasing b
reactionswerecarriedoutusing100ngofDNAastemplate, y
andthePCRreactionmixturealsoincluded2mlof10(cid:4)PCR Haplotypes were estimated with PHASE (available from: gu
e
buffer,25mMMgCl2,5mMeachofdGTP,dATP,dCTP,and http://stephenslab.uchicago.edu/software.html#phase, last st o
accessed October 3, 2013). For those pairs of SNPs whose n
dTTP,25pmolofeachPCRprimer,and1unitofAmpliTaq 1
phasecouldnotbestatisticallyinferredwithaPvalue>0.9, 2
Gold DNA polymerase (Life Technologies). Thermocycling A
fcoolnlodwiteiodnbsyi3n5clcuydceledsaofd9e5n(cid:7)aCtufroart3io0ns,s5t6ep(cid:7)Catfo9r54(cid:7)5Cs,faonrd87m2(cid:7)inC, w(TeOpPrOoc-TeeAdekdit,tIonvinitfreorgtehne)iranphdarseeseeqxupeenricminegntaaglalyinb.ycloning pril 20
1
for 1 min, plus a final extension cycle of 72(cid:7)C for 10min. 9
Amplified products were purified using ExoSAP-IT (USB Population Genetics Analysis of the MC1R Sequences
Corporation) according to the manufacturer’s instructions. Population diversity and neutrality tests parameters were
In order to obtain both forward and reverse sequences of obtained using DnaSP v5.10 (Librado and Rozas 2009)
theentireMC1Rexon,purifiedPCRproductsweresubjected (http://www.ub.edu/dnasp, last accessed October 3, 2013).
to four sequencing reactions using forward and reverse For the neutrality tests, P values were initially obtained
sequencing primers: MC1R-F1 50-CAG CAC CAT GAA CTA from simulations run under a neutral, constant-sized popu-
AGC AGG ACA CCT-30; MC1R-F2 50-TGG GTG GCC AGT lationmodelusingtheprogramms(Hudson2002).Twenty
GTCGTCTTCAGCA-30;MC1R-R150-AAGGGTCCGCGC thousand simulations were obtained fixing on (cid:3)^ . To take
w
TTC AAC ACT TTC AGAG-30; and MC1R-R2 50-CCA GCA into consideration the possible effect of demography on
TAG CCA GGA AGA AGA CCA CGA-30. Gaps were filled theneutralitytests,wecarriedoutadditionalsimulationsin-
using primers, MC1R-F3-50-ATGCCAGGAGGTGTCTGG corporatingoptimizedparametersthatattempttoreflectthe
ACTG-30, MC1R-F4-50-GATCACCTGCAGCTCCATGCTG-30, humanevolutionaryhistory.Theseparameterswereobtained
2661
MBE
Mart´ınez-Cadenasetal. . doi:10.1093/molbev/mst158
fromthemodelsofGutenkunstetal.(2009)andthatofLaval -Samples from 553 individuals from 14 Spanish prov-
etal.(2010).ThemscommandlinesfortheGutenkunstetal. inces from the Spanish National DNA Bank (www.
(2009)modelwereasfollows: bancoadn.org, last accessed October 3, 2013): Andalusia
(n=105),Murcia(n=27),Extremadura(n=25),Castile-
./ms314220000-t1.018-I3031420-n11.68202
La Mancha (n=50), Valencia (n=38), Castile and Leon
-n 2 3.73683 -n 3 7.29205 -g 2 116.010723 -g 3
(n=93), Madrid (n=6), Catalonia (n=55), Galicia
160.246047 -ma x 0.881098 0.561966 0.881098 x (n=69), Cantabria (n=18), Navarre (n=16), La Rioja
2.79746 0.561966 2.79746 x -ej 0.028985 3 2 -en (n=21),Asturias(n=14),andAragon(n=16).
0.028985 2 0.287184 -ema 0.028985 3 x 7.29314 x -Samples from 24 Spanish Basques satisfying that all
7.29314xxxxx-ej0.19796321-en0.30350111 theirfourparentalsurnameswereofBasqueorigin.
ThemscommandlinesfortheLavaletal.(2010)parameters In addition, data from 80 Canary Islanders taken from
wereasfollows: previous studies (Pino-Yanes et al. 2011), selected based on
D
thelowerancestrytonorthernAfrica,wereusedforcompar- o
w
./ms314220000-t1.922-I3031420-n11666.667 isons,alongwithdatafrom144Swedish,52Polish,69English, nlo
-n 2 2.261 -n 3 1.051 -g 1 372.278 -g 2 380.217 a
117Italians,64Greeks,and28FrenchfromPriceetal.(2008), d
e
-0g.713726509..7316766-mxa-ejx00.0.711673630.271-7en6 00..70147365 2x 00..27012796 aansdwewlleasstferronmE5u6roUpteah(rCeEsiUde),ntasllwoitfhtahnecmestraylsforomsatnisofyritnhgerna d from
-ej0.043521 genotypingcompletionratemorethan90%. http
s
Using these simulations, the P values for the Tajima’s D, ://a
EHH Analysis of the 1000G Data c
Fay’sandWu’sH(normalized),andFu’sandLi’sDandFtests a
d
wereobtainedusingthesoftwaremsstats(K.Thornton’slab, To assess the possibility of positive selection acting on em
http://www.molpopgen.org/software/msstats/, last accessed single alleles, we used our own implementation of the EHH ic.o
October3,2013).TheDHEWcombinedtest,whichincludes method(Sabetietal.2002),whichusessimulationsobtained up
Tajima’sD,Fay’sandWu’sH(normalized),andtheEwensand by means of Hudson’s ms program (Hudson 2002) to .co
m
Watersontestwereperformedwithsoftwarekindlyprovided estimatethePvalues.Thisprocesswasimplementedinthe /m
byK.Zengandmodifiedbyustoallowforlargersamplesizes. followingway: be
This compound test is more robust against the presence of We selected two 100kb regions, one on each side of the /artic
recombination, is insensitive to background selection and SNPwewantedtotestforrecentpositiveselection,and,for le
demography,anditisthusmorepowerfulindetectingpos- each region, we estimated divergence to the orthologous -ab
s
itiveselection(Zengetal.2007). region in the chimpanzee. From divergence (assuming a tra
c
The genealogical relationships among MC1R haplotypes 6 My divergence time and the Jukes and Cantor model; Li t/3
0
were inferred by means of Network v4.6 (http://www. 1997), we estimated the mutation rate (cid:4), which in its turn /1
2
fluxus-engineering.com, last accessed October 3, 2013). We was used to estimate the parameter (cid:3) for both the /2
6
used genetree v9.0 (by R.C. Griffiths, available from: http:// Gutenkunstetal.(2009)andLavaletal.(2010)models.The 54
www.stats.ox.ac.uk/~griff/software.html, last accessed recombinationparameter(cid:5)wasestimatedfromthegenome- /10
1
October3,2013)toinfertherootofthetree. wide recombination rate from Phase 2 HapMap estimated 46
4
from phased haplotypes in HapMap Release 22 (NCBI 36) 5
b
(International HapMap Consortium 2007). In each simula- y
Structure Analysis tion, the specific (cid:3) and (cid:5) values were randomly sampled gu
e
s
ToexplorethegeneticstructureintheSpanishpopulations, from a normal distribution with a standard deviation equal t o
datafrom93EuroAIMsrecoveringthelargestfeaturesofthe tothemean.Thedemographicmodelsusedforthesimula- n 1
2
European NNW to SSE axis of differentiation (Price et al. tionswerethesameasmentionedinGutenkunstetal.(2009) A
p
2as0s0e8s)s,afrPoCmanaalytositsablyomf e1a,1n8s7ofinEdIGivEidNuSaOlsF,Tw4a.2s (uPtriilcizeedettaol. aonfdsimLauvlaalteiotnasl.w(2a0s1a0s)s.uFmoreedactohrreepgiroense,nthtethfiersStNSNPPofininttheeressett. ril 20
1
2006).Datafinallyincludedinthisanalysiswerederivedfrom Only those simulations for which the frequency pattern 9
individuals known to have all grandparents from a single of this SNP, taking into account the ancestral or derived
specific province of Spain, and from reference populations, state, was equal (±10% allowance) to the observed values
satisfying a genotyping completion rate more than 90%. in the 1000G European population, were considered for
We restricted the genotyping and analyses to those 93 Pvalueestimations.Inaddition,onlythosesimulationsshow-
SNPs that overlapped with samples used as reference (see ingatotalnumberofSNPsper100kbequaltotheobserved
Pino-Yanes et al. [2011] for a full description of methods). value(±10%allowance)wereconsideredforEHHestimations.
Briefly, genotyping was conducted utilizing the iPLEX Gold EHH profiles from the set of simulated haplotypes carrying
assay on MassARRAY system (Sequenom, San Diego, CA) the (derived) allele under consideration (at the first SNP)
by the Spanish National Genotyping Center, Santiago de were then used to infer P values. Typically, for each region
Compostela Node (CeGen, http://www.cegen.org, last severaltensofthousandssimulationswereneededtoobtaina
accessed October 3, 2013), and was performed in the set of 400–500 simulations satisfying these conditions.
following: After this, the EHH values at points distant from the first
2662
MBE
VaryingSelectiononMC1RAlleles . doi:10.1093/molbev/mst158
SNP0.1,0.2,0.3,0.4,and0.5units(inascaleof0to1,being1 parameter) option is a random number from a
equivalentto100kb)wererecordedandtheindividualsim- normal distribution with mean 3.25 and variance
ulationsshowingtheEHHvalueatthe95thpercentileateach equaltothemean.
pointwereselected.Allsimulationselectedwerelumpedto- FortheLavaletal.(2010)model:
gether.Next,duplicatedpointswere removed,andthen we
./ms760252-ttbs-rtbs10000-I307600-n1
forced all points to satisfy the condition EHH (cid:5)EHH
i+1 i
1666.667-n22.261-n31.051-g1372.278-g2
(where i increases from SNP number 1 in any direction).
The resulting distribution of EHH values was used as the 380.217-g3259.366-max0.71760.71760.7176
95%cut-offline. x0.71760.71760.7176x-eg0.006320-eg0.0063
30-ej0.016332-eg0.019910-ej0.043521-en
Tajima’s D Test for Specific Haplogroups
0.043520.2029
AsanadditionaltestforpositiveselectiononV60L,wecarried
wherethe-t(theta)optionisarandomnumberfrom
outTajima’sDtestonthatsubsetofhaplotypesdefinedby D
thederivedstateTatrs1805005(V60L),usingtheinformation a normal distribution with mean 18.4 and variance ow
n
fromthe1000Gproject.Weinitiallytriedtocarryoutthistest equaltothemean,andthe-r(recombinationparam- lo
a
eter) option isa random number from a normaldis- d
onthesouthEuropean1000Gsamples,butwithinthehap- e
d
logroup defined by allele rs1805005-T, there is no variability tribution with mean 6.14 and variance equal to the fro
atleast10kbonbothsidesofrs1805005.Thus,wedecidedto mean. m
testthewhole1000GEuropeansample(GBR,CEU,FIN,TSI, b)Then,foreachsimulation,wecheckthatatacertain http
ainndmIoBrSepeoxotrleemd)e.NnoetgeattihvaetTitaijsimuan’lsikDelyvatlhuaets,psoootlihnigsiwsiallprerisourlit, panrecdeeofifne5d%p)otsoititohnosaelleolebsferervqeudenactierss1a8re05c0lo0s5ea(atlltohwat- s://aca
a conservative approach. Besides, there is some evidence population. This predefined position is normally the de
m
(althoughweaker)fromtheEHHtestthattheV60Lalleleis firstorthelastSNP,becausethisfacilitatesstepsb),c), ic
also under some positive selection in the 1000G north andd).Similarly,wecheckthatthenumberofsegre- .ou
p
European populations (fig. 3). Note that evidence for selec- gating sites is close to the observed value (allowance .co
tion on the pooled 1000G European population is only of 10%) for the whole sample. If these conditions m/m
apparent when we analyze the internal diversity of the are satisfied, we keep that simulation, otherwise we be
rgslo18b0a5l0d0iv5e-Trsithya,pthloagtroisu,pincalnuddinngobtowthhersn18w05e00a5n-aTlyaznedt-hGe c) dTihsceanr,dfoitr. each accepted simulation, we extract the /article
haplogroups. set of haplotypes that are of a desired allelic state at -ab
s
Forthishaplogroup-specificTajima’sDtest,wedecidedto the predefined position (0=ancestral, 1=derived), tra
analyze a region comprising 10kb of the upstream flanking emulating the rs1805005-T (derived) allele in that ct/3
region to V60L. Using just one flanking region of V60L in population. 0/1
which one SNP at the beginning or the end is the V60L, d)Wecheckthatthenumberofsegregatingsiteswithin 2/2
SNPismoreefficientwhencalculatingthePvalues(discussed this extracted set of haplotypes is similar (allowance 65
4
later). Even after pooling 1000G data from northandsouth 10%) to that observed for that haplogroup in the /10
Europeansamples,thereisjustonepolymorphismwithinthe 1000G population. For 10kb upstream, rs1905005-T 14
6
Thaplogroupinthedownstream10kbregiontoV60L,which alleles6segregatingsitesareobserved,sotheallowed 45
doesnotgiveuspowertoperformTajima’sDtest.However, rangewasbetween5and7segregatingsites.Notethat by
g
there are6 SNPsin the upstream10kbregionwithin theT severalmillionsimulationsareneededtoendupwith u
e
haplogroup. 100–200 simulations satisfying our conditions. For st o
ToestimatethePvaluesforthers1805005-Thaplogroup, eachacceptedsimulation,wecalculateitscorrespond- n 1
weproceededasfollows: ingTajima’sDvalue(forthespecifiedsetofhaplotypes 2 A
a)First,runcoalescentsimulationsusingparametersfrom fToarjmiminag’saDhvaapluloegsrfoourpm)s. Tohuer dseisttroibfutthioensofr-ocmalcuwlahtiecdh pril 20
thewholeEuropeanpopulation 1
to infer the P value for the observed Tajima’s D of 9
FortheGutenkunstetal.(2009)model:
the rs1805005-T haplogroup of the 1000G European
./ms760690-ttbs-rtbs10000-I307600-n1 sample.
1.68202 -n 2 3.73683 -n 3 7.29205 -g 2
116.010723 -g 3 160.246047 -ma x 0.881098
Association Analysis
0.561966 0.881098 x 2.79746 0.561966 2.79746
To detect the association of some MC1R variants to
x -ej 0.028985 3 2 -en 0.028985 2 0.287184
melanoma susceptibility, we sequenced MC1R from 127
-ema 0.028985 3 x 7.29314 x 7.29314 x x x x x
DNAsamplesfrommelanomapatientsattendingtheoutpa-
-ej0.19796321-en0.30350111
tient clinic at the Castello´n Province Hospital (Castello´n,
where the -t (theta) option is a random number Spain) for a 5-year period, between 2007 and 2011. In addi-
from a normal distribution with mean 9.75 and vari- tion, already published MC1R sequence data on 595 mela-
ance equal to the mean, and the -r (recombination noma patients from Madrid were also recruited (Fernandez
2663