This is a sel -a chi ed e sion o an o iginal a icle. This e sion
may di e om he o iginal in pagina ion and ypog aphic de ails.
Au ho (s):
Ti le:
Yea :
Ve sion:
Copy igh :
Righ s:
Righ s u l:
Please ci e he o iginal e sion:
CC BY 4.0
h ps://c ea i ecommons.o g/licenses/by/4.0/
Cold adap a ion d i es popula ion genomic di e gence in he ecological specialis ,
D osophila mon ana
© 2021 he Au ho s
Published e sion
Wibe g, R. A. W.; Tyukmae a, V.; Hoikkala, A.; Ri chie, M. G.; Kanka e, M.
Wibe g, R. A. W., Tyukmae a, V., Hoikkala, A., Ri chie, M. G., & Kanka e, M. (2021). Cold
adap a ion d i es popula ion genomic di e gence in he ecological specialis , D osophila
mon ana. Molecula Ecology, 30(15), 3783-3796. h ps://doi.o g/10.1111/mec.16003
2021
Molecula Ecology. 2021;30:3783–3796.
|
3783wileyonlinelib a y.com/jou nal/mec
Recei ed: 22 Ap il 2020
|
Re ised: 10 May 2021
|
Accep ed: 20 May 2021
DOI: 10.1111/mec.16003
ORIGINAL ARTICLE
Cold adap a ion d i es popula ion genomic di e gence in he
ecological specialis , D osophila mon ana
R. A. W. Wibe g1 | V. Tyukmae a2 | A. Hoikkala2 | M. G. Ri chie1 |
M. Kanka e2
This is an open access a icle unde he e ms o he C ea i e Commo ns A i bu ion License, which pe mi s use, dis ibu ion and ep oduc ion in any medium,
p o ided he o iginal wo k is p ope ly ci ed.
© 2021 The Au ho s. Molecula Ecology published by John Wiley & Sons L d.
1Cen e o Biological Di e si y, School
o Biology, Uni e si y o S And ews, S
And ews, UK
2Depa men o Biological and
En i onmen al Science, Uni e si y o
Jy äskylä, Jy äskylä, Finland
Co espondence
R. A. W. Wibe g, Cen e o Biological
Di e si y, School o Biology, Uni e si y
o S And ews, S And ews, KY16 9TH,
Sco land, UK.
Email: aw.wibe[email p o ec ed]
P esen add ess
R. A. W. Wibe g, Depa men o
En i onmen al Sciences, Zoological
Ins i u e, Uni e si y o Basel, Basel,
Swi ze land
V. Tyukmae a, Cen e d'Ecologie,
Fonc ionelle e E olu i e, CNRS,
Mon pellie , F ance
Funding in o ma ion
Ella ja Geo g Eh n oo hin Sää iö; Academy
o Finland, G an /Awa d Numbe : 267244,
268214 and 322980; Na u al En i onmen
Resea ch Council, G an /Awa d Numbe :
NE/L501852/1 and NE/P000592/1
Abs ac
De ec ing signa u es o ecological adap a ion in compa a i e genomics is challenging,
bu analysing popula ion samples wi h cha ac e ised geog aphic dis ibu ions, such as
clinal a ia ion, can help iden i y genes showing co a ia ion wi h impo an ecologi-
cal a ia ion. He e, we analysed pa e ns o geog aphic a ia ion in he cold- adap ed
species D osophila mon ana ac oss pheno ypes, geno ypes and en i onmen al condi-
ions and es ed o signa u es o cold adap a ion in popula ion genomic di e gence.
We i s de i ed he clima ic a iables associa ed wi h he geog aphic dis ibu ion o
24 popula ions ac oss wo con inen s o ace he scale o en i onmen al a ia ion
expe ienced by he species, and measu ed a ia ion in he cold ole ance o he lies
o six popula ions om di e en geog aphic con ex s. We hen pe o med pooled
whole genome sequencing o hese six popula ions, and used Bayesian me hods o
iden i y SNPs whe e gene ic di e en ia ion is associa ed wi h bo h clima ic a iables
and he popula ion pheno ypic measu emen s, while con olling o e ec s o de-
mog aphy and popula ion s uc u e. The op candida e SNPs we e en iched on he
X and ou h ch omosomes, and hey also lay nea genes implica ed in o he s udies
o cold ole ance and popula ion di e gence in his species and i s close ela i es. We
conclude ha ecological adap a ion has con ibu ed o he di e gence o D. mon ana
popula ions h oughou he genome and in pa icula on he X and ou h ch omo-
somes, which also showed highes in e popula ion FST. This s udy demons a es ha
ecological selec ion can d i e genomic di e gence a di e en scales, om candida e
genes o ch omosome- wide e ec s.
KEYWORDS
chill coma eco e y ime, cline popula ions, cold ole ance, CTmin, D. mon ana, en i onmen al
adap a ion, genomic di e gence
3784
|
WIBERG E al.
1 | INTRODUCTION
The geog aphic s uc u e o a species is a esul o i s phylogeo-
g aphic his o y, in luenced by pas and p esen dispe sal, popula ion
demog aphy, and selec ion. Ob aining genome- wide da a on gene ic
polymo phisms ac oss mul iple popula ions o a species is becom-
ing ela i ely easy, bu in e p e ing he pa e ns o geog aphic a i-
a ion in such da a and iden i ying genes which a y p ima ily due
o selec ion emains challenging. O en simple”ou lie ” app oaches
using genome scans which measu e gene ic di e en ia ion such as
FST o Dxy a e adop ed, bu esul s a e di icul o in e p e due o
con ounding e ec s o selec ion, d i and popula ion s uc u e, o
genomic ea u es such as in e sions and o he causes o a ia ion
in ecombina ion a e (C uickshank & Hahn, 2014; Noo & Benne ,
2009; Ra ine e al., 2017; Wol & Elleg en, 2016). I en i onmen al
da a a e a ailable, we can use associa ions wi h such ac o s o help
iden i y loci whe e gene ic di e en ia ion co a ies wi h his en i-
onmen al a ia ion. Some genome scan me hods can inco po a e
en i onmen al a ia ion and simul aneously i e ec s o co a-
iance wi h en i onmen al ac o s, while con olling o e ec s o
popula ion demog aphy (Foll & Gaggio i, 2008; de Villeme euil &
Gaggio i, 2015). This app oach has success ully iden i ied gene ic
a ia ion associa ed wi h al i ude in humans, among o he examples
(Foll e al., 2014; Gau ie , 2015; de Villeme euil & Gaggio i, 2015)
and has become a use ul app oach o in es iga e he ecological ad-
ap a ions unde lying popula ion di e gence.
Clinal pa e ns o a ia ion in pheno ypes o gene equencies
ha e a long his o y o being used o in e selec ion along eco ones,
and analyses o cline shape can some imes iden i y loci unde di-
ec selec ion om o he s showing clinal a ia ion o o he easons,
such as phylogeog aphic his o y (Ba on & Gale, 1993). Such s udies
can be e y powe ul, especially when independen pa allel clines
a e a ailable. Fo example, Kolaczkowski e al., (2011) sampled iso e-
male lines om ex emes o a cline in Aus alian popula ions o D.
melanogas e , and ound many genes implica ed in clinally a ying
pheno ypes o show highes di e en ia ion. Also, Be gland e al.,
(2014), and Kapun e al., (2016) sampled No h Ame ican clines in
D. melanogas e and D. simulans o e se e al yea s o unco e clinal
a ia ion a a genome scale. Be gland e al., (2014) also ound con-
sis en luc ua ions in allele equencies o a popula ion sampled
o e se e al seasons, indica ing a egula esponse o seasonally
a ying selec ion p essu es. On he o he hand, Machado e al.,
(2015) concluded ha mig a ion and gene low play a g ea e ole
han adap a ion in he o e all clinali y o genomic a ian s in D. sim-
ulans han D. melanogas e . While he wo species sha e a signi ican
p opo ion o he genes showing clinal a ia ion, hei di e ences
in o e win e ing abili y, mig a ion and popula ion bo lenecks p ob-
ably ac as addi ional d i e s o di e ences in pa e ns o a ia ion
be ween hem (Machado e al., 2015). Simila s udies o clinal a ia-
ion in pheno ypes and allele equencies ha e also been ca ied ou
in o he insec s (Paolucci e al., 2016), plan s (B adbu y e al., 2013;
Chen e al., 2012), mammals (Ca nei o e al., 2013; Hoeks a e al.,
2004), ish (Vines e al., 2016), and o he o ganisms (Endle , 1973,
1977; Takahashi, 2015). Howe e , mo e s udies a e needed ha si-
mul aneously compa e pa e ns o a ia ion in allele equencies and
po en ially causal en i onmen al a ia ion o ecologically impo an
ai s while con olling o popula ion s uc u e, as such s udies a e
necessa y o de e mine i adap a ion o clima e is di ec ly d i ing
pa e ns o gene ic di e en ia ion.
He e, we in es iga ed geog aphic a ia ion a bo h he pheno-
ypic and gene ic le els in D osophila mon ana samples om wo
con inen s. This species has sp ead a ound he no he n hemisphe e
(Th ockmo on, 1969) and bo h m DNA and mic osa elli e da a ha e
e ealed gene ically dis inc Finnish and No h Ame ican popula ions
(Mi ol e al., 2007). Mo eo e , mo e ecen modelling o genome-
wide SNP equencies sugges ha he Finnish- No h Ame ican spli
happened a ound 1.75 Mya and ha o he No h Ame ican pop-
ula ions sho ly a e ha (Ga lo sky e al., 2020). I is one o he
mos cold- ole an D osophila species (Kelle mann e al., 2012) and
he basic cold ole ance o D. mon ana lies can inc ease owa ds he
cold seasons h ough wo mechanisms, pho ope iodic ep oduc i e
diapause (Vesala & Hoikkala, 2011) and cold- acclima ion induced
by a dec ease in day leng h and/o empe a u e (Kau anen e al.,
2019; Vesala, Salminen, Kos al, e al., 2012; Vesala, Salminen, Laiho,
e al., 2012). D. mon ana popula ions ha e been ound o show clinal
a ia ion in he c i ical day leng h equi ed o diapause induc ion
(CDL; Lankinen e al., 2013; Tyukmae a e al., 2011). The e is also
a co ela ion be ween CDL and la i udinally co a ying clima ic ac-
o s such as he mean empe a u e o he coldes mon h (Tyukmae a
e al., 2020). In addi ion, D. mon ana popula ions om di e en geo-
g aphic egions show a ia ion in hei cou ship cues and ma e
choice (Klappe e al., 2007; Rou u e al., 2007), which has led o
pa ial ep oduc i e isola ion be ween some dis an popula ions
(Jennings e al., 2011, 2014). A he gene ic le el, di e en ial gene
exp ession s udies ha e iden i ied candida e genes unde lying dia-
pause (Kanka e e al., 2010, 2016), pe cep ion o day leng h (Pa ke
e al., 2016), and cold acclima ion (Pa ke e al., 2015). Fu he mo e,
a quasi- na u al selec ion expe imen o sho e CDL, accompa-
nied by a dec ease in cold- ole ance, induced widesp ead changes
in loci wi h po en ial oles wi h hese ai s (Kau anen e al., 2019).
Finally, popula ion genomic analyses ha e iden i ied se e al ou lie
loci when examining di e en ia ion be ween No h Ame ican and
Eu opean popula ions (Pa ke e al., 2018). All his makes D. mon ana
an in e es ing example o nascen specia ion, po en ially in luenced
by ecological adap a ion.
He e, we sough o ask o wha ex en pa e ns in he genomic
di e gence o D. mon ana popula ions ac oss wo con inen s a e
co ela ed wi h clima ic a ia ion and pheno ypic esponses o
cold adap a ion. We pe o med pooled whole- genome sequencing
(pool- seq) on six di e en popula ions and used Bayesian me h-
ods o examine he associa ion be ween genomic di e en ia ion
be ween popula ions and en i onmen al a iables ac oss bo h
con inen s. We also pheno yped popula ions o wo di e en
cold ole ance measu es; c i ical he mal minimum (CTmin), and
chill coma eco e y ime (CCRT), and in es iga ed he associa ions
be ween hem and he gene ic and clima ic da a. Ul ima ely, we
|
3785
WIBERG E al.
asked i he genomic loci showing an associa ion be ween ge-
ne ic and en i onmen al di e en ia ion also showed associa ion
wi h popula ion di e en ia ion in cold ole ance pheno ypes, and
examined he possible o e lap be ween he se o genes close o
candida e SNPs wi h se s o candida e genes om p e ious s udies
o cold adap a ion in D. mon ana. I popula ion di e en ia ion is
d i en by ecological selec ion hen we would p edic he ex eme
cold adap a ion o D. mon ana o ha e le a signa u e o genomic
di e gence associa ed wi h en i onmen al and pheno ypic di e -
en ia ion ac oss hese loci.
2 | MATERIALS AND METHODS
2.1  |  Sample collec ions and DNA ex ac ion
We collec ed samples o 49– 50 wild- caugh lies om six D. mon-
ana popula ions om a ange o la i udes om 66°N o 38°N in he
sp ing o 2013 o 2014. Fou o hese popula ions ep esen ed a
ange o la i udes in No h Ame ica (N.A.), and wo popula ions we e
om a ange in Finland (Figu e 1; Table 1). Samples o wild- caugh
lies om he six popula ions we e s o ed in e hanol ( he male/ e-
male a io a ied ac oss samples; Table 1) and DNA o indi idual
lies was ex ac ed using CTAB solu ion and phenol- chlo o o m-
isoamylalcohol pu i ica ions in 2016. Genomic DNA was ex ac ed
om indi idual lies and quan i ied using Qubi (The mo Fishe
Scien i ic), and an equal amoun o DNA om each indi idual (50 ng)
was pooled in o he inal sample. Sequencing was pe o med a he
Finnish Func ional Genomics Cen e in Tu ku, Finland (www.b k. i/
unc ional - genomics) on he Illumina HiSeq3000 pla o m (pai ed-
end eads, ead leng h = 150 bp, es ima ed co e age ~121x).
2.2  |  Pheno yping
We measu ed he c i ical he mal minimum (CTmin) and chill coma e-
co e y ime (CCRT) o lies om six popula ions. Fly samples o hese
es s we e collec ed o i e popula ions (Sewa d, Te ace, Ash o d,
C es ed Bu e, and Ko pilah i) om mass- b ed popula ions ha ha e
been main ained in he labo a o y since 2013– 2014. Fo he Oulanka
popula ion, lies we e collec ed om h ee iso emale s ains (es ab-
lished in 2014), because he mass- b ed popula ion had been con ami-
na ed by ano he species. The mass- b ed popula ions we e o iginally
es ablished om F4 p ogenies o 20 iso emale s ains each (400 lies)
and ha e been main ained in cons an empe a u e (19 ± 1℃) and
ligh (24 h o ligh , o p e en he lies om en e ing diapause) e-
gimes o abou 20– 25 gene a ions be o e he expe imen . All lies
we e supplied wi h esh mal medium in hal - pin bo les e e y week
(Lako aa a, 1969). The newly eme ged lies we e collec ed using ligh
CO2 anaes hesia wi hin 24 h a e eme gence, sepa a ed by sex and
placed in mal - ials in he same condi ions un il sexual ma u i y (20–
21 days old) and we e hen used in CTmin and CCRT es s. The same
indi idual lies we e i s assayed o CTmin and hen o CCRT, he
lies we e no anaes he ised be o e hese es s.
We assayed a o al o 328 emales and 302 males o CTmin
and CCRT. These assays we e done in ba ches o be ween 22 and
30 lies, and spli e enly by sex (21 ba ches in o al). Be ween 32
and 46 (mean 39) lies pe popula ion pe sex we e es ed ( o he
FIGURE 1 Maps o all D. mon ana popula ions in his s udy. Panels show popula ion loca ions om (a) Finland and (b) No h Ame ica
showing he loca ions o all popula ions sampled. Labelled, blue ci cles gi e he loca ions o popula ions sampled o pheno yping and
sequencing
(a) (b)
3786
|
WIBERG E al.
Oulanka iso emale s ains, be ween 15 and 39 lies pe s ain pe
sex we e used). CTmin es s a e based on de ec ing he empe a-
u e (CTmin) a which lies lose neu omuscula unc ion and en e
e e sible s a e o chill coma (Ande sen e al., 2015). In hese es s,
he lies we e placed in o ubes sealed wi h pa a ilm and subme ged
in o a 30% glycol- wa e mix u e in Julabo F32- HL chambe . The
empe a u e was dec eased a he a e o 0.5℃ pe min ( om 19℃
o – 6 ℃) and CTmin was de e mined by eye, as he empe a u e a
which a ly was unable o s and on i s legs. Immedia ely a e he
CTmin es , he empe a u e was se o – 6℃ and he lies we e le
in his empe a u e o 16 h. Vials we e hen quickly aken ou o he
glycol- wa e ba h and he lies’ CCRT was de e mined as he ime
equi ed o he lies o eco e om chill coma and s and on hei
legs. The ambien oom empe a u es was eco ded du ing ials bu
in an ini ial analysis including sou ce popula ion and oom empe a-
u e he e was no s a is ically signi ican e ec o oom empe a u e
on CCRT (F5,569 = 0.72, p = .4), his a iable was he e o e le ou o
u he analyses.
To in es iga e popula ion di e ences in CTmin and CCRT phe-
no ypes we i gene al addi i e mixed models (GAMMs) in ( .
3.6.3; R De elopmen Co e Team, 2020) using he “mgc ” pack-
age (Wood, 2004). In simple linea models including popula ion, sex,
and expe imen al ba ch as a ixed e ec s, expe imen al ba ch had
an e ec on CTmin (F20,533 = 2.8, p < .01), al hough compa ed o
popula ion (F5,533 = 12, p < .01), and sex (F1,533 = 11.5, p < .01) his
e ec was no s ong. Expe imen al ba ch had a s onge e ec on
CCRT (F20,533 = 2.83, p < .01), compa ed o popula ion (F5,533 = 4.5,
p < .01) and sex (F1,533 = 1.5, p < .22). We he e o e included i as
a andom e ec in GAMM analyses. The ull GAMMs included al-
i ude, la i ude, and sex as ixed e ec s and expe imen al ba ch
as a andom e ec . We used a cubic eg ession spline as he basis
smoo hing unc ion o bo h al i ude and la i ude. The aw da a o
all pheno yping a e gi en in Table S1. The ull models a e shown in
he Suppo ing In o ma ion (Supplemne a y Ma e ial).
2.3  |  Bioclima ic a iables and
popula ion geog aphy
We ob ained ep esen a i e clima e da a om he Wo ldClim da a-
base (Hijmans e al., 2005) o each D. mon ana popula ion sampled
o he pool- seq (see abo e), as well as o 18 addi ional popula ions
o his species (Table S1; Tyukmae a e al., 2020). We downloaded
he clima e da a and ex ac ed he alues co esponding o popula-
ion coo dina es using he R package “ as e ” ( e sion 2.5- 8; Hijmans
e al., 2016). In o al his amoun s o 55 bioclima ic a iables o each
popula ion (Table S1). To educe he numbe o a iables in he da a
se a p inciple componen s analysis (PCA) was pe o med using he
“PCA()” unc ion om he R package “Fac oMineR” ( e sion 1.28;
Lê e al., 2008). P inciple componen s we e kep o u he analysis
i hei eigen alues we e >1. PCA sco es o each popula ion we e
z- ans o med using he “scale()” unc ion in base R. Addi ionally,
CTmin and CCRT we e summa ised o a mean alue o each popula-
ion. In o al, his gi es ou “en i onmen al” a iables measu ed o
each popula ion (PC1, PC2, CTmin, and CCRT).
2.4  |  Mapping, SNP calling and genomic analysis
Quali y o aw eads was checked wi h as qc ( . 0.11.5; And ews,
2015) and eads we e immed using immoma ic ( . 0.32; Bolge
e al., 2014; see Suppo ing In o ma ion (Supplemen a y Ma e ial)
o ull imming pa ame e s). T immed eads we e mapped o he
D. mon ana e e ence genome (Pa ke e al., 2018) using bwa mem
( . 0.7.7; Li, 2013) wi h he de aul op ions bu keeping only align-
men s wi h a mapping quali y o >20 ollowing bes p ac ice guide-
lines o pool- seq (Schlö e e e al., 2014). Duplica e alignmen s
we e emo ed wi h sam ools mdup ( 1.3.1; Li e al., 2009) and
egions a ound indels we e ealigned using pica d ( e sion 1.118;
B oad Ins i u e), ga k ( . 3.2- 2; McKenna e al., 2010) and sam ools.
Sepa a e.bam iles o each o he sequenced popula ions we e i-
nally me ged using bam ools ( . 2.4.0; Ba ne e al., 2011).
O e 80% o eads we e p ope ly mapped and e ained in all sam-
ples. The mean co e age o Sewa d samples was nea ly wice ha
o he o he samples (Figu e S1 and Figu e S2). To emo e he po en-
ial o his di e ence o cause a e ac s in downs eam analyses,
he.bam iles o Sewa d we e downsampled o con ain 94.1 million
eads ( he a e age ac oss he i e emaining popula ions). Median
empi ical co e age was be ween ~62 and 88x (Table S2; Figu e S3)
and much less a iable among he popula ions, allowing common
maximum and minimum h esholds o be se based on he agg ega e
dis ibu ion. Allele coun s o each popula ion a each genomic po-
si ion we e ob ained wi h sam ools mpileup ( e sion 1.3.1; Li e al.,
TABLE 1 The sou ces o genomic samples (coo dina es and he
name o he nea es own), al i ude o he sampling si e, he yea
in which sampling was pe o med, and he numbe o males and
emales sampled (M/F) o each pool
Sou ce Sampling si e Yea M/F
USA, Alaska Sewa d
60°9’N; 149°27′W
Al i ude 35 m
2013 30/20
Canada, B i ish Columbia Te ace
54°27′N; 128°34′W
Al i ude 217 m
2014 22/27
USA, Washing on Ash o d,
46°45′N; 121°57′W
Al i ude 573 m
2013 16/34
USA, Colo ado C es ed Bu e
38°54′N; 106°57′W
Al i ude 2,900 m
2013 36/13
Finland Oulanka
66°22′N; 29°20′E
Al i ude 337 m
2013 25/25
Finland Ko pilah i
62°20′N; 25°34′E
Al i ude 133 m
2013 27/23

|
3787
WIBERG E al.
2009) using op ions o skip indel calling as well as igno ing eads
wi h a mapping quali y <20 and si es wi h a base quali y <15. This
was ollowed by he heu is ic SNP calling so wa e PoolSNP using a
minimum coun o 5 o call an allele, and a minimum co e age o 37
and a maximum co e age <95 h pe cen ile o he sca old- wide co -
e age dis ibu ion o call a SNP (Kapun e al., 2020). E en i all hese
il e s we e passed, an allele was no conside ed i i s equency
was <0.001. Finally, we only conside ed SNPs on sca olds >10 kb in
leng h. The inal se consis ed o 2,190,511 biallelic SNPs ha could
be placed on sca olds o de ed along he ch omosomes and we e
used in downs eam analyses.
To es o an associa ion be ween he ou en i onmen al a i-
ables and gene ic di e en ia ion we used bayescen ( . 1.1; Foll &
Gaggio i, 2008; de Villeme euil & Gaggio i, 2015). bayescen i s a
model o FST o popula ion di e en ia ion o each locus, inco po-
a ing en i onmen al di e en ia ion as a p edic o (included as he
pa ame e g) while i ing wo locus- speci ic e ec s, one o en i-
onmen al selec ion, and he o he o o he p ocesses (demog a-
phy, o he ypes o selec ion). I he e o e con ols o con ounding
e ec s o popula ion s uc u e/ ela edness in es ing o an associ-
a ion wi h en i onmen al a iables. bayescen was un wi h i e pilo
uns o 1000 i e a ions each, ollowed by a main chain o 4000 i e -
a ions o which 2000 we e disca ded as bu nin. Fou MCMC chains
we e un o each analysis o e alua e con e gence o he chains o
common pa ame e es ima es. Because o he unbalanced numbe
o males (and he esul ing a ia ion in a ios o X:Y ch omosomes)
in he pools, bayescen analyses we e pe o med sepa a ely on SNPs
ha could be assigned o he au osomal linkage g oups (ch omo-
somes) and he X ch omosome. Raw coun da a we e used o he
au osomal da a. Fo X linked SNPs, allele coun da a we e scaled
o he known numbe o X ch omosomes in he pool using ne , he
e ec i e sample size aking in o accoun he mul iple ounds o bi-
nomial sampling inhe en o a pool- seq design (Fede e al., 2012;
Kolaczkowski e al., 2011). Ne scales he allele coun s a each SNP
downwa ds based on he known numbe o ch omosomes in he
pool (see Table 1).
Chains we e assessed o con e gence wi h he “coda” package
( . 0.19- 1; Plumme e al., 2006). Con e gence was eached ac oss
he ou chains o mos analyses and pa ame e s (po en ial scale
educ ion ac o s (PSRFs) o ~1 in a Gelman- Rubin diagnos ic es ;
Figu es S4– S7), excep o analyses o au osomal SNPs and PC2 as
he en i onmen al a iable which showed mild con e gence p ob-
lems (PSRF = 1.71), al hough pa ame e es ima es ag eed well wi h
he o he chains. Thus, his i s chain (black line in Figu es S4– S7)
was discoun ed o all pa ame e s and only es ima es om he e-
maining h ee chains we e used. The union o signi ican SNPs (using
q- alues o he g pa ame e , desc ibing he associa ion be ween ge-
ne ic di e en ia ion a a locus and en i onmen al di e en ia ion, o
con ol he FDR a 0.05, i.e., SNPs wi h q- alues <0.05) ac oss hese
chains we e aken as he inal candida e SNPs.
Finally, popula ion gene ic s a is ics (π, and Tajima's D) we e
compu ed in windows o 10 kb wi h a s ep size o 5 kb using me h-
ods implemen ed o pool- seq da a (Kapun e al., 2020). Windows
con ained a mean o 488.2 ± 1.2 SNPs (mean ± SE). These s a is ics
we e only compu ed o sca olds wi h a leng h >50 kb. SNP- wise
FST was compu ed o each popula ion wi h he package “pool -
s a ” ( . 1.1.1; Hi e e al., 2018) by i s compu ing all pai wise
alues, and hen de i ing popula ion speci ic FST alues by a e -
aging ac oss all pai wise alues whe e a popula ion was included.
We compu ed 95% con idence in e als o mean SNP- wise FST on
each ch omosome om he dis ibu ion o mean FST alues ac oss
100 boo s ap samples o SNPs. See he Supplemen a y Ma e ial
(Suppo ing In o ma ion) o pseudocode commands o he key
pipeline s eps.
We pe o med GO e m en ichmen analyses wi h da id ( .
6.8; Huang e al., 2009a, 2009b) and gowinda ( . 1.12; Ko le &
Schlö e e , 2012), which accoun s o di e ences in gene leng h. Fo
da id, since he D. mon ana anno a ion con ains in o ma ion abou
o hologs in D. i ilis, we ex ac ed all genes wi hin 10 kb ups eam
o downs eam o a candida e SNP wi h an o holog in D. i ilis and
submi ed hem o da id. Fo analyses wi h gowinda, because he
gene- se s need o be gi en manually, we ob ained gene- se s om
FuncAssocia e3 (Be iz e al., 2009). As abo e, we conside ed SNPs
wi hin 10 kb o a gene (- - gene- de ini ion updowns eam 10,000) and
pe o med 1 million simula ions o ob ain he empi ical p- alue dis-
ibu ion. Regula o y elemen s, such as enhance s, and ansc ip ion
ac o binding si es, can occu up o 1 Mb up- o downs eam om a
a ge gene in o he species (Chan e al., 2010; Mas on e al., 2006;
Pennacchio e al., 2013; We ne e al., 2010) bu gene ally lie wi hin
2 kb o a gene egion in D. melanogas e (A nos i, 2003), hus 10 kb
ep esen s a comp omise.
3 | RESULTS
3.1  |  Cold ole ance measu es, bioclima ic a iables
and popula ion geog aphy
Ac oss indi iduals, he e was no e idence o an associa ion be ween
minimum c i ical empe a u e (CTmin; pooled ac oss sexes) and he
chill coma eco e y ime (CCRT; co = 0.08, p = .07). Nei he was he e
e idence o an associa ion be ween hese ai s ac oss popula ions
(co = 0.22 p = .68). CTmin was signi ican ly di e en be ween sexes,
wi h emales being mo e cold ole an ( = 3.14, p = .002). While la i-
ude appea s o ha e a nonlinea e ec on CTmin (Figu e 2a and
Figu e S8A), he e ec i e deg ees o eedom (ed = 1) o he pa ial
e ec o la i ude sugges ed a la gely linea e ec a e accoun ing
o al i ude (F = 40.4, p < .001). Al i ude had a mo e complex e ec
(ed = 1.262, F = 24.0, p < .001), bu wi h a s ong al i ude ou lie in
C es ed Bu e (Table 1), his esul should be aken wi h some cau-
ion. Mo eo e , he o e all adjus ed R- squa ed was low a only 0.09.
Fo CCRT, only la i ude had a signi ican e ec (ed = 1, F = 17.2,
p < .001; Figu e 2b). Taken oge he , CTmin was lowe o popula-
ions a highe la i udes as well as o popula ions a highe al i udes
(Figu e 2a and Figu e S8A). As one would expec , CTmin showed on
a e age lowe alues (Figu e 2a and Figu e S8A), and CCRTs we e
3788
|
WIBERG E al.
sho e (Figu e 2b and Figu e S8B), a highe la i udes, meaning ha
mo e no he n and highe al i ude popula ions show highe cold
ole ance.
We pe o med a p incipal componen analysis (PCA) o he
Wo ldClim clima e da a o a o al o 24 D. mon ana popula ions,
om which we had collec ed samples and whe e clima e da a we e
a ailable. We wan ed o examine whe he he popula ions chosen
o he cold ole ance es s and genome sequencing we e ep esen-
a i e o he ange o en i onmen al a ia ion expe ienced by he
species. These analyses iden i ied ou p incipal componen s (PCs)
ha oge he explained abou 98% o he a ia ion (Figu e 3a and
Figu e S9). The i s wo PCs sepa a ed he popula ions oughly by
a measu e o dis ance om he coas (PC1) and hen by la i ude and
al i ude (PC2). PC1 explained ~54% o he a ia ion (Figu e 3a and b)
and loaded hea ily on clima e and biological a iables associa ed wi h
p ecipi a ion and empe a u e such as “mean empe a u e o coldes
FIGURE 2 Popula ion pheno ypes. (a) and (b) show he a ia ion
in CTmin and CCRT ac oss popula ions and la i ude, espec i ely.
Solid and dashed lines show he p edic ed alues om he pa ial
e ec o la i ude om he bes model (see Resul s) o males
and emales, espec i ely. In (b) al hough he poin s a e plo ed
sepa a ely o males and emales, he bes model only included
la i ude as a co a ia e. In bo h panels he poin s and e o ba s gi e
he mean and s anda d e o s, espec i ely. Figu e S8 shows he
ull da a o each popula ion
(a)
(b)
FIGURE 3 P inciple componen analysis (PCA). (a) The
dis ibu ion o all popula ions along he wo i s PC axes. The
loadings o a iables on each axis can be ound in Table S3. Blue
ci cles gi e he popula ions ha ha e been pool- sequenced o his
s udy, ed squa es and iangles gi e he o he Finnish and No h
Ame ican popula ions, espec i ely. (b) and (c) gi e PCs 1 and 2
as a unc ion o la i ude and al i ude, espec i ely. The legend o
al i ude in bo h (b) and (c) is gi en in (c)
(a)
(b)
(c)
|
3789
WIBERG E al.
qua e ”, “p ecipi a ion o we es mon h”, “annual p ecipi a ion” (see
Figu e S9 and Table S3). Meanwhile, PC2 explained ~23% o he a i-
a ion (Figu e 3a and c) and loaded hea ily on biological a iables ha
a e associa ed wi h la i udinal clinali y, e.g., “mean diu nal [ empe a-
u e] ange,” and “iso he mali y” which is he diu nal ange di ided
by he mean “annual [ empe a u e] ange”. The emaining PCs (PC3
and PC4) explained abou 11.5% and 5% o he a ia ion, espec-
i ely and did no cap u e as much o he clima ic a ia ion, and we e
he e o e no conside ed u he . La i ude (Spea man's ank co ela-
ion: ho = – 0.50, p = .01) bu no al i ude ( ho = – 0.19, p = .39) co -
ela ed signi ican ly wi h PC1. Howe e , bo h al i ude ( ho = – 0.50,
p = 0.02) and la i ude ( ho = – 0.59, p = .003) co ela ed wi h PC2
(see also Figu e 3c). Impo an ly, hese pa e ns we e ai ly obus
also when pe o ming PCA using only he six popula ions o which
genomic da a we e collec ed. Loadings on PC1 we e highly compa-
able, and al hough co ela ions wi h la i ude and al i ude did no
achie e s a is ical signi icance, hey we e simila in magni ude and
di ec ion (la i ude s. PC1: ho = 0.43, al i ude s. PC1: ho = – 0.14,
p = .8). Fo PC2 he co ela ions we e also no s a is ically signi ican
and di e ed o la i ude bo h in magni ude and di ec ion ( ho = 0.6,
p = .2), while o al i ude, he co ela ion emained simila and ma -
ginally non- signi ican ( ho = – 0.82, p = .06). Mo eo e , while he
ange o he sequenced popula ions co e ed ~70% o he ange o all
popula ions o PC1, hey co e ed only ~50% o he ange o PC2
(Figu e 3a). This sugges s ha en i onmen al a ia ion ac oss he six
popula ions selec ed o sequencing e lec s ha expe ienced by all
24 popula ions a leas o he PC1 axis. The e o e, he ela ionship
be ween en i onmen al a iables and gene ic di e en ia ion in he
samples selec ed o pool- seq is likely o e lec ue pa e ns ac oss
popula ions o D. mon ana, hough powe migh be somewha e-
duced o PC2. Fo all subsequen analyses, we used he esul s o
he PCA using all 24 popula ions, and o examine he associa ion
be ween clima e and pheno ype, we compa ed hese ac oss he six
popula ions. CTmin was posi i ely co ela ed wi h PC1 (Pea son's
co ela ion coe icien (co ) = 0.94, p < .01) and had a ma ginally non-
signi ican associa ion wi h PC2 (co = 0.75, p = .09). Howe e , CCRT
showed no ela ionship wi h ei he PC1 (co = 0.09, p = .87) o PC2
(co = – 0.46, p = .36) al hough he small sample sizes (N = 6 in all
cases) make eliable conclusions di icul (see Figu e S10).
3.2  |  Genomic analyses
The numbe o SNPs wi h a signi ican (q- alue <0.05) associa ion be-
ween popula ion- based FST, he wo cold ole ance measu es (CCRT
and CTmin), and he wo PCs o he bioclima ic da a a ied om 312
(ch omosome 3, CCRT) o 2480 (ch omosome 4, PC2) ac oss he
ch omosomes (Figu e 4). Using PC1 as an en i onmen al a iable
wi h bayescen ga e a o al o 2976 and 1528 SNPs wi h a q- alue
<0.05 on he au osomes and on he X ch omosome, espec i ely.
In e es ingly, he dis ibu ion ac oss he ch omosomes was no an-
dom. Using he dis ibu ion o all SNPs o ob ain expec ed coun s,
he e was a signi ican de ia ion om expec a ion (Χ2 = 2906.4,
d. . = 4, p < .001). The e we e many mo e SNPs han expec ed on
ch omosome 4 (1432 s. 954) and on he X ch omosome (1528 s.
526; Figu e 4). Resul s we e simila o PC2 wi h 6607 and 1861
SNPs wi h a q- alue <0.05 on he au osomes and X ch omosomes,
espec i ely. Again, he e was a signi ican de ia ion om he ex-
pec ed dis ibu ion o SNPs ac oss he ch omosomes (Χ2 = 1681.9,
d. . = 4, p < .001) wi h an o e ep esen a ion on he ou h (2480 s.
1794) and he X ch omosomes (1861 s. 989; Figu e 4).
We also used a e age CTmin alues pe popula ion in simila anal-
yses and ound a o al o 2668 and 1272 SNPs wi h a q- alue <0.05
on he au osomes and X ch omosomes, espec i ely. The pa e n
o signi ican de ia ions om expec ed dis ibu ions (Χ2 = 2,526.2,
d. . = 4, p < .001) was also due o an excess on he ou h (1383
s. 835) and he X ch omosomes (1272 s. 460; Figu e 4). Simila
esul s we e ound o CCRT wi h a o al o 2240 and 1228 SNPs
wi h a q- alue <0.05, espec i ely. Once again, he e was a signi i-
can de ia ion om he expec ed dis ibu ion o SNPs (Χ2 = 2825.6,
d. . = 3, p <.001) wi h and excess on he ou h (1252 s. 735) and he
X ch omosomes (1228 s. 405; Figu e 4). Manha an plo s o he dis-
ibu ion o SNPs ac oss ch omosomes a e gi en in Figu es S11– S14.
To mo e closely examine he loci implica ed in he ou BayScEn
analyses, we iden i ied genes wi hin 10 kb o , o con aining he can-
dida e SNPs (Table S4). O e all, he e is qui e a la ge o e lap among
hese genes wi h ~39% (1102 in o al) o hem being sha ed by all he
ou analyses (Table S4, Figu e S15). This is a mo e han would be
expec ed (mean ± SE: 3.6 ± 0.06) by andomly esampling (wi hou
eplacemen ) a numbe o genes equi alen o ha associa ed wi h
op SNPs o each BayScEn analysis om he D. mon ana anno a-
ion, hen compu ing he ou - way in e sec ion. This la ge obse ed
o e lap he e o e p obably e lec s a genuine o e lap among he
unde lying op SNPs (da a no shown). Some (147, ~13%) o hese
common genes a e no el o D. mon ana (i.e., no anno a ed in D. i-
ilis o in o he D osophila spp.) and he e o e ha e no anno a ion,
bu 955 ha e an iden i iable D. i ilis o holog (Table S4). The sec-
ond la ges se a e hose genes unique o he analysis o PC2 as
an en i onmen al co a ia e (Figu e S15). We pe o med unc ional
en ichmen analyses wi h da id keeping all anno a ion clus e s wi h
an en ichmen sco e >1.3 (co esponding o an a e age co ec ed
p- alue o .05). This e ealed se e al common ca ego ies o genes
associa ed wi h he clima ic a iables and popula ion pheno ypes
(Table S5). Fo example, e ms associa ed wi h memb ane and
ansmemb ane s uc u es, immunoglobulins, HAD hyd olase and
nucleo ide binding we e en iched in mos o he a iables (Table
S5). In e es ingly, he e we e also se e al gene on ology ca ego ies
ha we e only en iched in one o he a iables, such as glycoside
and ATPase hyd olase in CCRT and ion channels and anspo , as
well as me al binding in PC1 (Table S5). gowinda analyses e ealed
signi ican en ichmen o GO e ms (a e accoun ing o a ia ion
in gene leng hs and co ec ing o mul iple es ing) only o genes
nea SNPs associa ed wi h a ia ion in CCRT ac oss popula ions.
In e es ingly, he e m “ca bohyd a e de i a i e binding” was iden i-
ied among he mos en iched e ms (Table S6), which ag ees closely
wi h some e ms iden i ied using da id o he same a iable (Table
3790
|
WIBERG E al.
S5). Simila ly, o PC2, “nucleo ide binding” was among he op sco -
ing en iched e ms in bo h gowinda and da id analyses, hough his
was no signi ican a e con olling o mul iple es ing in gowinda
(Table S5, Table S6).
We hen compa ed loci nea candida e SNPs wi h genes impli-
ca ed in p e ious s udies o clima ic adap a ion in D. mon ana in-
cluding gene exp ession s udies o ai s connec ed o diapause
and cold- ole ance (Kanka e e al., 2010, 2016; Pa ke e al., 2015,
2016). Addi ionally, se e al candida e genes ha e been iden i ied
nea he mos signi ican ly di e en ia ed SNPs among D. mon ana
popula ions om Oulanka (Finland), and om No h Ame ican pop-
ula ions in Colo ado and Vancou e (Pa ke e al., 2018). Finally,
quasi- na u al selec ion expe imen s iden i ied se e al genes wi hin
10 kb o SNPs esponding o selec ion o a sho e CDL o diapause
induc ion (Kau anen e al., 2019). We es ed o an o e lap be ween
he o al se o genes wi hin 10kb o ou lie SNPs om all o he
BayScEn uns (N = 2694) and he candida e gene se s iden i ied
in ea lie s udies (see Table S7 o he gene se s and s udies used).
We compu ed a boo s ap dis ibu ion o o e laps by sampling 2694
andom genes om he D. mon ana anno a ion (N = 13,683). Fo
each o he gene se s om p e ious s udies his was done 100 imes
and he dis ibu ion compa ed o he empi ical o e lap. Resul s a e
gi en in Table 2. In all he cases he empi ical o e lap was g ea e
han expec ed by chance wi h empi ical p- alues <.05 ( anging om
<.0001 o .01; Table 2). The only gene ha was ound in all i e o he
p e ious s udies used and in he compa ison he e is called sides ep
II (side- II; Table S8). Un o una ely, he e is no in o ma ion a ailable
abou he biological p ocesses o molecula unc ions connec ed o
i . Mo eo e , om 44 o he genes ha we e common o ou o ou
p e ious s udies and his s udy (Table S8) mos (27) ha e an o holog
in D. melanogas e . These genes ha e molecula unc ions such as
ansmemb ane signaling o anspo e , ace ylcholines e ase, ATP
binding, p o ein se ine/ h eonine kinase, ca boxylic es e hyd olase,
o Rho guanyl- nucleo ide exchange ac o ac i i y (Thu mond e al.,
2019). Many o he genes a e also connec ed o me al ion, nucleid
acid o zinc ion binding (Table S8). A e iden i ying in o ma ion on
molecula o biological unc ion and In e p o domains, e en ually
only i e genes emained o which he e was no in o ma ion a ail-
able (Table S8).
Examina ion o popula ion gene ic pa ame e s iden i ied he
C es ed Bu e popula ion as anomalous. The dis ibu ion o Tajima's
D is cen ed close o ze o in mos popula ions, being sligh ly mo e
nega i e in No h Ame ican popula ions (Figu e S16A). Howe e ,
C es ed Bu e is an ou lie wi h a g ea ly educed genome- wide
Tajima's D (Figu e S16A). Simila ly, di e si y (π) is also lowe in his
popula ion han in o he popula ions. The e is no o e all ela ion-
ship be ween la i ude and π (Spea man's ho = 0.14, N = 6, p = .8;
Figu e S16B) bu he e is a s ong co ela ion be ween la i ude and
Tajima's D which is in luenced by his popula ion (wi h C es ed Bu e:
ho = 0.88, N = 6, p = .03; wi hou C es ed Bu e: ho = 0.8, N = 5,
p = .13). Al hough C es ed Bu e occu s a a much highe al i ude
(>2800 m) han o he popula ions nei he Tajima's D no π co ela ed
signi ican ly wi h al i ude (Tajima's D: ho = – 0.6, N = 6, p = .24, π:
ho = – 0.6, N = 6, p = .24). Fu he mo e, FST was simila ac oss all
popula ions and ch omosomes wi h he excep ion o C es ed Bu e
which emained an ou lie wi h high FST (Figu e 5). Boo s apped
95% con idence in e als a e p esen ed o gi e a guide o s a is i-
cal signi icance o di e ences among popula ions and con i m ha
C es ed Bu e has high FST (Figu e 5). Finally, FST was always highes
FIGURE 4 The obse ed coun s o candida e SNPs (i.e., SNPs
wi h a q- alue <0.05) ac oss ch omosomes o each en i onmen al
a iable (see Sec ion 2). The o al numbe o SNPs, and he
p opo ion o all SNPs, on each ch omosome a e gi en below each
se o ba s. The expec ed coun s on each ch omosome, ob ained
om he p opo ions o all SNPs ac oss ch omosomes, a e shown
as poin s aligned wi h each ba
Re e ence N genes N o e laps
Mean (95% CI)
esampled N o e laps
Empi ical
p- alue
Kanka e e al., (2010) 14 82.77 (0, 6) .0004
Pa ke e al., (2015) 42 13 6.28 (2, 11) .001
Kanka e e al., (2016) 3,929 946 773.14 (732, 814) <.0001
Pa ke e al., (2016) 130 27 18.72 (11, 26) .01
Pa ke e al., (2018) 2,114 629 414.7 (382, 448) <.0001
Kau anen e al., (2019) 1,751 402 344.64 (314, 375) <.0001
Shown a e he ci a ions o he o iginal s udy, he numbe o genes iden i ied om he o iginal
s udy, he numbe o o e lapping genes wi h he candida e gene se om his s udy, and he mean
and 95% CI om 10,000 esampled se s o candida e genes.
TABLE 2 The o e lap o he union o
genes wi hin 10 kb o SNPs associa ed wi h
PC1, PC2, CTmin, and CCRT and p e ious
candida e gene se s (see Table S7)