gene-09-00728 Janua y 7, 2019 Time: 16:56 # 1
ORIGINAL RESEARCH
published: 11 Janua y 2019
doi: 10.3389/ gene.2018.00728
Edi ed by:
Ino Cu ik,
Uni e si y o Zag eb, C oa ia
Re iewed by:
Luiz Lehmann Cou inho,
Uni e si y o São Paulo, B azil
E eline M. Ibeagha-Awemu,
Ag icul u e and Ag i-Food
Canada (AAFC), Canada
*Co espondence:
Juha Kan anen
[email p o ec ed]
Special y sec ion:
This a icle was submi ed o
Li es ock Genomics,
a sec ion o he jou nal
F on ie s in Gene ics
Recei ed: 15 Augus 2018
Accep ed: 22 Decembe 2018
Published: 11 Janua y 2019
Ci a ion:
Weldenegodguad M, Popo R,
Pokha el K, Ammoso I, Ming Y,
I ano a Z and Kan anen J (2019)
Whole-Genome Sequencing o Th ee
Na i e Ca le B eeds O igina ing F om
he No he nmos Ca le Fa ming
Regions. F on . Gene . 9:728.
doi: 10.3389/ gene.2018.00728
Whole-Genome Sequencing o Th ee
Na i e Ca le B eeds O igina ing
F om he No he nmos Ca le
Fa ming Regions
Melak Weldenegodguad1,2, Ruslan Popo 3, Kisun Pokha el1, Innoken yi Ammoso 4,
Yao Ming5, Zoya I ano a3and Juha Kan anen1*
1Depa men o P oduc ion Sys ems, Na u al Resou ces Ins i u e Finland (Luke), Helsinki, Finland, 2Depa men
o En i onmen al and Biological Sciences, Uni e si y o Eas e n Finland, Kuopio, Finland, 3Yaku ian Resea ch Ins i u e
o Ag icul u e (FGBNU Yaku skij NIISH), Yaku sk, Russia, 4Boa d o Ag icul u al O ice o E eno-By an aj Region,
Ba agay-Aly a, Russia, 5BGI-Genomics, BGI-Shenzhen, Shenzhen, China
No he n Fennoscandia and he Sakha Republic in he Russian Fede a ion ep esen
he no he nmos egions on Ea h whe e ca le a ming has been adi ionally p ac iced.
In his s udy, we pe o med whole-genome sequencing o gene ically cha ac e ize
h ee a e na i e b eeds Eas e n Finnca le, Wes e n Finnca le and Yaku ian ca le
adap ed o hese no he n Eu asian egions. We examined he demog aphic his o y,
gene ic di e si y and un olded loci unde na u al o a i icial selec ion. On a e age,
we achie ed 13.01- old genome co e age a e mapping he sequencing eads on
he bo ine e e ence genome (UMD 3.1) and de ec ed a o al o 17.45 million single
nucleo ide polymo phisms (SNPs) and 1.95 million inse ions-dele ions (indels). We
obse ed ha he ances al species (Bos p imigenius) o Eu asian au ine ca le
expe ienced wo no able p ehis o ical declines in e ec i e popula ion size associa ed
wi h d ama ic clima e changes. The mode n Yaku ian ca le exhibi ed a highe le el o
wi hin-popula ion a ia ion in e ms o numbe o SNPs and nucleo ide di e si y han he
con empo a y Eu opean au ine b eeds. This esul is in con as o he esul s o ma ke -
based ca le b eed di e si y s udies, indica ing asso men bias in p e ious analyses. Ou
esul s sugges ha he e ec i e popula ion size o he ances al Asia ic au ine ca le
may ha e been highe han ha o he Eu opean ca le. Al e na i ely, ou indings could
indica e he hyb id o igins o he Yaku ian ca le ances ies and possibly he lack o
in ensi e a i icial selec ion. We iden i ied a numbe o genomic egions unde selec ion
ha may ha e con ibu ed o he adap a ion o he no he n and suba c ic en i onmen s,
including genes in ol ed in disease esis ance, senso y pe cep ion, cold adap a ion and
g ow h. By cha ac e izing he na i e b eeds, we we e able o ob ain new in o ma ion on
ca le genomes and on he alue o he adap ed b eeds o he conse a ion o ca le
gene ic esou ces.
Keywo ds: adap a ion, demog aphic his o y, Finnca le, indels, selec i e sweeps, SNPs, WGS, Yaku ian ca le
Abb e ia ions: ∂a∂i, di usion app oxima ion o demog aphic in e ence; CI, con idence in e al; CLR, composi e
likelihood a io; Gb, gigabases; GO, gene on ology; NA, ances al popula ion size; Ne, e ec i e popula ion size; nsSNPs,
non-synonymous SNPs; PCA, p incipal componen analysis; PSMC, pai wise sequen ially Ma ko ian coalescen ; SFS, si e
equency spec um.
F on ie s in Gene ics | www. on ie sin.o g 1Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 2
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
INTRODUCTION
Du ing hei 8,000–10,000 yea s o domes ica ion, au ine ca le
(Bos au us) ha e adap ed o a wide a ie y o biogeog aphic
zones and sociocul u al en i onmen s as a esul o na u al and
human-de i ed selec ion (Felius, 1995). Fennoscandia along wi h
no hwes e n Russia and he egion o Sakha (Yaku ia) in eas e n
Sibe ia, a e he no he nmos e i o ies whe e ca le a ming
has had a ela i ely long adi ion as he li elihood o local
people (Kopo e a and Pa anen, 2009;Bläue and Kan anen,
2013;C amp e al., 2014;Ego o e al., 2015). In p ehis o ic
and his o ic imes, animal husband y aced se e al challenges in
hese no he n clima ic condi ions, such as sho summe s and
limi ed ege a ion esou ces o eeding du ing he long win e s,
and his p ac ice equi ed well-adap ed animals ha we e sui ed
o he a ailable en i onmen al esou ces and socioeconomic and
cul u al condi ions (Kan anen e al., 2009a;Bläue and Kan anen,
2013;Ego o e al., 2015).
Ca le b eeds such as Eas e n Finnca le, Icelandic ca le,
Swedish Moun ain ca le, Yaku ian ca le and o he no he n
na i e ca le b eeds a e assumed o ha e hei o igins in he
nea -eas e n domes ica ed au ine ca le ha once sp ead o
hese no he n egions (Kan anen e al., 2000, 2009a;Li e al.,
2007). He d books, pedig ee egis e s and b eeding associa ions
we e es ablished in he la e 19 h and ea ly 20 h cen u ies. Ea ly
na i e b eeds had a pi o al socioeconomic ole in dai y and
bee p oduc ion in he no he n Eu asian egions bu ha e
been almos exclusi ely eplaced by comme cial in e na ional
ca le popula ions b ed o high-inpu , high-ou pu a ming
sys ems. Excep ions o his end a e Yaku ian ca le in Sibe ia
and Icelandic ca le, which con inue o ha e high egional
impo ance in ood p oduc ion (Kan anen e al., 2000, 2009a).
The conse a ion o he gene ic esou ces o na i e, ypically low-
p o i b eeds is o en mo i a ed by he ac ha hese b eeds may
possess aluable gene ic a ia ions o u u e animal b eeding
and o add ess he challenges ha animal p oduc ion will ace
du ing adap a ion o u u e condi ions, b ough abou by ac o s
such as clima e change (Odegå d e al., 2009;Boe che e al.,
2010;Kan anen e al., 2015). In addi ion, b eeds such as Yaku ian
ca le exhibi adap a ion in demanding en i onmen s and may
be ex emely use ul o enabling animal p oduc ion in ma ginal
egions (Kan anen e al., 2015).
P e ious s udies on he cha ac e iza ion o ca le gene ic
esou ces in no he n Eu asian b eeds ha e used a ious me hods
o s udy wi hin-b eed gene ic di e si y, popula ion s uc u e,
demog aphic ac o s and in e b eed ela ionships, e.g., au osomal
and Y-ch omosomal mic osa elli es, mi ochond ial D-loop and
whole-genome SNP-ma ke scans (Li e al., 2007;Kan anen e al.,
2009b;Iso-Tou u e al., 2016). These s udies ha e indica ed,
o example, he gene ic dis inc i eness o he na i e no he n
Eu opean ca le b eeds (e.g., he Finnish na i e b eeds and
Yaku ian ca le) om mode n comme cial dai y b eeds (such
as he Finnish Ay shi e and Hols ein b eeds). In addi ion,
a whole-genome SNP geno yping analysis de ec ed genomic
egions a ge ed by selec ion, which, o example, con ain
immune and en i onmen al adap a ion ela ed genes (Iso-Tou u
e al., 2016;Yu chenko e al., 2018). Whole-genome sequencing
(WGS)-based app oaches p o ide addi ional possibili ies o
in es iga ion o he gene ic di e si y o li es ock b eeds adap ed
o a ious biogeog aphic egions and p oduc ion en i onmen s.
Mo eo e , ecen ad ancemen s in bioin o ma ics and s a is ical
ools ha e enhanced ou unde s anding o he demog aphic
e olu ion o domes ic animal species, he possible ole o genomic
s uc u al a ia ions in he adap a ion o li es ock b eeds in
he cou se o domes ica ion and selec ion and he biological
unc ions o hese genomic a ia ions (Gu enkuns e al., 2009;
Li and Du bin, 2011;Alachio is e al., 2012;Pa lidis e al., 2013;
Wang G.-D. e al., 2014; Wang M. e al., 2014;Lib ado e al.,
2015).
To expand ou knowledge o genomic a ia ions in
no he n Eu asian au ine ca le, we pe o med whole-genome
sequencing o i e animals om each o h ee no he n
na i e b eeds, namely, Eas e n Finnca le, Wes e n Finnca le,
and Yaku ian ca le (Figu e 1). We examined he gene ic
di e si y and popula ion s uc u es o he b eeds and iden i ied
ch omosomal egions and genes unde selec ion p essu e.
We also s udied he demog aphic his o y o he no he n
Eu asian au ine ca le by using he whole-genome sequence
da a.
MATERIALS AND METHODS
E hics S a emen
Blood samples o animals o DNA ex ac ion we e collec ed
by using a p o ocol app o ed by he Animal Expe imen Boa d
o MTT Ag i ood Resea ch Finland (cu en ly he Na u al
Resou ces Ins i u e Finland, Luke) and he Boa d o Ag icul u al
O ice o E eno-By an aj Region, Sakky y , Sakha, Russia.
DNA Sample P epa a ion and
Sequencing
DNA ex ac ed om blood samples was a ailable o he wo
Finnish ca le b eeds (Eas e n Finnca le and Wes e n Finnca le)
and one Sibe ian b eed (Yaku ian ca le) om a p e ious s udy
(Li e al., 2007). Fi e un ela ed indi iduals om each b eed (14
emales and one Yaku ian ca le bull) we e examined. Genomic
DNA was ex ac ed using a s anda d phenol/chlo o o m-based
p o ocol (Malke, 1990). Fo sequencing lib a y p epa a ion
ollowing he manu ac u e ’s speci ica ions, he genomic DNA o
each indi idual was agmen ed andomly. A e elec opho esis,
DNA agmen s o desi ed leng h we e gel pu i ied. One ype
o lib a y was cons uc ed o each sample (500 bp inse
size); 15 pai ed-end DNA lib a ies we e cons uc ed o he 15
samples. Adap e liga ion and DNA clus e p epa a ion we e
pe o med, and he DNA was subjec ed o Illumina HiSeq 2000
sequencing using he 2 ×100 bp mode a Beijing Genomics
Ins i u e (BGI, Shenzhen, China). Finally, pai ed-end sequence
da a we e gene a ed. To ensu e quali y, he aw da a was
modi ied by he ollowing wo s eps using SOAPnuke (Chen
e al., 2018a,b): i s , he con amina ing adap e sequences om
he eads we e dele ed, and hen, he eads ha con ained
mo e han 50% low-quali y bases (quali y alue ≤5) we e
emo ed.
F on ie s in Gene ics | www. on ie sin.o g 2Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 3
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
FIGURE 1 | Th ee No h Eu asian na i e ca le b eeds a e included in his s udy. (A) Eas e n Finnca le a e ypically ed-sided and polled. Ca le b eeding in Finland
was s a ed wi h his b eed, and he b eed’s he d book was es ablished in 1898. The b eed was h ea ened wi h ex inc ion in he 1970s and 1980s. The cu en
census size is 1,600 cows, and he annual milk yield on a e age 4,000 Kg. (B) Wes e n Finnca le a e solid ligh o da k b own and polled. The b eed is one o he
mos p oduc i e na i e ca le b eeds: he a e age annual milk yield is abou 7,000 Kg. (C) The Yaku ian ca le a e cha ac e ized by being pu eb ed abo iginal na i e
ca le om Sakha. Adul Yaku ian cows weigh ypically 350–400 Kg and hei heigh a he wi he s is 110–112 cm on a e age. The animals a e well adap ed o
Sibe ian ha sh condi ions whe e he empe a u e alls below –50◦C in long win e s. The a e age annual milk yield is 1,000 Kg. Rep oduced wi h pe mission om Ulla
Rams adius, Na u al Resou ces Ins i u e Finland (A), Ki s i Hassinen (B) and Anu Os a (C).
Sho Read Alignmen and Mapping
Fo sho ead alignmen , he bo ine e e ence genome
(UMD 3.1), including egions ha we e no assembled in o
ch omosomes (Zimin e al., 2009), we e downloaded om he
Ensembl da abase elease 71 (Flicek e al., 2013) and indexed
using SAM ools 0.1.19 (Li e al., 2009). Pai ed-end 100-bp sho
eads om each indi idual sample we e mapped agains he
bo ine e e ence genome assembly UMD 3.1 using BWA 0.7.5a
wi h he de aul pa ame e s. A e mapping, o downs eam
SNP and inse ion-dele ion (indel) de ec ion, he SAM iles ha
we e gene a ed om BWA we e con e ed o he co esponding
bina y equi alen BAM iles and so ed simul aneously using
So Sam.ja in Pica d ools 1.1021. We used Pica d ools o
emo e PCR duplica es om he aligned eads and hen used he
uniquely mapped eads o a ian calling.
SNP and Indel De ec ion
We used he Genome Analysis Toolki (GATK) 2.6-4
acco ding o he GATK bes p ac ices pipeline (McKenna
e al., 2010;DeP is o e al., 2011;Van de Auwe a A. G.
e al., 2013) o downs eam SNP and indel calling. We used
Realigne Ta ge C ea e o iden i y poo ly mapped egions
(nea by indels) om he alignmen s and ealigned hese
egions using IndelRealigne . Nex , he Uni iedGeno ype
was used o call SNPs and indels wi h a Ph ed scale quali y
g ea e han 30. A e SNP calling, we used Va ian Fil a ion
o disca d sequencing and alignmen a i ac s om he
SNPs wi h he pa ame e s “MQ0 ≥4 && ((MQ0/(1.0 ∗
DP)) >0.1)”, “SB ≥ −1.0, QUAL <10,” and “QUAL <30.0
| | QD <5.0 | | HRun >5 | | SB >−0.10” and om he
indels wi h he pa ame e s “QD <2.0,” “FS >200.0,” and
“ReadPosRankSum <−20.0.” All he a ian s ha passed he
abo e il e ing c i e ia we e used in he downs eam analysis and
compa ed o he ca le dbSNP150 (Van de Auwe a A. G. e al.,
2013) o iden i y no el a ian s.
1h p://pica d.sou ce o ge.ne /
SNP and Indel Anno a ion and Gene
On ology Analysis
ANNOVAR (Wang e al., 2010) was used o anno a e he
unc ions o he a ian s (exonic, in onic, 50and 30UTRs,
splicing, in e genic) using Ensembl elease 71. SNPs ha we e
iden i ied in he exonic egions we e classi ied as synonymous o
nsSNPs. We pe o med GO analysis o genes con aining nsSNPs
and indels using he GO Analysis Toolki and Da abase o
Ag icul u al Communi y (Ag iGO) (Du e al., 2010). Following
he app oaches by Kawaha a-Miki e al. (2011),Liao e al. (2013),
and Li e al. (2014), we selec ed genes con aining >5 nsSNPs o
each b eed. The signi ican ly en iched GO e ms we e assessed
by Fishe ’s exac es wi h he Bon e oni co ec ion using de aul
pa ame e s (P- alue, 0.05; a leas 5 mapping en ies). Ou o ou
indel classes ( ameshi , non- ameshi , s opgain, and s oploss),
we anno a ed ameshi indels in exonic egions using de aul
pa ame e s in ANNOVAR. F ameshi indels may change amino
acid sequences and he eby a ec p o ein unc ion.
Iden i ica ion and Anno a ion o
Selec i e Sweeps
We in es iga ed o all iden i ied SNPs he signa u es o selec ion
using SFS-based αs a is ics in SweeD (Pa lidis e al., 2013)
wi h de aul pa ame e s, excep se ing he g id as he only
pa ame e . SweeD de ec s he signa u e o selec ion based on he
CLR es using SFS-based s a is ics. SweeD was un sepa a ely
o each ch omosome by se ing he g id pa ame e a 5-
kb equidis an posi ions ac oss he ch omosome (size o he
ch omosome/5 kb). We used BEAGLE p og am e .4 (B owning
and B owning, 2007) o impu e missing alleles and in e he
haplo ype phase o all indi idual Wes e n Finnca le, Yaku ian
ca le, and Eas e n Finnca le simul aneously (among he Eas e n
Finnca le, we excluded one inb ed animal; see sec ion “Resul s”).
The BEAGLE p og am in e s he haplo ype in o ma ion o each
ch omosome, which is equi ed o αs a is ics. Following he
app oaches desc ibed in p e ious s udies (Wang M. e al., 2014;
F on ie s in Gene ics | www. on ie sin.o g 3Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 4
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
McManus e al., 2015), we selec ed he ou lie s alling wi hin he
op 0.5% o he CLR dis ibu ion. The cu o alue o αs a is ics
was aken as he 99.5 pe cen ile o he empi ical dis ibu ion
o he 5-kb equidis an posi ions ac oss he genome o each
ch omosome. Anno a ion o he candida e si es ha exhibi ed
a signal o selec ion was pe o med using Ensembl BioMa
(Kinsella e al., 2011) by conside ing a 150-kb sliding window
on he ou lie si es. Candida e genes exhibi ing signa u es o
selec ion we e subjec ed o GO analysis wi h same pa ame e s
applied in he a ian anno a ion using Ag iGO.
Popula ion Gene ics Analysis
The a e age pai wise nucleo ide di e si y wi hin a popula ion
(π) and he p opo ion o polymo phic si es (Wa e son’s θ)
we e compu ed using he Bio::PopGen::S a is ics package in
BioPe l ( 1.6.924) (S ajich e al., 2002). PCA was conduc ed using
sma pca in EIGENSOFT3.0 so wa e (Pa e son e al., 2006) on
biallelic au osomal SNPs ha we e geno yped in all indi iduals.
Signi ican eigen ec o s we e de e mined using T acy–Widom
s a is ics wi h he ws a s p og am implemen ed in he same
EIGENSOFT package.
Demog aphic His o y In e ence
We used he PSMC model (Li and Du bin, 2011) o cons uc
he demog aphic his o y o he h ee b eeds. Fo he analysis,
one indi idual pe b eed wi h he highes sequence dep h was
selec ed o explo e changes in local densi y o he e ozygous
si es ac oss he ca le genome. The ollowing de aul PSMC
pa ame e s we e se : −N25, − 15, − 5 and −p ‘4+25∗2+4+6’.
To scale he PSMC ou pu o eal ime, we assumed a neu al
mu a ion a e o 1.1 ×10−8pe gene a ion and an a e age
gene a ion ime o 5 yea s (Kuma and Sub amanian, 2002;
Mu ay e al., 2010;MacLeod e al., 2013). As he powe o
he PSMC app oach o econs uc ecen demog aphic his o y
is no eliable (Li and Du bin, 2011;MacLeod e al., 2013;
Zhao e al., 2013), we econs uc ed a mo e ecen demog aphic
his o y o he Finnish and Yaku ian popula ions using he ∂a∂i
p og am (dadi-1.6.3) (Gu enkuns e al., 2009). We used he
in e genic si es om he iden i ied SNPs in he 15 indi iduals o
compu e he olded SFS. We me ged he esul s o he Eas e n
and Wes e n Finnca le b eeds, as hese b eeds exhibi ed simila
gene ic di e si y measu es (Supplemen a y Figu e S4). Since we
had 10 Finnca le and 5 Yaku ian samples, we downscaled he
Finnca le sample size o be equal o ha o he Yaku ian ca le.
We an he ∂a∂i algo i hm mul iple imes o ensu e con e gence
and selec ed he op imal pa ame e s wi h he highes likelihood
as he inal esul . As ∂a∂i equi es NA, we calcula ed NA using
he o mula NA = θ/4µL, whe e θwas he obse ed numbe o
seg ega ing si es di ided by he sum o he expec ed SFS using
he bes - i pa ame e s o ou model, L was he e ec i e sequence
leng h, and µwas he mu a ion a e pe gene a ion pe si e. We
used a mu a ion a e o 1.1 ×10−8mu a ions pe gene a ion
assuming ha one gene a ion was equal o 5 yea s (Kuma
and Sub amanian, 2002), and he e ec i e sequence leng h
(in e genic egions) was 10,836,904. We calcula ed popula ion
size and di e gence ime be ween he Finnish and Yaku ian
popula ions based on NA. Finally, using he pa ame e s desc ibed
p e iously, we gene a ed he demog aphic model using ∂a∂i
as shown in Supplemen a y Figu e S5. The op imal model
iden i ied he change om he NA o he e ec i e popula ion size
(nua) om he ime Ta o he ime Td. Ta was he ime pe iod
when he change in NA s a ed and Td was he ime when he
di e gence be ween he Finnish and Yaku ian ca le occu ed.
nu1F and nu2Y we e he e ec i e popula ion sizes du ing he
spli . To calcula e he s a is ical con idence in he es ima ed
pa ame e alues, we es ima ed he pa ame e unce ain ies using
he Hessian me hod (a.k.a. he Fishe in o ma ion ma ix).
RESULTS
Sequence Da a
A o al o 521 Gb o pai ed-end DNA sequence da a was ob ained
a e emo ing adap e sequences and low-quali y eads (Table 1
and Supplemen a y Table S1). On a e age, each sample had
347.4 million (M) eads, 98.45% o which we e success ully
mapped o he bo ine e e ence genome UMD3.1 (Table 1 and
Supplemen a y Table S1), ep esen ing 13.01- old co e age.
Iden i ica ion and Anno a ion o Va ian s
A o al o 17.45 M SNPs we e de ec ed in he mapped eads ac oss
all 15 samples, wi h Yaku ian ca le exhibi ing he highes numbe
o SNPs (Table 2,Figu e 2A and Supplemen a y Table S2).
The a e age numbe o SNPs de ec ed pe indi idual wi hin he
b eeds was 5.73, 6.03, and 7.12 M in Eas e n Finnca le, Wes e n
Finnca le and Yaku ian ca le, espec i ely (Supplemen a y
Table S2). A o al o 6.3 M (36.1%) SNPs we e sha ed by he h ee
b eeds, and as expec ed, he Finnish b eeds sha ed he highes
numbe (n= 8.06 M, 46.2%) o SNPs (Figu e 2A). Mo eo e , we
TABLE 1 | Summa y o sequencing and sho ead alignmen esul s.
Eas e n
Finnca le
Wes e n
Finnca le
Yaku ian
ca le
O e all
sample
Numbe o
indi iduals
5 5 5 15
Pai ed-end leng h
(bp)
100 100 100 100
A e age eads
pe indi idual
352.73 M 347.14 M 342.33 M 347.4 M
A e age
sequence dep h
pe indi iduala
13.21X 13.00X 12.82X 13.01X
A e age map
eads pe
indi idual
348.42 M 340.12 M 337.50 M 342.02 M
A e age unique
map eads pe
indi idual
316.89 M 312.58 M 309.19 M 312.88 M
A e age ead
mapping a e
98.78% 97.97% 98.59% 98.45%
A e age
co e age a e
98.42% 98.22% 98.46% 98.37%
aA e age sequence dep h pe indi idual was compu ed by di iding he clean eads
by he e e ence genome size.
F on ie s in Gene ics | www. on ie sin.o g 4Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 5
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
TABLE 2 | Func ional anno a ion o he de ec ed SNPs and indels.
Eas e n
Finnca le
Wes e n
Finnca le
Yaku ian
ca le
SNP To al numbe o
SNPs
11,017,215 10,543,290 12,242,166
In e genic 7,998,914 7,662,604 8,845,911
In onic 2,764,951 2,643,821 3,114,622
Exonica
Non-
synonymous
30,982 28,733 32,782
S op gain 294 284 310
S op loss 23 18 19
Synonymous 41,111 38,137 46,950
Ups eam 74,552 69,369 82,170
Downs eam 74,198 70,815 83,549
Ups eam;
downs eamb
1,636 1,534 2,036
UTRc25,545 23,342 28,273
Splicing 417 388 459
ncRNA 4,593 4,256 5,086
Indel To al numbe o
indels
1,275,128 1,188,892 1,374,577
In e genic 942,143 878,007 1,012,733
In onic 332,504 310,089 363,783
Exonica
Non- ameshi 397 327 427
S op gain 24 22 27
S op loss 1 1 0
F ameshi 1,045 972 1,148
Ups eam 9,269 8,286 9,861
Downs eam 10,611 9,609 11,377
Ups eam;
downs eamb
250 218 282
UTRc3,380 3,101 3,767
Splicing 248 233 268
ncRNA 406 371 437
aExonic = “exonic” and “exonic; splicing” as anno a ed by ANNOVAR.
bUps eam; downs eam = a ian loca ed in downs eam and ups eam egions.
cUTR = “UTR3” and “UTR5” as anno a ed by ANNOVAR.
ound ha 1.85 M SNPs (16.83%) in Eas e n Finnca le, 1.60 M
(15.15%) in Wes e n Finnca le and 3.96 M (32.33%) in Yaku ian
ca le we e p i a e SNPs in ou b eed se (Figu e 2A). The
ansi ion- o- ans e sion (TS/TV) a ios we e 2.20 and 2.23 in
he Finnca le and Yaku ian ca le, espec i ely (Supplemen a y
Table S2). The obse ed Ts/T a ios we e consis en wi h hose
obse ed in p e ious s udies in mammalian sys ems (Lachance
e al., 2012;Choi e al., 2013, 2014), indica ing he quali y o ou
SNP da a.
O he SNPs iden i ied in ou analysis, 1.07 M (6.13%)
SNPs we e ound o be no el when compa ed o NCBI dbSNP
bo ine build 150. A he b eed le el, 2.8, 2.6, and 4.5% o
he o al SNPs in he Eas e n Finnca le, Wes e n Finnca le
and Yaku ian ca le, espec i ely, we e no el. Fu he mo e, ou
o he no el SNPs iden i ied o each b eed, 258,409 (83.3%),
219,234 (81.5%), and 528,763 (95%) we e b eed-speci ic SNPs
in Eas e n Finnca le, Wes e n Finnca le, and Yaku ian ca le
(Figu e 2B), espec i ely. A summa y o he homozygous and
he e ozygous SNPs is gi en in Supplemen a y Tables S2, S3. One
Eas e n Finnca le cow (sample_3 in Supplemen a y Table S3)
exhibi ed excep ionally low di e si y, wi h only 1.66 M (32.58%)
he e ozygous and 3.44 M (67.42%) homozygous SNPs. This
animal o igina ed om an isola ed, inb ed he d and ep esen ed
one elic Eas e n Finnca le line (he d) ha passed h ough he
b eed’s demog aphic bo leneck (Kan anen e al., 2000). A e
excluding his sample, he a e age numbe o SNPs de ec ed
pe Eas e n Finnca le indi idual was 5.88 M, and he Eas e n
Finnca le animals exhibi ed 2.63 M (44.83%) homozygous
and 3.24 M (55.17%) he e ozygous SNPs, wi h a a io o
1:1.23 (homozygous:he e ozygous). Appa en ly, he numbe o
homozygous SNPs in he Eas e n Finnca le was highe han ha
in he o he wo b eeds.
In o al, we de ec ed 2.12 M indels, 79.72% o which we e
ound in he dbSNP build 150, wi h 20.28% being no el
(Figu e 2C and Supplemen a y Table S2). A he b eed le el,
12.9, 11.6, and 16% o he o al indels in he Eas e n Finnca le,
Wes e n Finnca le, and Yaku ian ca le, espec i ely, we e no el.
In ou da a, on a e age, 0.65% o he SNPs we e de ec ed in
exonic egions, 25.1% in in onic egions, 72.6% in in e genic
egions, and 1.65% in UTRs and in egions ups eam and
downs eam o genes (Table 2 and Supplemen a y Table S4).
In gene al, all he h ee b eeds exhibi ed simila dis ibu ions
o SNPs in a ious unc ional ca ego ies. A o al o 76,810,
71,256, and 84,927 exonic SNPs we e iden i ied in he Eas e n
Finnca le, Wes e n Finnca le and Yaku ian ca le, espec i ely.
O he exonic SNPs in he Eas e n Finnca le, Wes e n Finnca le,
and Yaku ian ca le, 31,299, 29,035 and 33,111, espec i ely, we e
nsSNPs (Table 2) and we e ound in 10,309, 9,864, and 10,429
genes, espec i ely.
The unc ional ca ego ies o he indel mu a ions a e p esen ed
in Table 2. In o al, 1,045, 927, and 1,148 o he indels we e
ameshi indels ha we e associa ed wi h 808, 770, and 895
genes in Eas e n Finnca le, Wes e n Finnca le, and Yaku ian
ca le, espec i ely (Supplemen a y Da as S1–S3).
We compu ed he geno ype conco dance be ween he SNPs
de ec ed by ou SNP calling pipeline and he p e ious SNP
geno ype s udy (Iso-Tou u e al., 2016); Illumina Bo ineSNP50
BeadChip .1.0 (Ma ukumalli e al., 2009). On a e age, we ound
85.34% o SNPs de ec ed by sequencing we e conco dan wi h he
SNP50 BeadChips geno ypes sugges ing ha he SNPs de ec ed
in ou s udy wi h he cu en co e age (on a e age 13.01- old
co e age/sample) yielded su icien geno ypic accu acy.
GO Analysis o he SNPs and Indels
The GO en ichmen analysis o 1,331, 1,170, and 1,442 genes
con aining >5 nsSNPs (Supplemen a y Da as S4–S6), iden i ied
111, 113, and 95 signi ican ly en iched GO e ms in Eas e n
Finnca le, Wes e n Finnca le and Yaku ian ca le, espec i ely
(Supplemen a y Da as S7–S9). A o al o 38, 43 and 38 GO
e ms we e associa ed wi h biological p ocesses in Eas e n
Finnca le, Wes e n Finnca le, and Yaku ian ca le, espec i ely
(Supplemen a y Da as S7–S9).
F on ie s in Gene ics | www. on ie sin.o g 5Janua y 2019 | Volume 9 | A icle 728

