evidence of adaptation from ancestral variation in young populations of beach mice

15 Pages·2012·1.01 MB·English

ORIGINAL ARTICLE
doi:10.1111/j.1558-5646.2012.01669.x
EVIDENCE OF ADAPTATION FROM
ANCESTRAL VARIATION IN YOUNG
POPULATIONS OF BEACH MICE
VeraS.Domingues,1,2,3,4 Yu-PingPoh,1,5 BrantK.Peterson,1,2,3 PleuniS.Pennings,1 JeffreyD.Jensen,6
andHopiE.Hoekstra1,2,3
1DepartmentofOrganismic&EvolutionaryBiology,HarvardUniversity,Cambridge02138,Massachusetts
2DepartmentofMolecular&CellularBiology,HarvardUniversity,Cambridge02138,Massachusetts
3MuseumofComparativeZoology,HarvardUniversity,Cambridge02138,Massachusetts
4E-mail:[email protected]
5PrograminBioinformatics&IntegrativeBiology,UniversityofMassachusettsMedicalSchool,Worcester,Massachusetts
6SchoolofLifeSciences,EcolePolytechniqueFe´de´raledeLausanne,Lausanne,Switzerland
ReceivedOctober18,2011
AcceptedMarch26,2012
Tounderstandhoworganismsadapttonovelhabitats,whichinvolvesbothdemographicandselectiveevents,werequireknowl-
edgeoftheevolutionaryhistoryofpopulationsandalsoselectedalleles.Therearestillfewcasesinwhichtheprecisemutations
(andhence,definedalleles)thatcontributetoadaptivechangehavebeenidentifiedinnature;oneexceptionisthegeneticbasis
of camouflaging pigmentation of oldfield mice (Peromyscus polionotus) that have colonized the sandy dunes of Florida’s Gulf
Coast.Toquantifythegenomicimpactofcolonizationaswellasthesignatureofselection,weresequenced50001.5-kbnoncoding
lociaswellasa160-kbgenomicregionsurroundingthemelanocortin-1receptor(Mc1r),agenethatcontributestopigmentation
differences,inbeachandmainlandpopulations.Usingagenome-widephylogeneticapproach,werecoveredasinglemonophyletic
groupcomprisedofbeachmice,consistentwithasinglecolonizationeventoftheGulfCoast.Wealsofoundevidenceofasevere
founderevent,estimatedtohaveoccurredlessthan3000yearsago.Inthisdemographiccontext,weshowthatallbeachsub-
speciesshareasinglederivedlightMc1r allele,whichwaslikelyselectedfromstandinggeneticvariationthatoriginatedinthe
mainland.Surprisingly,wewereunabletoidentifyaclearsignatureofselectionintheMc1rregion,despiteindependentevidence
thatthislocuscontributestoadaptivecoloration.Nonetheless,thesedataallowustoreconstructandcomparetheevolutionary
historyofpopulationsandallelestobetterunderstandhowadaptiveevolution,followingthecolonizationofanovelhabitat,
proceedsinnature.
KEY WORDS: Adaptation,colonization,demography,Mc1r,naturalselection,Peromyscus.
Adaptationfollowingcolonizationofnovelhabitatsisinfluenced ural populations (reviewed in Stinchcombe and Hoekstra 2007;
by the interplay of demographic and selective forces. Thus, to Ellegren and Sheldon 2008; Mackay et al. 2009; Stapley et al.
obtainacompletepictureoftheevolutionarychangeinvolvedin 2010). However, the precise genetic variants that define alleles
the process of adaptation, we must reconstruct the evolutionary responsibleforalteredgenefunctionand/orexpressionhaveyet
historyofallelesunderselectionaswellasthehistoryofpopula- to be identified in most cases (but see Chan et al. 2010) and/or
tions.Overthelastyears,considerableprogresshasbeenmadein theirdirectlinktofitnessisoftenmissing(BarrettandHoekstra
identifying genes that contribute to phenotypic variation in nat- 2011), thus limiting our ability to define and thereby study the
(cid:2)C 2012TheAuthor(s).Evolution(cid:2)C 2012TheSocietyfortheStudyofEvolution.
3209 Evolution66-10:3209–3223
VERA S. DOMINGUES ET AL.
evolution of adaptive alleles. Therefore, fundamental questions et al. 2010) place us in a strong position to accurately infer the
such as the primary source of beneficial mutations (i.e., pre- historyofnaturalpopulationsandrigorouslydocumenttheeffects
existing mutations segregating as standing genetic variation or ofselectionduringlocaladaptation.
denovomutationsthatappearafteranenvironmentalchange[but Here, we capitalize on a system—camouflaging pigmenta-
see Colosimo et al. 2005; Linnen et al. 2009]), the number of tion in the oldfield mouse (Peromyscus polionotus)—for which
independent origins of alleles (see Feldman et al. 2009; Steiner allelesthatcontributetophenotypicvariationhavebeenidentified
etal.2009;Songetal.2011;van’tHofetal.2011),andthetiming (Hoekstra et al. 2006; Steiner et al. 2007) and the link between
andstrengthofselectionactingonthesealleles(seeLinnenetal. phenotypeandfitnessisclear(Vignierietal.2010)toanswerques-
2009)remainlargelyunanswered. tionsregardingtheevolutionaryhistoryoffitness-relatedtraitsin
Answeringthesequestionsalsorequiresplacingthehistory ademographiccontext.Inthedarkloamysoilsofthesoutheastern
offitness-relatedallelesinademographiccontext.Forexample, United States, these mice have dark dorsal coats, whereas con-
only by comparing the evolutionary relationships of beneficial specificsthatsubsequentlycolonizedthepalesanddunesandbar-
alleles with the phylogenetic relationships of populations can riersislandsofFlorida’sGulfCoasthavepalecoatswithreduced
one infer the number of independent origins of beneficial pigmentation. Experimental work showed that natural selection
mutations. Moreover, the ability to distinguish signatures of forcrypsis,viapredationbyvisualhunters,actsondifferencesin
selection from purely demographic effects requires rigorous pigmentationbetweenmainlandandbeachforms(Vignierietal.
estimation of null demographic models (Thornton and Jensen 2010).AlongFlorida’sGulfshoreline,therearefivesubspecies
2007).Specifically,strongorrecentdemographicevents,suchas ofbeachmice;althoughallbeachmicehavelighterpelagethan
populationbottlenecks,canreduceoreveneliminateourability theirmainlandcounterparts,eachsubspeciesshowsauniquecolor
toidentifystatisticallysignificantdeparturesfromneutrality(i.e., pattern(Sumner1926;Mullenetal.2009).Asinglemutationin
selectivesweeps;seereviewofThorntonetal.2007).Inaddition thecodingregionofthemelanocortin-1receptor(Mc1r)hasbeen
to providing a framework for gaining insight into the evolution identifiedandshowntocontributetodifferencesbetweenSanta
of adaptive loci, several aspects of the demographic history RosaIslandbeachmice(SRIBM)(P.p.leucocephalus),thelight-
of colonization have intrinsic value to understand the forces estofthebeachsubspecies,andmainlandmice(P.p.subgriseus)
involved in genetic and phenotypic divergence. Specifically, (Hoekstra et al. 2006). Moreover, variation in the frequency of
when organisms colonize new areas, they typically experience Mc1r alleles among beach mouse subspecies is correlated with
foundereffectsassociatedwithreducedpopulationsizeandlow pigmentdifferences(Mullenetal.2009).
genetic diversity, which may limit opportunities to adapt to a Togainamorecompletepictureoftheevolutionofcryptic
new local environment (Nei et al. 1975). Moreover, ancestral coloration in P. polionotus following the colonization of the
and derived populations may still exchange migrants at least beach habitat, we used a targeted resequencing approach to
in the initial stages of the colonization process—because gene reconstruct the demographic signature of colonization, to infer
flow constrains genetic and phenotypic differentiation between theevolutionaryrelationshipsofadaptivealleles,andtoidentify
populations, this may limit adaptive divergence (Riechert 1993; genomic footprints of selection. Specifically, we evaluated the
Nosil and Crespi 2004; Nosil 2009). Therefore, a complete competing hypotheses of single or multiple colonization events
understandingofadaptiveevolutionnecessitatesinferenceofthe ofthebeachhabitat.Inaddition,weestimatedthetimingofthe
parametersofbothpopulationsizechangeandmigration. colonizationandthenumberofindependentoriginsofaderived
Because patterns of genetic diversity and relationships mutationinMc1rthatcontributestothereducedpigmentationin
among populations vary widely among loci, especially for re- beachmice,andweevaluatedourabilitytodetectsignaturesof
cently diverged species or populations (Tajima 1983), an accu- selectioninarecentlycolonizedpopulationwhilecontrollingfor
rate estimation of the demographic history calls for data from demographiceffects.
multipleunlinkedloci(TakahataandNei1985).Recentdevelop-
mentsinhigh-throughputsequencingtechnologiesandmultiplex
Methods and Materials
strategiesnowmakeitpossibletogenotypehundredsoflocidis-
tributedthroughoutthegenomeinwildpopulations.This,together MOUSE SAMPLES
with novel statistical and computational methods for genome- We collected mice from 13 locations across the range of P.
wide phylogenetic reconstruction (e.g., Guindon and Gascuel polionotus,includingthefiveGulfCoastbeachmousesubspecies
2003;Stamatakis2006),fordemographicmodeling(e.g.,Thorn- andeightmainlandpopulations(Fig.1A).Micewerecollectedus-
tonandAndolfatto2006;Gutenkunstetal.2009;Naduvilezhath ingShermanlivetraps,andatissuesample(tailtiporliver)was
et al. 2011), and for detecting targets of selection (e.g., Sabeti taken and stored in 95% ethanol. We extracted genomic DNA
etal.2002;Nielsenetal.2005;Williamsonetal.2007;Pavlidis fromeachtissuesampleusinganAutogenkitinanAutoGenPrep
3210 EVOLUTION OCTOBER2012
EVOLUTION OF A BENEFICIAL ALLELE IN PEROMYSCUS
Figure 1. (A)MapshowingtheapproximatedistributionofPeromyscuspolionotus(lightgrayshade)inAlabama(AL),Georgia(GA),
South Carolina (SC), and Florida (FL), and sample sites of mainland populations (filled dots) and beach subspecies (open triangles)
surveyedinthisstudy.PiechartsrepresenttheestimatedfrequencyofMc1rancestral“dark”alleles(inblack)andderived“light”alleles
(ingray)inthebeachsubspecies(ABM,Alabamabeachmice;PKBM,PerdidoKeybeachmice;SRIBM,SantaRosaIslandbeachmice;CBM,
Choctawhatcheebeachmice;SABM,SaintAndrewsbeachmice).(B)TopologyofP.polionotuspopulationsgeneratedusinggenomicdata
(211genomicregions,totaling161kb).Bootstrapsupportsbasedon100replicatesareshownonbranches.Taxonlabelscorrespondto
locationsshownin(A).
965 instrument (Autogen, Holliston, MA). Mainland samples genes,weusedsequencesobtainedfromaMc1r-containingBAC
wereaccessionedasmuseumspecimensintheHarvardMuseum clone(aSanger-sequenced160-kbclonefromP.m.rufinusBAC
ofComparativeZoologyMammalDepartment. library at 20× coverage) to design the capture array. We tiled
120merprobes(at5×coverage)acrossintervalsthatencompass
CAPTURE ARRAY DESIGN, TARGET ENRICHMENT, exonsandintronicregionsincludedintheBAC.
AND NEXT-GENERATION SEQUENCING WepreparedbarcodedDNAlibrariesforeachofoursamples.
We designed a custom SureSelect capture array (Agilent We used a multiplex strategy with the above-described array to
Technologies, Santa Clara, CA) to enrich templates for both enrichpoolsof12sampleseach.Captureprobeswerehybridized
random “neutral” regions and for the Mc1r locus (fol- to target DNA in solution, and targeted regions were then se-
lowing Gnirke et al. 2009). To obtain genome-wide data lectedusingmagneticbeadsandamplifiedwithuniversalprimers
for demographic and phylogenetic analysis, we targeted (Gnirkeetal.2009).Wethensequencedourcapturelibrariesand
5114 regions averaging 1.5 kb in length (5.2 Mb of non- generated100-bppaired-endreadsonanIlluminaHiSeq2000(Il-
repetitive sequence) at random locations throughout the Per- luminaInc.,SanDiego,CA).Rawsequencedataareavailableat
omyscus genome at 3× tiling. Probes were designed from an theSequenceReadArchive(accessionnumber:SRA050092.2).
in-house P. maniculatus bairdii ∼1x draft genome assembly
(3.3millionreadsfromBaylorCollegeofMedicine/NCBITrace SHORT-READ SEQUENCE ANALYSES AND
Archive; WGS assembler [Myers et al. 2000]). To maximize VALIDATION
recovery of unique sequence reads, we masked repetitive re- WefirstprocessedrawIlluminadatatogeneratequality-filtered
gionsusingexistingdataformouse/rat(RepeatMaskeropen3.2.7; fileswithsequencesfromeachbarcodedindividual.Second,we
http://www.repeatmasker.org)andcustomrepeatlibrariesforPer- used the Burroughs-Wheeler Alignment (BWA) tool (Li and
omyscusconstructedusingPiler(EdgarandMyers2005)andRe- Durbin 2009) to map reads to a composite reference sequence
peatScout (Price et al. 2005). To resequence an approximately setconsistingofP.m.bairdiigenomescaffoldsandP.m.rufinus
160-kbregioncontainingthesingleMc1rexonandneighboring BAC sequences. We performed initial mapping and alignment
EVOLUTION OCTOBER2012 3211
VERA S. DOMINGUES ET AL.
using default parameters values in BWA (details of final condi- optimumof0.54inpasstwo).Wethenremovedregionsinwhich
tions are explained below). Variant discovery and diploid geno- the average fraction of heterozygote individuals across sites in
typeinferencewereperformedusingtheGATKsoftwarepackage, thatregionexceeded0.54.
whichsimultaneouslyestimatesindividualgenotypesandpopu-
GENERATING FINAL DATASETS
lationallelefrequencies(DePristoetal.2011).
To generate final datasets, we performed mapping using the se-
To assess accuracy and completeness of genotype determi-
lectedBWAparameters(l=55,k=4,R=10)andremovedpar-
nation, we performed a receiver-operating characteristic (ROC)
alogousregions.BasedonROCcurveinflections,weselectedtwo
analysistocomparegenotypesofeightmicegeneratedfromour
thresholdsetconditionsforvariantcallinginGATKthatproduced
short-readsequencedatawiththoseobtainedfromatotalof690
reliable genotypes. We generated a “dataset A” using QD = 6,
kb of Sanger sequence for the same individuals. Regions sam-
QUAL=200,andGQ=12(99.6%concordancewithSangerdata
pled in Sanger sequencing were representative with regards to
and false-positive rate of 0.0019), and “dataset B” using QD =
GC content and short read coverage of the Mc1r locus overall
20 and GQ = 20 (concordance 99.8% and false-positive rate
(targetedintergenicregionshadequivalentGCcontentbutlower
0.0004). Although parameters in dataset A are more permissive
meancoverageasexpectedfromdifferentialprobetillingdepth).
thanindatasetB,theystillgeneratereliablegenotypesandhave
Wecalculatedareaunderthecurve(AUC)forthresholdmodels
the advantage of providing more power for analysis based on
offourgenotypequalitymetrics(describedinSupportinginfor-
patternsofgeneticvariability.Finally,weimputedmissinggeno-
mation): two individual-level metrics (mapped read depth, DP
types and estimated haplotype phase using a statistical model
andgenotypequality,GQ)andtwopopulation-levelmetrics(site
based on flexible clustering of patterns of genetic variation and
quality,QUALandvariantquality-by-depth,QD).Comparisons
linkagedisequilibriuminnaturalpopulations,implementedinthe
of AUC showed that population-level metrics outperform indi-
software package fastPHASE (Scheet and Stephens 2006), and
vidual metrics, confirming the utility of simultaneous genotype
obtained SNP-specific linkage disequilibrium-based genotyping
inferenceacrossmultipleindividuals(Fig.S1).
errorrates(followingScheetandStephens2008).
OPTIMIZING MAPPING AND FILTERING OUT
ANALYSIS OF POPULATION GENETIC DIVERSITY
PARALOGOUS REGIONS
AND DIFFERENTIATION
Havingfoundthebestmetricstoassessaccuracyofgenotyping,
Using genome-wide SNP data, we first estimated average het-
we returned to the question of optimizing mapping parameters.
erozygosityforeachpopulationasameasureofgeneticdiversity.
Because our reference assemblage consists of sequences from
To include only unlinked markers in our analysis, we used one
a different (although closely related) species, we tested perfor-
sitepergeneticregionrecoveredindatasetA,totaling2236SNPs.
manceofmappingwithdifferentparametersvalues(seedlength,
Next, we estimated Fst in the program SPAGeDi (Hardy and
l;seededitdistance,k;andnumberofsuboptimalalignmentsto
Vekemans 2002) to assess the extent of differentiation among
sample,R)inBWAusingROCcurvesofpopulation-levelmetrics
beachmiceandmainlandpopulations.Inaddition,weestimated
asdescribedabove.WetestedarangeofparametersaroundBWA
population structure among beach mouse subspecies using a
defaultssinglyandincombination(l=35,55,75;k =2,4,6;
Bayesian model based cluster method implemented in STRUC-
R=10,100)andconcludedthattheparametersetthatmaximized
TUREversion2.3(Pritchardetal.2000)thatassignsindividuals
numberofreadsmappedandconcordanceoffinalgenotypeswas
topopulationsbasedontheirgenotypefrequencies.Wealsoesti-
l=55,k=4,R=10.
matedthemaximum-likelihood(ML)numberofsubpopulations
Asthereferencegenomeisincomplete,weexpectsomeinci- (k=1–6)usinganadmixturemodelandassuminguncorrelated
denceofparalogousmappinginouralignedreads.Onesignature
allelefrequenciesbetweenclusters.First,werantheanalysis10
of paralogous mapping is an excess of heterozygous genotypes
timesforeachkusingaburn-inof100,000MarkovChainMonte
(which are actually divergent sites between paralogs) compared
Carlo (MCMC) steps followed by 100,000 iterations. Then, we
tofrequenciesofthetwohomozygousgenotypeclasses.Tofilter
used the Greedy algorithm in CLUMPP (Jakobsson and Rosen-
out paralogous regions, we computed the fraction of heterozy-
berg2007)tohandletheresultsfromreplicateanalyses.Finally,
gousindividualsateachvariablesiteforthetwovariantcalling
we compared posterior probabilities of each k to determine the
models used to generate the final datasets (see below). In brief,
numberofsubpopulationsthatbestfitourdata.
we tested independently the two final models with the addition
of a fixed threshold for mean fraction heterozygous genotypes ESTIMATION OF THE P.POLIONOTUS POPULATION
rangingfrom0to1instepsof0.1(passone),andsubsequently TREE
inintervalsof0.01(passtwo),betweenthebestpass-onemodels To reconstruct the colonization history of mice along the
(mean heterozygous fractions of 0.5 and 0.6 in pass one, final Floridacoast,weestimatedapopulationtopologyusinggenomic
3212 EVOLUTION OCTOBER2012
EVOLUTION OF A BENEFICIAL ALLELE IN PEROMYSCUS
sequencing data generated with the most restrictive of our pa-
rameter sets (dataset B) to ensure that only extremely reliable mainland
genotypes were used for phylogenetic construction. We used n2
the Maximum Pseudo-Likelihood for Estimating Species Tree
(MP-EST) method (Liu et al. 2010) to estimate a “species”
tree of P. polionotus under the coalescence model by maxi- s m12
mizing a pseudolikelihood function using a set of gene trees,
T
which accounts for gene tree heterogeneity. First, we se- m21
beach
lected sequences from regions across the genome that con-
tained at least five polymorphic sites (211 regions, total of
161 kb) and produced 100 bootstrap replicates for each region
n1
usinganonparametrictechnique(Efron1981).Second,webuilt
genetreesforeachbootstrapsampleusingaML-basedmethod
time
with a GTRGAMMA model of nucleotide substitution imple-
mented in the program RAxML (Stamatakis 2006). Gene trees
wererootedusingP.maniculatus,thesisterspeciesofP.poliono- Figure 2. Illustrationoftheisolation-with-migrationmodelused
tus,asanoutgroup.Next,weusedtherootedgenetreestocon- toestimatesixdemographicparameters.Inthismodel,afraction
of the ancestral mainland population(s) gives rise to a derived
struct100MP-ESTtrees.Finally,aconsensustreewasbuiltfrom
beachpopulationattimeT inthepast.Mainlandandbeachpop-
the 100 MP-EST trees using Majority-Rule-extension (MRe) in
ulations differentiate with bidirectional migration occurring be-
CONSENSEfromthePHYLIPpackage(Felsenstein1989).The
tweenthem(m21andm12).Finalbeachandmainlandpopulation
tipsoftheresultingtreerepresenttheeightmainlandpopulations
sizesaren1andn2,respectively.
andfivebeachsubspeciesusedinthisstudy.
ulation, while the two populations exchange migrants. The fol-
DEMOGRAPHY OF BEACH COLONIZATION lowingparameterswereestimated:sizeofthebeachpopulations
Weinferredpopulationhistoryofbeachmiceusinggenome-wide immediatelyfollowingthesplit(s,givenasafractionoftheorig-
SNPdata(datasetA)tomaximizethenumberofvariablesitesand inal population), final size of beach (nu1) and mainland (nu2)
thuspower.WeusedthemainlandpopulationfromAlabama(AL- populations, time of the split (T), migration from mainland to
2),whichwasmostcloselyrelatedtothebeachsubspeciesfrom beach (m21), and migration from beach to mainland (m12). To
the phylogenetic analysis described above, to represent the an- estimateparameteruncertainty,weusedconventionalbootstrap-
cestralpopulation.Forcomparison,wealsorepeatedtheanalysis ping(fitting100datasetsresampledoverloci)andtoensurethat
usingasecondcloselyrelatedpopulationfrommainlandFlorida, wereachedconvergence,weraneachdatasetatleastfivetimes.
FL-2.Demographicparameterswereestimatedfromthejointal- Wethendeterminedtheconfidenceintervalofeachparameteras
lelefrequencyspectrum(AFS;thedistributionofpolymorphism the 2.5 and 97.5 percentiles of the parameter distribution in the
frequenciesinthepopulations)byapplyingadiffusion-basedap- bootstrappeddata.
proachimplementedintheprogram∂a∂i(Gutenkunstetal.2009). Wenexttestedaless-pronouncedfoundereventbycompar-
WepolarizedSNPsusingP.maniculatusasanoutgroup.Because ingthelikelihoodsoftheoptimizedmodelwithamodelofequal
the sample size needs to be the same for each SNP when con- populationsizeofbeachandmainlandpopulationsfollowingthe
structingtheAFS,weusedtheprojectionoptionin∂a∂itoadjust populationsplit.Toexplicitlytestdifferentmigrationscenarios,
the number of individuals included to maximize the number of weperformedlikelihoodratiotestsbetweentheoptimizedmodel
segregatingsitesused.Thefinaldatasetconsistedof:onemain- andthemodelwithalternativemigrationrates(nomigrationand
landpopulation(AL-2;N=4),SRIBM(N=15),Alabamabeach symmetrical migration with low and high migration rates) and
mice(ABM;N =5),PerdidoKeyBeachmice(PKBM;N =4), determinedP-valuesforeachmodelcomparisonassumingaχ2
Choctawhatchee beach mice (CBM; N = 4), and St. Andrews distributionoftheteststatistics.
beach mice (SABM; N = 5). The number of segregating sites Finally,weconvertedthetimingofthepopulationsplit(given
ranged from 163 to 495 depending on the subspecies analyzed. inunitsof2Ngenerations)toyears.WefirstobtainedNfromθ=
Weestimateddemographicparametersinapairwisefashionus- 4NμL , where μ is the mutation rate per site per generation
eff
ing the mainland population with each of the beach subspecies. estimated for Mus domesticus (μ = 3.7 × 10−8, Lynch 2010),
TheAFSofthemainlandpopulationandeachofthebeachsub- and L is the effective length of the genomic region used (L
eff eff
specieswerefittedtoanisolation-with-migrationmodel(Fig.2). = 349,390 bp). Peromyscus polionotus achieves sexual matu-
Inthismodel,themainlandpopulationgivesrisetoabeachpop- rityatabout30days,andgestationlastsapproximately30days
EVOLUTION OCTOBER2012 3213
VERA S. DOMINGUES ET AL.
(Clark 1938). In the wild, these mice reproduce year-round, al- positiveselectioninMc1r,weusedseveralapproaches.First,we
thoughthereisadeclineinbreedingactivityduringthesummer. performedasliding-windowanalysisofdivergenceandpolymor-
WeestimatedtimeinyearsusingourestimateofNandaconser- phism using a modification of the Hudson–Kreitman–Aguade
vativetwotothreegeneration-per-yearestimate. (HKA) test (Hudson et al. 1987). This analysis was conducted
in the SRIBM population (which is fixed for the derived Mc1r
PHYLOGENETIC ANALYSIS OF PIGMENTATION 65Cys mutation) and the ABM population (fixed for the ances-
ALLELES tral65Argstate).Bycomparingthesetwobeachpopulations,we
Followingtheapproachforestimatingpopulationtrees,weused hopedtominimizedifferencesingeneticdiversityduetoeffec-
themostrestrictiveofourparametersets(datasetB)togenerate tive population size (reduced following colonization) alone. To
genealogiesofMc1r alleles.WefirstconfirmedMc1r genotype identify regions that deviate from neutralexpectations, we used
at the Arg65Cys site for all samples using a TaqMan assay (de- χ2 teststatisticsinwhichtheobservednumbersofpolymorphic
scribed in Steiner et al. 2007). To determine the number of in- andfixeddifferencesarecontrastedwiththeirexpectednumbers
dependentoriginsofthederivedMc1rallelesthatcontainedthe in 10-kb sliding windows along the 160-kb region surrounding
65Cys mutation (“light allele”), we estimated a ML genealogy Mc1r.Second,wefocusedontheSRIBMpopulationandapplied
usinga4-kbfragmentincludingtheMc1rexonanditsneighbor- a composite likelihood method, SweepFinder, to locate regions
ing sequence (comprising 35 parsimony informative sites). We withskewedsitefrequencyspectrum(i.e.,deviantpatternsinal-
usedtheGTRGAMMAmodelofnucleotidesubstitutionandap- lele frequency) in Mc1r and surrounding regions (Nielsen et al.
plied the rapid bootstrapping algorithm in RAxML (Stamatakis 2005).Third,wecharacterizedpatternsoflinkagedisequilibrium
2006).Inaddition,wegeneratedMPandBayesianphylogenies using the ω statistics (Kim and Nielsen 2004), via a recently
using the software Paup∗ version 4.0b10 (Swofford 2002) and proposed modification of SweepFinder (Pavlidis et al. 2010)—
MrBayesversion3.2(RonquistandHuelsenbeck2003),respec- an addition that was shown to greatly increase discriminatory
tively.Nodesupportwasevaluatedwith1000bootstrapreplicates power—to localize selective sweeps. Finally, we used XP-CLR
(Felsenstein1985)inthecaseofMLandMPmethodsandpos- (Chen et al. 2010) to compare patterns of allele frequency dif-
terior probabilities in the Bayesian analysis.Totest if theMc1r ferentiationbetweenSRIBMandABM,alongthe160-kbregion
treetopologywasstatisticallydifferentfromthepopulationtree surroundingMc1r,thatcouldbeindicativeofselectivesweeps.
topology, we conducted an SH test (Shimodaira and Hasegawa
1999)tocompareMLscoresforatopologyinwhichbeachmice
Results
wereforcedtobemonophyleticandthebestMLtopologyforthe
Mc1r fragment. The SH test was performed in Paup∗ using the TARGETED ENRICHMENT (SURESELECT)
RELLapproximationwith10,000bootstrapreplicates. PERFORMANCE
Because the Mc1r gene topology is expected to reflect se- A total of 601 million paired-end reads were generated in
lection,butthesignatureofselectionshouldbelostaswemove three HiSeq lanes, of which 99.6% could be confidently as-
away from the target of selection due to recombination, we es- signed to individual barcodes. Overall, 25–45% (mean 33%)
timated genealogies using sequences at varying distances from of reads mapped to target regions. For the regions targeted
the Mc1r gene. To first annotate the Peromyscus Mc1r BAC, at random locations throughout the Peromyscus genome, af-
we aligned Peromyscus sequences to the Rattus genome using ter removing putative multicopy regions, a total of 3754 and
GenomeVISTA (Bray et al. 2003; Couronne et al. 2003). Next, 3884 regions were retained applying the thresholds used in
we chose fragments at different distances from Mc1r that con- datasets A (QD = 6, QUAL = 200, and GQ = 12) and B
tained a similar number of informative sites. All phylogenetic (QD = 20 and GQ = 20), respectively. In dataset A, a total of
reconstructionswereperformedasdescribedabove. 4billionbasesmappedon-targetandaverageindividualcoverage
of target regions ranged from 3× to 38× (average 11×), corre-
DETECTING SIGNATURES OF POSITIVE SELECTION IN spondingtoamean46-foldenrichment.IndatasetB,atotalof8
MC1R billion bases mapped on-target and average individual coverage
Usingthedatasetthatprovidesincreasedpowerforgeneticdiver- of target regions ranged from 5× to 70× (average 20×), corre-
sitybasedanalysis(datasetA),wefirstmeasuredgeneticdiversity spondingtoamean84-foldenrichment.Thedifferenceinnumber
(e.g.,numberofsegregatingsites,averagenumberofpairwisedif- oftotalbasesmappedbetweendatasetAandBreflectsthedif-
ferences,andWaterson’sθ)andcharacterizedtheskewinthefre- ferenceinthestringencythresholdsusedtofilteroutparalogous
quencyspectraofthe160-kbresequencedregionincludingMc1r, regions.However,becausetheadditional4billionbasesindataset
bycalculatingTajima’sD(Tajima1989)andFuandLi’sD∗ (Fu Bcorrespondtodenselycoveredregions,theycorrespondtoonly
andLi1993)foreachbeachpopulation.Totestforsignaturesof 4%ofthetotalvariantsitesinthedataset.
3214 EVOLUTION OCTOBER2012
EVOLUTION OF A BENEFICIAL ALLELE IN PEROMYSCUS
Table 1. Demographic parameters inferred from each of the five beach mouse subspecies. Values are ML parameter estimates for
anisolation-with-migrationmodelfittedtothejointallelefrequencyspectrumofeachofthefivebeachmousesubspecies,usingthe
mainlandpopulationmostcloselyrelatedtothebeach(AL-2)astheancestralpopulation.Bootstrap95%confidenceintervalsareshown
inparentheses.
ABM PKBM SRIBM CBM SABM
s 0.020 0.005 0.001 0.016 0.090
(0.000–0.213) (0.000–0.238) (0.001–0.225) (0.000–0.423) (0.000–0.249)
n1 1.877 1.455 0.413 41.487 0.185
(0.400–6.848) (0.307–198.878) (0.130–0.509) (1.526–275.570) (0.185–2.912)
n2 117.170 61.320 222.080 64.274 698.040
(46.787–996.389) (23.636–997.544) (57.580–960.281) (21.946–987.185) (36.078–996.447)
T 1.399 1.428 1.225 1.733 2.537
(1.121–4.292) (1.085–3.080) (0.883–1.769) (1.034–2.832) (1.392–3.617)
T (years) 2787 3287 3040 2771 3,177
(2233–8548) (2499–7090) (2191–4391) (1654–4530) (1,744–4,528)
m21 0.072 0.124 0.203 0.148 0.236
(0.021–0.122) (0.075–0.289) (0.121–0.495) (0.001–0.356) (0.015–0.236)
m12 0.017 0.027 0.018 0.000 0.016
(0.000–0.034) (0.000–0.078) (0.000–0.029) (0.000–0.025) (0.000–0.023)
θ 102.983 119.043 128.322 82.708 64.762
lnL –86.591 –85.107 –128.170 –77.808 −92.860
s,fractionoftheancestralpopulationthatgaverisetothebeachpopulation;n1,finalsizeofthebeachpopulationrelativetotheancestralpopulation
(reference);n2,finalsizeofthemainlandpopulationrelativetotheancestralpopulation;T,timeinthepastofsplit,inunitsof2Nrefgenerationsandin
years;m12,effectivemigrationratefrombeachtomainland,inunitsof2Nrefpergeneration;m21,effectivemigrationratefrommainlandtobeach,inunits
of2Nrefpergeneration;θ,4Neμ(Ne,effectivepopulationsize;μ,mutationrate);lnL,logarithmofthemodelLikelihood.
Atotalof189millionbaseswererecoveredinthetargeted To reconstruct the evolutionary history of beach mice, we
160-kb region including Mc1r and surrounding sequences. Av- estimatedapopulationtopologyofP.polionotususingagenome-
erage individual coverage of target regions ranged from 5× to wide phylogenetic approach, which accounts for gene tree het-
49× (average 17×), corresponding to a mean 72-fold enrich- erogeneity (Fig. 1B). Beach mice form a well-supported mono-
ment.GenotypecallaccuracyisshowninTableS1.SNP-specific phyleticgroup,withtheclosestrelatedpopulationfromsouthern
linkage disequilibrium-based genotyping error rates associated Alabama. These data suggest a single colonization event of the
with phasing genotypes ranged from 0 to 0.1823 (median = coastfromnorthernmainlandpopulations,followedbydifferen-
0.0001). tiationamongthefivebeachmousesubspecies.
POPULATION GENETIC DIVERSITY AND DEMOGRAPHY OF BEACH COLONIZATION
DIFFERENTIATION Weestimateddemographicparametersforthefivebeachmouse
Wefirstcharacterizedpatternsofgeneticdiversityanddifferen- subspecies. Parameter estimates are similar across each of the
tiation among populations. Not unexpectedly, beach mice show beach subspecies (Table 1). In all cases, we obtained a recent
lower levels of genetic diversity (average heterozygosity, H = splitbetweenbeachandmainlandpopulations(1800–3300years;
0.14–0.19)thanmainlandpopulations(H=0.22–0.32).Wecon- 95%confidenceinterval:1600–8500years).Moreover,estimates
firmedwithgenomicdatathatbeachmiceontheGulfCoastare oftheeffectivepopulationsizeofbeachpopulationsimmediately
welldifferentiatedintofivesubspecies.Fstvaluesrangedbetween afterthesplitareextremelysmallinallcases,rangingfrom0.01%
0.38and0.58,andtheSTRUCTUREanalysissupportedfivedis- to9%oftheancestralpopulation.Fouroffivebeachsubspecies
tinctgeneticclusters,correspondingtotherecognizedsubspecies. rejected the null hypothesis of equal size beach and mainland
All individuals were correctly assigned to the cluster with their populations immediately following the split (Table 2). The one
subspeciesdesignation(Fig.S2).Fstvalueswerealsohighamong exceptionisSABM,thesubspecieswiththelargestestimateofs
mainlandpopulations(0.09–0.31),suggestingthatmainlandpop- (i.e.,estimatedfractionoftheancestralpopulationthatgaverise
ulationsarealsogeneticallydifferentiated. tothebeachpopulation),whichshowedaverysmalldifference
EVOLUTION OCTOBER2012 3215
VERA S. DOMINGUES ET AL.
Table 2. ML ratio tests of four demographic models. For each 23 Mc1r alleles that contain the 65Cys mutation form a derived
beachmousesubspecies,thedifferencein–lnLbetweentheML monophyleticclade,indicatingacommonancestryofMc1rlight
model and the four models tested (listed in the first column) is alleles.ThisMc1r“selected”genealogycontrastswiththepopula-
given.
tion“neutral”topology(Fig.1B),asnotallallelesfoundinbeach
miceclustertogether.Indeed,atopologyforcingallbeachmice
ABM PKBM SRIBM CBM SABM
alleles to be monophyletic is significantly worse than the Mc1r
Largefounder 83.5∗ 72.5∗ 92.1∗ 60.9∗ −1.78 MLtopology(SHtest,P<0.05).TheMc1rtopologyisrecovered
populationsize
in genealogies that are 2-kb, 10-kb, and 20-kb downstream and
(s=0.5)
upstream of Mc1r (data not shown). However, trees 35-kb and
Nomigration 163.7∗ 116.1∗ 204.8∗ 215.8∗ 47.3∗
40-kb downstream of Mc1r (Fig. 3, blue and green panels) are
(m=0)
Symmetrical 8.8∗ 12.1∗ 27.9∗ 24.7∗ 65.3∗ similartothepopulationtreetopologyasthevastmajority(more
migration than90%)oftheallelesfoundinthebeachpopulationsfallina
(m21=m12; singleclade.
higherestimate)
Symmetrical 7.2∗ 16.0∗ 60.6∗ 196.5∗ 19.3∗ DETECTING SIGNATURES OF POSITIVE SELECTION IN
migration
MC1R
(m21=m12;
By comparing patterns of nucleotide variation between the
lowerestimate)
SRIBM, in which Mc1r light alleles are fixed, and ABM, in
whichancestraldarkallelesarefixed,wecouldtestforpatterns
Allbutonecomparisonissignificant(∗P-values<0.001).
consistent with selection acting on Mc1r. Overall, the complete
160-kbMc1rsequencesinbeachmiceshowedlowlevelsofpoly-
in –ln L between the ML model and the model with larger s, morphismandnegativevaluesforTajima’sDandFu’sandLi’s
limitingourpowertorejectthe“equalsize”model.Overall,these D∗duetoanexcessofrarepolymorphisms(Table3).Comparison
results support a strong founder event associated with the colo- ofpatternsofgeneticvariationamongthefivebeachsubspecies
nizationofthebeachhabitat. reveals that SRIBM shows the most skewed frequency spectra,
We detected evidence of asymmetrical migration between whereasABMbestfitstheneutralequilibriumexpectation.This
mainlandandbeachpopulationsfollowingcolonization.Migra- pattern is also evident when comparing haplotypes (Fig. 4), as
tionfromthemainlandintothebeachpopulationsishigherthan most polymorphisms among light Mc1r alleles are singletons,
in the opposite direction for all five cases. The models with no whereasvariantsamongancestralallelesaresegregatingathigher
migration or with symmetrical migration rates were rejected by frequency.Theplotofnucleotidepolymorphismanddivergence
likelihoodratiotests(Table2). alongthe160-kbregionshowshigherpolymorphismthandiver-
Finally,demographicresultsobtainedusingasecondmain- genceintheMc1rregioninSRIBM,whereastheoppositeistrue
landpopulation,FL-2,werelargelyconsistentwithestimatesob- forABM(Fig.4).Althoughthispatternisconsistentwithpositive
tainedwiththefirstmainlandpopulation,AL-2(TableS2).Any selection,theincreaseddivergenceattheMc1rgeneisnotpartic-
differencesinmigrationratecanbeattributedtodifferentiationof ularlystrikingwhencomparedwiththeentireregion,asseveral
themainlandpopulationsfollowingcolonization. regions show biased polymorphism-to-fixed differences ratios,
likelyduetodemographiceffects.Inaddition,methodsbasedon
EVOLUTIONARY HISTORY OF PIGMENTATION skewsinAFS(SweepFinder),patternsoflinkagedisequilibrium
ALLELES (ω statistics) or allele frequency differentiation did not show a
ToreconstructtheevolutionaryhistoryofthelightMc1rallele(de- clearsignatureofselectioninMc1r.Thus,weidentifiedpatterns
finedbythederived65Cysmutation),wefirstdetermineditsfre- ofvariationconsistentwithpositiveselectionactingonthelight
quencyamongsubspeciesandthengeneratedgenealogiesbased Mc1ralleles,althoughnonewerestatisticallysignificant.
ongenomicfragmentsatvaryingdistancefromMc1r.Frequen-
ciesofthe65CysMc1rmutationinthefivebeachmicesubspecies
Discussion
wereconsistentwithpreviousindependentresults(Mullenetal.
2009). The derived mutation was fixed in SRIBM, in high fre- Usinghigh-throughputsequencingoftargetedregions,wetakea
quency in PKBM (60%) and SABM (75%), low frequency in comprehensiveapproachtogainnovelinsightintotheevolution
CBM(10%),andabsentinABM(Fig.1A). ofcamouflagingcolorationofbeachmice,awell-knowncaseof
ThegenealogybuiltusingthefragmentthatincludesMc1r adaptationinthewild.Werecreatedthehistoryofbeachmouse
and neighboring sequences (Fig. 3, red panel) shows that the populations in parallel with the evolutionary history of Mc1r
3216 EVOLUTION OCTOBER2012
EVOLUTION OF A BENEFICIAL ALLELE IN PEROMYSCUS
Figure 3. Schematic representation of a genomic region containing Mc1r and neighboring sequences and ML trees estimated from
sequencesalongthatregion.Allelesfoundinmainlandaredepictedascircles;allelesfoundinthebeachastriangles.Redpanel:Mc1r
genealogyof71alleles,estimatedusinga4-kbsequenceincludingthe954-bpMc1rexon(notehowthistopologydiffersfrompopulation
treetopologyshowninFig.1B).Lightalleles,definedbythederived65Cysmutation(redtriangles),formamonophyleticclade,indicating
theircommonorigin.Blueandgreenpanels:Genealogiesestimatedfromregionsapproximately35-and40-kbdownstreamof Mc1r,
showingtopologiessimilartothepopulationtree,inwhichallelesfoundinbeachpopulations(regardlessoftheirgenotypeatMc1r
position65)fallinthesameclade.
alleles, a pigment gene that contributes to cryptic coloration in Table 3. Genetic diversity indices of Mc1r and its neighboring
these mice. We first show that mice likely colonized the novel regions(∼160kb)inthefivebeachmousesubspecies.
beachhabitatlessthan3000yearsagoinasinglefounderevent.
ABM PKBM SRIBM CBM SABM
Inthisdemographiccontext,wedemonstrateasingleoriginofa
No.ofsegregating 806 940 876 981 481
derived Mc1r mutation that contributes to light pigmentation in
sites
severalsubspeciesofbeachmice,andthatthederivedalleleisold,
Averagepairwise 279.293 318.039 79.901 338.168 113.349
likelyfirstarisinginthemainlandpopulations.Further,although differences
we showpatterns ofgeneticvariation inthe regionsurrounding Watterson’sθ 284.909 332.277 208.493 346.770 170.027
Mc1rconsistentwithrecentpositiveselection,wedidnotdetect Tajima’sD −0.099 −0.215 −2.332 −0.124 −1.670
FuandLi’sD∗ −0.331 −0.682 −4.082 −0.723 −2.087
astatisticallysignificantsignatureofselection,likelyduetothe
demographichistoryofthesepopulations.Togethertheseresults
retracetheevolutionaryhistoryofabeneficialalleleinademo- flow.Giventhisresult,thecolonizationoftheFloridacoastiscon-
graphiccontext,contributingtoourunderstandingoftheadaptive sistentwithtwohypotheses:(1)mainlandmicecolonizedthecoast
processinthewild. in a single event, followed by subsequent population differenti-
ation;or(2)beachmousesubspeciesresultedfromindependent
invasions, possibly from multiple mainland source populations.
POPULATION HISTORY Totestthesealternativehypotheses,weinferredthephylogenetic
With genome-wide polymorphism data, we confirmed previous relationshipsof13P.polionotuspopulationsusing211genomic
findings(Mullenetal.2009)thatthefivebeachmousesubspecies regions (161 kb) resequenced in a total of 71 individuals. Be-
are highly differentiated and show no evidence of current gene causegenetreescanshowdiscordanttopologiesduetovariation
EVOLUTION OCTOBER2012 3217
VERA S. DOMINGUES ET AL.
Figure 4. GeneticvariationinMc1randflankingregionsforSRIBM(derivedmutationinMc1risfixed,n=12alleles)andABM(derived
mutationinMc1risabsent,n=12alleles).Ontheleft,haplotypesfora10-kbregionsurroundingMc1r(gray-shadedregioninthegraphs
totheright).Rowsareobservedhaplotypes;columnsarevariablenucleotidesites(blackindicatesancestral,redindicatesderivedstate,
grayindicatesunknowngenotype).NucleotidestateswerepolarizedusingthemainlandpopulationAL-2.Graphsontherightareplots
ofpolymorphismanddivergenceacrossa50-kbregionincludingMc1r.
incoalescence,especiallyindivergingpopulationsorincasesof previous age estimate based on mtDNA, suggesting that beach
recent speciation (Tajima 1983), we used a coalescence model micedivergedmorethan200,000yearsagofrommainlandmice
to account for gene tree heterogeneity (Liu et al. 2010). This (VanZantandWooten2007).Thispreviousanomalousresultis
methodallowedustoconfidentlydetermineifbeachmiceform likelyduetothelow-resolutionphylogeneticreconstructionbased
a monophyletic group, or instead, fall into distinct clades. Our onasinglemtDNAmarker,andspecificallyfromestimatingadi-
genome-widepopulationtopologyclearlyshowsastronglysup- vergencetimefora“beachclade”thatincludedmainlandmice,
ported monophyletic clade of beach mice (Fig. 1B). Based on leading to a large overestimate of divergence time. Our recent
theseresults,themostlikelyscenarioisthatmainlandmicecol- splittingtimeestimatesuggeststhatdifferencesincolorationbe-
onizedthe sandy beachesof Florida’s Gulf Coast once andthat tweenmainlandandbeachmicemighthaveevolvedquiterapidly.
thisfoundingpopulationsubsequentlydifferentiatedtoyieldthe Thus,theevolutionofcamouflagingpigmentationinbeachmice
fivecontemporarysubspecies. adds to the growing evidence of natural selection and adaptive
divergence occurring on ecological time scales in different taxa
DEMOGRAPHY OF BEACH COLONIZATION suchasinfish(Belletal.2004;Elmeretal.2010),lizards(Losos
Usingpolymorphismsobtainedfromtargetedresequencingofap- et al. 2004; Herrel et al. 2008), and birds (Hendry et al. 2006).
proximately5000genomicregionsandapowerfulmethodbased Second,weinferredastrongfoundereventassociatedwithbeach
onthemultipopulationallelefrequencyspectrum(AFS),wewere colonization.Accordingtoourmodel,only0.1–9.0%ofthege-
abletoestimateseveraldemographicparametersassociatedwith neticvariationpresentinthemainlandpopulationwascapturedin
thecolonizationofthebeachhabitat.First,weestimatedthetime thebeachpopulationatthetimeofdivergence.Suchanextreme
ofthecolonizationevent(i.e.,timeofthepopulationsplit)tobe reduction in diversity is consistent with multilocus estimates of
between1800and3300yearsago,suggestingtherecentcoloniza- bottlenecks associated with colonization of novel areas in a va-
tionofthecoast.Thisisaftertheformationofthebarrierislands, riety of other species (e.g., Rosenblum et al. 2007; Peters et al.
estimatedtobeapproximately6000yearsoldbasedongeological 2008; Elmer et al. 2010), suggesting that in nature new pop-
data(McNeil1950).However,thisfindingstronglycontradictsa ulations often are established and adapt to novel environments
3218 EVOLUTION OCTOBER2012