gene-09-00728 Janua y 7, 2019 Time: 16:56 # 6
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
FIGURE 2 | Venn diag am showing o e lapping and unique SNPs/indels be ween he h ee b eeds (Eas e n Finnca le, Wes e n Finnca le, and Yaku ian). The
numbe s in pa en heses ou side he ci cles a e he o al numbe o de ec ed SNPs om each b eed. The numbe s in he ci cle componen s show speci ic SNPs o
each b eed o o e lapping SNPs/indels be ween any wo b eeds o among h ee b eeds. (A) The iden i ied sha ed and speci ic SNPs o each b eed, (B) he
iden i ied sha ed and speci ic no el SNPs o each b eed, and (C) he iden i ied sha ed and speci ic indels o each b eed.
A de ailed compa ison o he biological p ocesses associa ed
wi h genes wi h >5 nsSNPs wi h he bo ine Ensembl gene se
(n= 25,160) is shown in Supplemen a y Figu e S1. The GO
en ichmen analysis e ealed ha a majo i y o he signi ican ly
en iched GO e ms we e sha ed by he h ee ca le b eeds.
“Response o s imulus, GO:00050896” was associa ed wi h
app oxima ely 50% o he genes in Eas e n Finnca le (n= 611),
Wes e n Finnca le (n= 544), and Yaku ian ca le (n= 629) (see
Supplemen a y Figu e S1). In addi ion, his analysis showed ha
in each b eed, a la ge numbe o genes we e associa ed wi h
immune unc ions, such as “Immune esponse, GO:0006955,”
“De ense esponse, GO:0006952,” “An igen p ocessing and
p esen a ion, GO:0019882,” and “Immune sys em p ocess,
GO:0002376.” Among he h ee b eeds, he Yaku ian ca le had
mo e en iched genes associa ed wi h immune unc ions han
he wo Finnca le b eeds. On he o he hand, in he Finnca le
b eeds, a la ge numbe o genes we e associa ed wi h senso y
pe cep ion unc ions, such as “Senso y pe cep ion, GO:0007600,”
“Senso y pe cep ion o smell, GO:0007608,” and “De ec ion o
chemical s imulus in ol ed in senso y pe cep ion, GO:0050907.”
In Yaku ian ca le, none o he GO e ms associa ed wi h senso y
pe cep ion we e en iched. Howe e , 55 genes associa ed wi h
“De elopmen al g ow h, GO: 0048589” we e en iched in only
Yaku ian ca le.
We u he iden i ied he op genes, namely, TTN, PKHD1,
GPR98, and ASPM, ha had a leas 40 nsSNPs in all he b eeds.
These genes ha e la ge sizes; TTN is 274 kb in size, PKHD1 is
455 kb, GPR98 is 188 kb and ASPM is 64 kb. Among he genes
wi h nsSNPs, TTN con ained he highes numbe o nsSNPs: 68,
63, and 87 nsSNPs in Eas e n Finnca le, Wes e n Finnca le,
and Yaku ian ca le, espec i ely. The TTN gene is p esen on
ch omosome 2 and is associa ed wi h mea quali y (Sasaki e al.,
2006;Wa anabe e al., 2011).
A o al o 709, 675, and 772 genes associa ed wi h ameshi
indels in hese b eeds we e linked o a leas one GO e m
(Supplemen a y Figu e S2 and Supplemen a y Da as S10–S12).
The esul s indica ed ha a majo i y o he signi ican ly en iched
GO e ms we e sha ed by he b eeds. The GO e ms “De ense
F on ie s in Gene ics | www. on ie sin.o g 6Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 7
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
esponse, GO:0006952” and “Female p egnancy, GO:0007565”
we e en iched exclusi ely in Yaku ian ca le. In o al, 96 genes
we e en iched in “De ense esponse, GO:0006952.”
Selec ion Signa u es
We iden i ied 2,528 si es exhibi ing signa u es o selec ion in
each b eed, o which 58, 61, and 53% mapped o gene egions
in Eas e n Finnca le, Wes e n Finnca le, and Yaku ian ca le,
espec i ely (Supplemen a y Figu e S3). In o ma ion ega ding
he SNPs ound in selec i e sweep egions in each b eed is shown
in Supplemen a y Table S5.
Ch omosome 1 exhibi ed he highes (n= 159) numbe
o selec ion signals and ch omosome 25 he lowes (n= 43).
Conside ing a 150-kb window cen e ed on he candida e
si e, Wes e n Finnca le exhibi ed he highes numbe
(n= 371) o candida e genes wi h selec ion signa u es,
ollowed by Eas e n Finnca le (n= 331), while Yaku ian
ca le exhibi ed he lowes numbe (n= 249) (Supplemen a y
Da as S13–S15). Appa en ly, 36 (Eas e n Finnca le), 35
(Wes e n Finnca le), and 28 (Yaku ian ca le) candida e
genes con ained >5 nsSNPs (Supplemen a y Da as S16–
S18). Se en genes wi h g ea e han 5 nsSNPs in Eas e n
Finnca le (CCSAP, CEP72, GBP5, LOC100297846, GBP2,
LOC613867, and ENSBTAG00000045571), Wes e n Finnca le
(CDH23, PCDHB4, PCDHB6, PCDHB7, SIRPB1, LOC783488,
and ENSBTAG00000012326), and Yaku ian ca le (FER1L6,
GBP5, ENSBTAG00000015464, ENSBTAG00000025621,
GBP2, ENSBTAG00000039016, and LOC101902869) exhibi ed
he s onges signa u es o selec ion. O he genes wi h
he s onges signa u es o selec ion, one gene each om
Eas e n (ENSBTAG00000045571) and Wes e n Finnca le
(ENSBTAG00000012326) and h ee genes om Yaku ian
ca le (ENSBTAG00000015464, ENSBTAG00000025621,
and ENSBTAG00000039016) lacked gene desc ip ions
(Supplemen a y Table S6).
A o al o 28, 67, and 13 GO e ms we e signi ican ly
en iched in Eas e n Finnca le, Wes e n Finnca le, and Yaku ian
ca le, espec i ely (Supplemen a y Da as S19–S21). We ound
only one signi ican ly en iched GO e m (“GMP binding,
GO:0019002”) ha was sha ed by he h ee ca le b eeds. The
GO e ms “Homophilic cell adhesion, GO:0007156,” “Calcium-
dependen cell–cell adhesion, GO:0016339,” and “Mul icellula
o ganism ep oduc ion, GO:0032504” we e sha ed by he
Finnca le b eeds. Mos o he signi ican ly en iched GO e ms
(23, 62, and 12 in Eas e n Finnca le, Wes e n Finnca le, and
Yaku ian ca le, espec i ely) we e ‘b eed-speci ic’ in ou da a.
In addi ion, we examined he signi ican ly en iched GO e ms
ha we e po en ially in ol ed in cold adap a ion by assuming
ha in ex emely cold en i onmen s, ene gy equi emen is high
and a and lipids a e he main sou ces o ene gy (Liu e al.,
2014). The le els o a y acids, lipids and phospholipids ypically
inc ease wi h dec easing empe a u es (Pu aæ e al., 2011).
The signi ican ly en iched GO e ms associa ed wi h Wes e n
Finnca le included “Lipid localiza ion, GO:0010876,” “Lipid
diges ion, GO:0044241,” “Unsa u a ed a y acid biosyn he ic
p ocess, GO:0006636,” and “Unsa u a ed a y acid me abolic
p ocess, GO:0033559.” Howe e , no signi ican ly en iched GO
e ms associa ed wi h a y acid and lipid me abolism and
biosyn hesis we e iden i ied in Eas e n Finnca le and Yaku ian
ca le.
We examined he candida e selec i e sweep genes in each
b eed. A numbe o genes po en ially associa ed wi h cold
adap a ion (Ca dona e al., 2014) we e p esen in Eas e n
Finnca le (DNAJC28, HSP90B1, AGTRAP, TAF7, TRIP13,
NPPA, and NPPB), Wes e n Finnca le (CD14, COBL, JMJD1C,
KCNMA1, PLA2G4, SERPINF2, SRA1, and TAF7), and Yaku ian
ca le (DNAJC9, SOCS3, TRPC7, SLC8A1 GLP1R, PKLR, and
TCF7L2).
Among he selec i e sweep genes, he e we e se e al genes ha
ha e been p e iously shown o be associa ed wi h domes ica ion-
ela ed changes, such as changes in disease esis ance, neu onal
and b ain de elopmen , g ow h, mea quali y, pigmen a ion,
senso y pe cep ion and milk p oduc ion (Gu ié ez-Gil e al.,
2015). Fo example, he ch omosomal egions exhibi ing
selec i e sweeps in Eas e n Finnca le included genes associa ed
wi h disease esis ance (IFNAR1,IFNAR2,IL10RB, and NOD2),
neu onal and b ain de elopmen (OLIG1), g ow h (ACTA1) and
mea quali y (IGFBP5, NRAP, PC, and S1PR1) (Supplemen a y
Da a S13). In Wes e n Finnca le, selec i e sweeps we e de ec ed
in genes associa ed wi h pigmen a ion (ULBP3), senso y
pe cep ion (LOC521946, LOC783558, and LOC783323), mea
quali y (COX5B, KAT2B, and ITGB3) and disease esis ance
(CD96, CD14, GZMB, and IL17A) (Supplemen a y Da a S14).
Simila ly, selec i e sweep-in luenced genes in Yaku ian ca le
we e associa ed wi h disease esis ance (PFKM, ADAM17,
and SIRPA), senso y pe cep ion (OR13C8, LOC100336881,
LOC101902265, LOC512488, LOC617388, LOC783884,
LOC788031, and LOC789957), mea quali y (ALDH1B1,
CAPNS1, COX7A1, PFKM, SLC8A1, SOCS3, and THBS3) and
milk p oduc ion (MUC1) (Supplemen a y Da a S15).
Popula ion Gene ics Analysis
The o e all genome-wide gene ic di e si y, as measu ed by
Wa e son’s θand pai wise nucleo ide di e si y (π), we e
highe in he Yaku ian ca le (0.001588 and 1.728 ×10−3,
espec i ely) han in Eas e n Finnca le (0.001445 and
1.559 ×10−3, espec i ely) and Wes e n Finnca le (0.001398
and 1.512 ×10−3, espec i ely), and hese esul s we e
inconsis en wi h hose o p e ious s udies based on au osomal
mic osa elli e and SNP da a se s, which showed ha Finnca le
we e mo e di e se han he Yaku ian ca le (e.g., Li and Kan anen,
2010).
We also applied PCA o examine he gene ic ela ionships
among he h ee ca le b eeds. In he PCA plo , he Finnca le and
Yaku ian ca le we e g ouped in he i s eigen ec o s, indica ing
clea gene ic di e en ia ion (Supplemen a y Figu e S4). The
inb ed Eas e n Finnca le animal g ouped sepa a ely om he
o he Finnca le animals.
Demog aphic Popula ion Size His o y
The PSMC p o iles o he con empo a y Finnish and Sibe ian
na i e ca le we e used o cons uc he demog aphic p ehis o y
and e olu ion o ances al popula ions o no he n Eu asian
ca le. As shown in Figu e 3, he empo al PSMC p o iles o he
F on ie s in Gene ics | www. on ie sin.o g 7Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 8
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
FIGURE 3 | Demog aphic his o y o he no he nmos ca le b eeds econs uc ed om h ee ca le genomes, one om each b eed, by using PSMC. The X-axis
shows he ime in housand yea s (Ky ), and he Y-axis shows he e ec i e popula ion size.
h ee ca le genomes ollowed a simila pa e n. The ances al
species o no he n Eu asian au ine ca le, he nea -eas e n
au ochs (Bos p imigenius) (Kan anen e al., 2009a), expe ienced
wo popula ion peaks s a ing a ∼1 Mya and ∼40 kya and wo
bo lenecks a ∼250 and ∼12 kya (Figu e 3). A e he i s
popula ion expansion, he popula ion size declined g adually.
The second popula ion expansion o he ances al wild species
began a ound ∼80 kya and s a ed o decline a ound ∼30 kya,
leading o a second bo leneck.
We also used he ∂a∂i p og am o econs uc he ecen
no he n Eu opean ca le demog aphic his o y ( om 418 kya o
he p esen ). The pa ame e s Ta, Td, nua, nu1F, and nu2Y in he
demog aphic model a e shown and explained in Supplemen a y
Figu e S5 and Supplemen a y Table S7. Based on his model, we
es ima ed ha he e e ence NA was 43,116. The op imal model
i o each pa ame e and CI a e shown in Supplemen a y Table
S7 by ixing NA a 43,116 and gene a ion ime a 5 yea s. Ou
bes - i model indica ed ha he ances al popula ion unde wen
a size change o 51,883 (95% CI, 51,658–52,108) a 418 kya
(95% CI, 413.96–409.47 kya) (Supplemen a y Table S7). This
esul is consis en wi h he PSMC p o ile (Figu e 3). In addi ion,
ou model sugges ed ha he di e gence o No h Eu opean
na i e ca le and Eas Sibe ian Tu ano-Mongolian ype o ca le
occu ed 8,822 yea s ago (95% CI, 8,775–8,869 yea s ago).
DISCUSSION
To ou knowledge, his is he i s whole-genome sequence-
based epo on he gene ic di e si y o Eu asian na i e ca le
(B. au us) b eeds ha ha e adap ed o he no he nmos ca le
a ming egions, e en suba c ic egions. The con empo a y
gene ic esou ces o he Eas e n Finnca le, Wes e n Finnca le
and Yaku ian ca le b eeds s udied a e he esul o a complex
p ocess o gene ic and demog aphic e en s ha occu ed du ing
he domes ica ion and selec ion and e en he e olu ion o he
ances al species o no he n Eu asian au ine ca le, namely, he
nea -eas e n au ochs (B. p imigenius).
Demog aphic E olu ion o
Bos p imigenius
As shown in Figu e 3, he au och species (B. p imigenius)
expe ienced wo no able p ehis o ical popula ion expansions,
a e which he popula ion size declined g adually. The i s
ma ked decline in he e ec i e popula ion size (Ne) occu ed
du ing he Middle Pleis ocene pe iod s a ing a e ∼1 Mya,
which may ha e been associa ed wi h educ ion in global
empe a u es and e en wi h nega i e ac ions o humans
on he au och popula ion (Ba nosky e al., 2004;Hughes
e al., 2007). The second ma ked decline in Ne p io o
domes ica ion was ob iously caused by d ama ic clima e changes
du ing he las glacial maximum (Yokoyama e al., 2000).
Al hough he sequencing dep h a ained in his s udy was
no ideal o PSMC analysis ( ypically >20- old co e age),
ou obse a ions ega ding he empo al changes in he Ne
o he au ochs du ing he Pleis ocene pe iod (Mei e al.,
2018) ollowed he pa e n obse ed o ances al popula ions
o se e al o he domes ic mammalian species, such as pig
[Sus sc o a; (G oenen e al., 2012)], ho se [Equus caballus;
(Lib ado e al., 2016)], and sheep [O is a ies; (Yang e al.,
2016)]. The ∂a∂i esul s con i med he pas luc ua ions in he
p ehis o ical Ne o B.p imigenius (Supplemen a y Table S7),
and he compa ison o he cu en SNP-based es ima ed Ne
o he p esen ca le b eeds [∼100; (Iso-Tou u e al., 2016)]
o he Ne o he co esponding ea ly domes ica ed ances al
popula ions showed ha he e was a d ama ic decline in he
Ne du ing domes ica ion and b eed o ma ion. In addi ion,
ou demog aphic analysis (Supplemen a y Figu e S5) p o ided
new knowledge o he p ehis o y o no he n Eu asian na i e
F on ie s in Gene ics | www. on ie sin.o g 8Janua y 2019 | Volume 9 | A icle 728
gene-09-00728 Janua y 7, 2019 Time: 16:56 # 9
Weldenegodguad e al. Sequencing o No h-Eu asian Na i e Ca le
ca le. As sugges ed by a p e ious s udy (Kan anen e al.,
2009b), bo h he Finnish and Yaku ian na i e ca le descended
om he nea -eas e n au ochs domes ica ed 8,000–10,000 yea s
ago. He e, ou esul s ha e shown ha he wo no he n
Eu asian na i e ca le lineages may ha e al eady di e ged in
he ea ly s age o au ine ca le domes ica ion, mo e han
8,000 yea s ago.
High Gene ic Va iabili y in he Yaku ian
Ca le
The o al numbe o sequence a ian s iden i ied on a e age
in Eas e n Finnca le and Wes e n Finnca le animals (e.g.,
5.88 and 6.03 M SNPs, espec i ely, exhibi ing a mino allele
equency >0.05) co esponded well o numbe s ound ypically
in Eu opean au ine animals. In con as , we ound ha he
Yaku ian ca le exhibi ed a highe numbe o SNPs on a e age pe
indi idual (7.12 M SNPs) han he numbe o SNPs de ec ed in
Eu opean and Asia ic humpless ca le o da e (Tsuda e al., 2013;
Choi e al., 2014;Szyda e al., 2015). Acco ding o (Szyda e al.,
2015) and s udies ci ed he ein, a Eu opean au ine animal may
exhibi on a e age 2.06–6.12, 5.89–6.37, 5.85–6.40, and 5.93 M
SNPs, while (Choi e al., 2014) de ec ed 5.81 M SNPs in a Ko ean
Hols ein ca le indi idual, a b eed ha o igina ed om wes e n
Eu ope and No h Ame ica. Typically, i may be possible o
de ec addi ional SNPs by inc easing he sequencing dep h (Szyda
e al., 2015). In addi ion o he a e age numbe o SNPs pe
indi idual, o al numbe o SNPs and numbe o indels, he
Yaku ian ca le exhibi ed he highes numbe o exonic SNPs
and nsSNPs among he h ee no he n na i e b eeds s udied.
Howe e , al hough he Yaku ian ca le had he highes numbe o
nsSNPs and genes wi h >5 nsSNPs, he unc ional anno a ion o
he exonic SNPs by GO analysis indica ed ha he lowes numbe
o signi ican ly en iched GO e ms was ob ained o he Yaku ian
ca le.
Ou es ima es o he popula ion-le el di e si y o he
Eas e n Finnca le, Wes e n Finnca le, and Yaku ian ca le [ he
nucleo ide di e si y (π) alues we e 1.559 ×10−3, 1.512 ×10−3,
and 1.728 ×10−3, espec i ely, and he p opo ions o
polymo phic si es (θ) we e 0.001445, 0.001398, and 0.001588,
espec i ely] exceed hose ypically ound in Eu opean au ine
ca le b eeds (Kim e al., 2017;Chen e al., 2018a,b;Mei e al.,
2018). We obse ed ha Yaku ian ca le such as he Asia ic
au ine ca le b eeds exhibi high le els o genomic di e si y
in e ms o πand θes ima es. The ypical nucleo ide di e si y
alues o he Eu opean au ine ca le a e >1.0 ×10−3, while
hose o he Asia ic au ine b eeds a e close o ∼2.0 ×10−3
han o 1.0 ×10−3(Kim e al., 2017;Chen e al., 2018a,b;Mei
e al., 2018). We obse ed highe wi hin-popula ion di e si y o
he Yaku ian ca le han ha obse ed o se e al o he au ine
ca le b eeds, which di e s om p e ious es ima es based on
au osomal mic osa elli es and whole-genome SNP da a (Li e al.,
2007;Iso-Tou u e al., 2016), whe e lowe le els o a ia ion
we e obse ed in Yaku ian ca le, indica ing ha he gene ic
a ia ion in Yaku ian ca le has been unde es ima ed. The se o
au osomal mic osa elli es ecommended by FAO ( he Food and
Ag icul u al O ganiza ion o he Uni ed Na ions) o biodi e si y
analysis o ca le b eeds and he design o comme cial SNP
BeadChips used in ca le whole-genome geno yping we e de i ed
mainly om he gene ic da a o wes e n b eeds, causing a bias
in he di e si y es ima es o clea ly gene ically dis inc ca le
b eeds, such as Yaku ian ca le (Li e al., 2007;Iso-Tou u e al.,
2016).
The e could ha e been di e ences in he pas e ec i e
popula ion sizes o he Eu opean and Asia ic au ine ca le, and
he p esen ele a ed genomic di e si y o he Asia ic au ine
ca le b eeds may e lec he highe “ancien ” e ec i e sizes o
he ances al popula ions o he Asia ic au ine b eeds (Chen
e al., 2018a,b). Howe e , he p ehis o y o domes ica ed ca le
in Eas Asia appea s o be mo e complex han p e iously
hough (Zhang e al., 2013;Gao e al., 2017;Chen e al.,
2018a,b), and an addi ional specula i e explana ion o he
ele a ed genomic di e si y in he Yaku ian ca le and se e al
o he Asia ic au ine ca le b eeds (o hei ances al popula ions)
could be ancien in og ession wi h he Eas Asian au ochs
(B. p imigenius) ha li ed in he Eas Asian egion du ing he
a i al o nea -eas e n au ine ca le (Chen e al., 2018a,b). The
p e ious m DNA and Y-ch omosomal di e si y s udy indica ed
he nea -eas e n o igins o he ances al popula ion o he
Yaku ian ca le (Kan anen e al., 2009b). The possible hyb id
o igins o he Yaku ian ca le ances ies may ha e inc eased he
gene ic a ia ion in he ances al popula ion o Yaku ian ca le
seen e en in he cu en popula ion and may ha e played a
pi o al ole in he p ocess o adap a ion o he Yaku ian ca le
o he suba c ic en i onmen in he Sakha Republic, eas e n
Sibe ia.
The high numbe o SNPs and high genomic di e si y
ound in he Yaku ian ca le may be due pa ly o he b eed’s
selec ion his o y: he a i icial selec ion by humans has no been
in ensi e (Kan anen e al., 2009b). The Yaku ian ca le b eed is
an abo iginal au ine popula ion, he gene pool o which has
been shaped by na u al and a i icial selec ion. Howe e , he
cen u ies-old “ olk selec ion” me hods and adi ional knowledge
o he selec ion o he mos sui able animals o he challenging
suba c ic en i onmen ollowed he me hods used by local
people a he han he b eeding implemen ed by o ganiza ions
o ins i u ions (Kan anen e al., 2009a). When compa ed wi h
he Wes e n Finnca le and Eas e n Finnca le in he p esen
s udy, he Yaku ian ca le exhibi ed dis inc ly low numbe s o
candida e genes ha exhibi ed selec ion signa u es (n= 371,
n= 331, and n= 249, espec i ely). Among hese h ee b eeds,
Wes e n Finnca le ha e been subjec ed o he mos in ensi e
a i icial selec ion o milk p oduc ion cha ac e is ics, while
he p oduc ion selec ion p og am o Eas e n Finnca le was
s opped in he 1960s, when he census popula ion size o
his na i e b eed declined apidly. Cu en ly, in i o and in
i o conse a ion ac i i ies a e being implemen ed o Eas e n
Finnca le (and o Wes e n Finnca le and Yaku ian ca le). In
addi ion, al hough Yaku ian ca le had he highes numbe o
genes con aining SNPs (also nsSNPs) among he h ee b eeds,
he GO analysis indica ed ha his b eed had he lowes
numbe o signi ican ly en iched GO e ms (Eas e n Finnca le,
111; Wes e n Finnca le, 113; and Yaku ian ca le, 95). This
di e ence be ween he na i e Finnish ca le and Yaku ian ca le
F on ie s in Gene ics | www. on ie sin.o g 9Janua y 2019 | Volume 9 | A icle 728