Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants

RESEARCH ARTICLE Open Access Genome-wide associa ion analysis o milk yield ai s in No dic Red Ca le using impu ed whole genome sequence a ian s T. Iso-Tou u 1* , G. Sahana 2 , B. Guldb and sen 2 , M. S. Lund 2 and J. Vilkki 1 Abs ac Backg ound: The No dic Red Ca le consis ing o h ee di e en popula ions om Finland, Sweden and Denma k a e unde a join b eeding alue es ima ion sys em. The long his o y o eco ding o p oduc ion and heal h ai s o e s a g ea oppo uni y o s udy p oduc ion ai s and iden i y causal a ian s behind hem. In his s udy, we used whole genome sequence le el da a om 4280 p ogeny es ed No dic Red Ca le bulls o scan he genome o loci a ec ing milk, a and p o ein yields. Resul s: Using a genome-wise signi icance h eshold, egions on Bos au us ch omosomes 5, 14, 23, 25 and 26 we e associa ed wi h a yield. Regions on ch omosomes 5, 14, 16, 19, 20 and 25 we e associa ed wi h milk yield and ch omosomes 5, 14 and 25 had egions associa ed wi h p o ein yield. Signi ican ly associa ed a ia ions we e ound in 227 genes o a yield, 72 genes o milk yield and 30 genes o p o ein yield. Ingenui y Pa hway Analysis was used o iden i y ne wo ks connec ing hese genes displaying signi ican hi s. When compa ed o p e iously mapped genomic egions associa ed wi h e ili y, signi ican ly associa ed a ia ions we e ound in 5 genes common o a yield and e ili y, hus linking hese wo ai s ia biological ne wo ks. Conclusion: This is he i s ime when whole genome sequence da a is u ilized o s udy genomic egions a ec ing milk p oduc ion in he No dic Red Ca le popula ion. Sequence le el da a o e s he possibili y o s udy quan i a i e ai s in de ail bu s ill canno unambiguously e eal which o he associa ed a ia ions is causa i e. Linkage disequilib ium c ea es di icul ies o pinpoin he causa i e genes and a ia ions. One solu ion o o e come hese di icul ies is he iden i ica ion o he unc ional gene ne wo ks and pa hways o e eal impo an in e ac ing genes as candida es o he obse ed e ec s. This in o ma ion on a ge genomic egions may be exploi ed o imp o e genomic p edic ion. Keywo ds: Milk ai s, No dic Red Ca le, Whole genome sequence, Associa ion s udy Backg ound The numbe o dai y cows in he No dic coun ies has been dec easing du ing he 21 s cen u y [1]. Howe e , o al milk p oduc ion le els ha e emained s able, as milk yield pe cow has inc eased. Fo example in Finland (including all dai y b eeds) he a e age p oduc ion pe cow pe yea has inc eased om 6786 l (2000) o 8201 l (2014), while a and milk con en s ha e emained ai ly cons an [2]. Global yea ly milk consump ion pe capi a is inc easing, and global demand o animal based oods is expec ed o be doubled by 2050 [3], d i en by bo h popula ion g ow h and inc eased consume p e e ences o mea and milk p oduc s. Ruminan s a e unique in hei capaci y o diges ib e and con e non-edible esou ces in o high quali y human nu i ion, making hem highly ele an o mee ing he inc easing global demand o ood. While animal b eede s ha e achie ed conside able imp o emen s in p oduc ion ai s, cow e ili y has been declining [4–6]. Howe e , du ing he ecen yea s, he dec ease in cow e - ili y in No dic coun ies has been slowing down and e en is e ac ed, due o he weigh ing o e ili y ai s in he b eeding p og am [7]. Many emale e ili y ai s in dai y ca le show an agonis ic gene ic co ela ions wi h milk p oduc ion ai s [8] bu wi h low o mode a e co ela- ions [9]. This implies ha simul aneous gene ic selec ion o inc eased milk yield and ep oduc i e pe o mance is * Co espondence: [email p o ec ed] 1 Animal Genomics, G een Technology, Na u al Resou ces Ins i u e Finland (Luke), Jokioinen, Finland Full lis o au ho in o ma ion is a ailable a he end o he a icle © 2016 Iso-Tou u e al. Open Access This a icle is dis ibu ed unde he e ms o he C ea i e Commons A ibu ion 4.0 In e na ional License (h p://c ea i ecommons.o g/licenses/by/4.0/), which pe mi s un es ic ed use, dis ibu ion, and ep oduc ion in any medium, p o ided you gi e app op ia e c edi o he o iginal au ho (s) and he sou ce, p o ide a link o he C ea i e Commons license, and indica e i changes we e made. The C ea i e Commons Public Domain Dedica ion wai e (h p://c ea i ecommons.o g/publicdomain/ze o/1.0/) applies o he da a made a ailable in his a icle, unless o he wise s a ed. Iso-Tou u e al. BMC Gene ics (2016) 17:55 DOI 10.1186/s12863-016-0363-8 possible [9]. Simul aneous b eeding o bo h p oduc i e and e ile cows would bene i subs an ially om knowing he gene ic and physiological links be ween p oduc ion and heal h o disen angle he e ec s on hese ai s. Recen esul s in Hols ein and Je sey b eeds indica e li le o no o e lap be ween genomic egions associa ed wi h milk yield and e ili y [10, 11]. Genome wide associa ion s udies (GWAS) ha e bene- i ed om he apid de elopmen o single nucleo ide polymo phism (SNP) geno yping echnologies, bu des- pi e o he ela i ely high densi y o he a ailable SNP chips, inding he causa i e mu a ion is no s aigh o - wa d. The high le el o linkage disequilib ium in dai y ca le esul s in long quan i a i e ai loci (QTL) egions wi h se e al possible candida e genes. Using whole genome le el sequence a ian s o associa ion analyses would be an ul ima e choice, because hen he causa i e a ian is mos likely included among he s udied a ian s. Po en ially his helps o pinpoin he causa i e mu a ions hus leading o a be e unde s anding o bio- logical mechanisms behind he QTL [12] and imp o e he e iciency o genomic selec ion [13]. Using sequence le el SNPs will also enable iden i ica ion o SNPs ha explain a small ac ion o he ai a ia ion because ei he he causal SNP and/o SNP(s) wi h high linkage disequilib ium (LD) wi h he causal a ian a e included in he analysis [13]. His o ically sepa a ed h ee dai y b eeds Finnish Ay shi e om Finland, Danish Red om Denma k and Swedish Red om Sweden a e a p esen unde a join b eeding alue es ima ion sys em, known as he No dic Ca le Gene ic E alua ion [14]. P e ious QTL s udies o milk ai s in No dic Red Ca le (NRC) ha e been done wi hin he subpopula ions wi h mic osa elli e ma ke s and ai ly small sample sizes (e.g. [15, 16]). The objec i e o his s udy was o use a ia ions a he genome sequence le el o ca y ou associa ion s udy o milk, a and p o ein yields in NRC; o iden i y po en ial causal a ian s and unde s and he gene ic a chi ec u e o hese ai s. In addi ion, he da a p o ides he possibili y o compa e he esul s o simila s udies o e ili y ai s in he NRC [17], o e eal po en ial QTL wi h an agonis- ic e ec s o milk p oduc ion and e ili y ai s. Me hods No animal expe imen s we e pe o med in his s udy, and, he e o e, app o al om he e hics commi ee was no equi ed. Semen samples we e collec ed o b eeding pu poses by local o ganiza ions wi h app op ia e pe mi s. Milk, a and p o ein yields’ ai de ini ions a e s an- da dized ac oss he No dic coun ies. Pheno ypic eco ds o dai y ca le a e housed in a cen alized da abase [14]. B eeding alues o milk, a and p o ein yield (MY, FY and PY) a e based on p oduc ion igu es exp essed in kilog ams aken om ou ine milk eco ds and hen com- bined in o an index o each ai . Fo de ails on gene ic e alua ion o milk yield ai s in No dic coun ies see [18]. The b eeding alues used o associa ion analysis we e de- eg essed b eeding alues om he ou ine gen- e ic e alua ion by NAV (No dic ca le gene ic e alua ion) and we e a ailable o 4280 p ogeny es ed NRC bulls (2127 om Finland, 1217 om Sweden, 915 om Denma k and 21 om o he coun ies). The eliabili ies o he de eg essed b eeding alues we e in he ange o 0.67 o 0.99 wi h a mean o 0.95 and he i s qua ile a 0.94. SNP a ay geno yping All 4280 NRC bulls wi h de eg essed b eeding alues we e geno yped using Bo ineSNP50 BeadChip SNP a ay e sion 1 o 2 (Illumina Inc., San Diego, CA). DNA was ex ac ed using s anda d p ocedu es om semen samples. Chip ypings we e done by GenoSkan A/S, Tjele, Denma k o labs belonging o Aa hus Uni e si y. The quali y pa ame e s used o selec ion o SNPs we e minimum call a es o 85 % o indi iduals and 95 % o loci. Ma ke loci wi h mino allele equencies below 5 % and de ia ion om Ha dy-Weinbe g p opo ion (P< 0.00001) we e excluded. The minimal accep able GC sco e was 0.60 o indi idual ypings, and indi iduals wi h a e age GC sco es below 0.65 we e excluded. The numbe o SNP emaining a e quali y con ol was 43,415 in he geno ypes ob ained om Bo ineSNP50 BeadChip SNP a ay (50 K da a se ). The genome posi ions o he SNPs we e acco ding o he UMD3.1 Bo ine genome assembly [19]. Impu a ion o whole genome sequences The 50 K geno ypes o hese bulls we e impu ed o whole genome sequence da a using a wo-s ep app oach [20]. Geno ypes om 50 K chip o each bull we e i s impu ed o a high-densi y SNP a ay (HD) using a mul i-b eed e e ence o 3383 animals (1222 Hols ein, 1326 NRC and 835 Danish Je sey indi iduals) which had been geno yped wi h he Illumina Bo ineHD chip (Illumina Inc., San Diego, CA). The numbe o SNPs, a e impu a ion o he Bo ineHD chip, was 648,219. These impu ed HD geno ypes we e subsequen ly impu ed o he whole genome sequence le el using a mul i-b eed e e ence panel o 1228 animals om Run4 o he 1000 Bull Genomes P ojec [13, 21] and addi ional whole genome sequences om Aa hus Uni e si y [22] including 368 Hols ein, 86 RDC, 88 Je sey and es om numbe o ca le b eeds. Da ase s (SNP a ay ypes and whole sequence) we e p e-phased wi h BEAGLE4 1274 [23] and geno ype impu a ion we e done using Minimac2 [24]. The impu a ion accu acy o his da a is epo ed ea lie [25], bu wi h a smalle whole genome sequence e e - ence popula ion. Sequence a ian s ha ing impu a ion accu acy 2 ( a io o empi ically obse ed a iance o Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 2 o 12 he allele dosages o he expec ed binomial a iance a Ha dy-Weinbe g equilib ium and was ob ained om Minimac2 so wa e ou pu ) less han 0.5 we e il e ed away. The mean accu acy o he a ian s wi h 2 >0.5 was 0.94. Associa ion analysis The associa ion analysis o each o he impu ed sequence a ian s (mino allele equency, MAF > 0.005 and de i- a ion om Ha dy-Weinbe g p opo ion > 0.00001) was ca ied ou using a wo-s ep a iance componen s-based app oach o accoun o popula ion s a i ica ion imple- men ed in he EMMAX so wa e ool [26]. In a i s s ep, he polygenic and e o a iances a e es ima ed using ollowing a iance componen model: y¼1μþaþe whe e yis a ec o o de- eg essed b eeding alues, 1is a ec o o ones, μis he in e cep , Gis he kinship ma ix buil based on high-densi y SNP geno ypes using EMMAX so wa e, ais a ec o o b eeding alues assumed o ha e a mul i a ia e no mal dis ibu ion a~N(0,Gσ a 2 ), eis a ec o o andom esiduals assumed o ha e a mul i a ia e no mal dis ibu ion e~N(0,Iσ e 2 ), whe e Iis an iden i y ma ix, σ a 2 is he addi i e gene ic a iance and σ e 2 is he e o a iance. In a second s ep, he SNP e ec is ob ained using a linea eg ession model: y¼1μþxb þη; whe e xis a ec o o impu ed geno ype dosages ( anged be ween 0 and 2), bis he allele subs i u ion e ec and ηis a ec o o andom esidual de ia es wi h (co) a iance s uc u e Gσ a 2 +Iσ e 2 . Sea ch o mul iple QTL in a genomic egion To es i mul iple QTL a e seg ega ing in a genomic egion we included he mos signi ican o he known causal a ian as co ac o in he model and check o addi ional QTL in a genomic egion o a yield on ch omosomes 14, 25 and 26 and o milk yield on ch omosome 14. We i ed he SNPs (Addi ional ile 1) as ixed e ec o a linea mixed model. The s a is ical model is desc ibed by he o mula: y¼1μþqsnp op þxgþZu þe whe e y,1,μa e desc ibed as in he EMMAX model, snp op is he e ec o he SNP i ed as co- ac o in he model, and gis headdi i egene ice ec o he heSNP unde s udy, qand xa e ec o s o SNP geno ype dosages ( anging om 0 o 2), and uis a ec o o andom poly- genic e ec s, which a e no mally dis ibu ed u~ N(0, Aσ u 2 ), whe e Ais he pedig ee-based addi i e ela ionship ma ix, σ u 2 is he polygenic a iance, Zis an incidence ma ix ela ing pheno ypes o he co esponding andom poly- genic e ec s, and eis a ec o o esidual e ec s, which a e no mally dis ibu ed e~N(0,Dσ e 2 ), whe e Dis a diagonal ma ix wi h elemen s d ii =(1− DRP 2 )/ DRP 2 o accoun o he e ogeneous esidual a iances due o di e en eliabil- i ies o DRP ( DRP 2 ), and σ e 2 is he esidual a iance. Analyses we e pe o med using he DMU package [27]. Signi i- cance es ing o SNP e ec s was pe o med using a wo-sided - es . The null hypo hesis was g = 0. A e ha , he Bon e oni co ec ion was applied same as in he EMMAX analysis o con ol o alse posi i e associa ions. The genome-wise signi icance h eshold co esponding o an e o a e o 0.05 was se a 3.16 × 10 −9 a e co ec- ion o mul iple es ing using a Bon e oni co ec ion o 15,679,852, 15,679,853 and 15,679,844 independen es s o a , milk and p o ein yield espec i ely. Only SNPs wi h he p alue less han 3.16 × 10 −9 (−log 10 (p)≥8.50) we e anno a ed wi h he a ian e ec p edic o (VEP) ool using he Ensembl da abase, Release 82 [28]. The p edic ion whe he an amino acid subs i u ion caused by missense a ia ion a ec s p o ein unc ion was es ima ed by SIFT analysis [29] implemen ed in VEP ool [28]. The SIFT p edic ion is based on sequence homology and he physical p ope ies o amino acids. Manha an plo s we e c ea ed wi h he qqman .0.1.2 R package [30]. In addi ion, we compa ed ou indings o esul s ob ained om s udy by Höglund e al. [17] whe e a simila genome-wide associa ion s udy o emale e ili y in No dic Red ca le was conduc ed. Ingenui y pa hway analysis and en ichmen analysis Lis s o genes wi h signi ican hi s om he QTL peak egions associa ed wi h each milk ai we e uploaded in o he Qiagen’s Ingenui y® Pa hway Analysis IPA® [31]. Fo his pu pose, also SNPs signi ican ly associa ed wi h Fe ili y index (FI) in he No dic Red Ca le [17] we e anno a ed wi h he VEP ool [28] and genes ha ing one o mo e signi ican SNP we e analyzed wi h IPA® [31]. Bioma ool [32] embedded in Ensembl da abase [33], was used o sea ching human homologs o cow Ensembl IDs o he genes. In case he e was mo e han one, all epo ed human o hologs we e kep . The human homologue lis s o each ai included 214, 69, 66 and 263 genes o a , milk, p o ein yields and e ili y index, espec i ely. A e unning co e analysis o each ai , ne wo ks based on he in o ma ion o gene connec i i y in Ingenui y Knowledge Da abase wi h highes sco e- alues we e conside ed. Sco e- alue ep esen s he nega i e log o he p- alue o he likeli- hood ha he molecules would be ound oge he by chance. Gene on ology (GO) e m en ichmen analysis wi h genes ound wi hin he op SNPs was pe o med wi h a Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 3 o 12 Singula En ichmen Analysis (SEA, Fishe ’s exac es , FDR < 0.01) p o ided by Ag iGO webpage [34]. Resul s Depending on he ai we could iden i y se e al hou- sands (3594, FY), less han a housand (755, MY) o less han a hund ed (85, PY) signi ican ly associa ed SNPs (Bon e oni co ec ed h eshold o signi icance -log 10 (p)≥8.50; Addi ional iles 2, 3 and 4). No common signi ican ly associa ed SNPs we e ound be ween his s udy and wi h hose ound o emale e ili y ai s [17]. Howe e , signi ican ly associa ed SNPs we e ound in i e common genes be ween e ili y and a yield. These i e genes a e loca ed on ch omosomes 25 (ENSBTAG00000034643) and 26 (GBF1,TMEM180, ACTR1B, and b a-mi -146b). The summa y o he anno a ions o signi ican ly associ- a ed SNPs a e p esen ed in Table 1. Among anno a ed SNPs, in on a ian s a e he mos common ype o each o he ai s. SNPs changing amino acid (missense a ia- ions) a e a e, 21 o FY, 1 o MY, 3 o PY and 18 sha ed be ween FY and MY. Se en missense a ia ions we e p edic ed by SIFT analysis [29] o be dele e ious, i.e. po en ially leading o changes in he unc ion o he p o ein (Addi ional ile 5). Splice egion a ian s a e a e o ming only one pe cen o less om he o al amoun o signi ican ly associa ed SNPs (Addi ional ile 5). Resul s o a yield wi hin he la ge associa ed a eas on ch omosomes 14, 25 and 26 we e u he examined by ixing he e ec o op SNP(s). Only peaks ha emained a e ixing he o he op SNPs we e consid- e ed as po en ial QTL in hese ch omosomes. Se en, eigh and ou sepa a e QTL egions (“peaks”) we e de ined o a , milk and p o ein yield, espec i ely (Fig. 1, Addi ional iles 6, 7, 8, 9, 10, 11, 12 and 13). The peaks we e de ined as con inuous egions con aining SNPs ha ing –log 10 (p)≥8.50. Top SNPs wi h he conse- quences o each de ined QTL egion pe ai a e lis ed in Table 2. The highes peak was obse ed on Bos au us ch omosome 14 (BTA14) (Addi ional ile 7) spanning he egion om 1,448,510 bp o 2,271,832 bp o a yield (509 SNPs ha ing –log 10 (p)≥8.50) and om 1,448,510 bp o 2,271,832 bp o milk yield (455 SNPs ha ing –log 10 (p)≥8.50). The highes –log 10 (p) alues wi hin hese egions we e ob ained o SNPs s136783505 (bp 1,807,140) o a yield and s133033480 (bp 1,743,939) o milk yield. The highes peak o p o ein yield was loca ed on BTA25 (Addi ional ile 11) wi hin he egion 3,306,363–3,516,671 which con ained 40 SNPs, he highes p alue being o SNP s110749311 (bp 3,498,960). Fa yield We iden i ied se en di e en QTL egions on i e di e en ch omosomes a ec ing a yield (Table 2, Fig. 1, Addi ional ile 2). The s onges associa ion ound o a yield was lo- ca ed on BTA14 (Addi ional ile 7) in he DGAT1 (Diacyl- glyce ol O-acyl ans e ase 1) gene egion. In ou da a he s onges associa ion is o he a ia ion s136783505 (bp 1,807,140) wi h no unc ional anno a ion, loca ed 2578 bp downs eam o he DGAT1 gene (Table 2). Howe e , se - e al o he a ian s loca ed nea by, including he p e iously iden i ied causa i e a ian K232A [35] a bp 1,802,266, show simila ly high -log 10 (p) (Addi ional ile 2). To in es i- ga e he signi icance o o he SNPs in he egion, we i ed he a ia ion K232A as ixed e ec . None o he o he SNPs emained signi ican (Addi ional ile 14) a e he ixa ion o he K232A a ia ion. O he SNPs wi h s ong associa ions wi h a yield we e ound on BTA5 (Addi ional ile 6), BTA23 (Addi ional ile 11), BTA25 (Addi ional ile 12) and BTA26 (Addi ional ile 13). The QTL egion on BTA5 is loca ed be ween 92,372,732 bp and 94,425,668 bp and he a ia ion wi h he s onges associa ion ( s209818856, pos. 93,945,694) is loca ed in an in on o he gene MGST1 (Table 2). On BTA23 he associa ion signal o a yield comes om he egion 28,567,796–28,591,530, he op a ia ion being lo- ca ed a bp 28,567,796 in he in on o gene TRIM26 (Table 2). On BTA25 and BTA26, complex pa e ns o associ- a ion we e seen (Addi ional iles 12 and 13). To cla i y he numbe o independen QTL wi hin hese egions we in es iga ed he signi icance o he SNPs by i ing Table 1 The numbe o signi ican SNPs (- log 10 (p)≥8.50) o each ai and how SNPs a e di ided in o di e en consequences. SNPs we e anno a ed wi h he a ian e ec p edic o – ool [28]. One SNP can ha e mo e han one anno a ion T ai Numbe o signi ican SNPs In on a ian In e genic a ian Downs eam gene a ian Ups eam gene a ian Synonymous a ian Missense a ian 3′UTR a ian Splice egion a ian , in on a ian 5′UTR a ian Splice egion a ian , synonymous a ian Fa yield 3594 1641 1307 600 551 94 40 39 9 6 2 Milk yield 755 322 195 301 240 42 20 16 5 1 2 P o . yield 85 50 14 32 12 3 3 1 1 Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 4 o 12 he op SNPs ( i e o BTA25 and eigh o BTA26, Addi ional ile 1) as ixed e ec s alone and in di e en combina ions. On BTA 25, signi ican associa ions emained a bp 9,870,005 (in onic egion o he CLEC16A gene) and a bp 36,226,978 (in e genic egion) (Table 2). Two QTL emained also on BTA26, one in he egion o he NEURL1 gene (posi ion 24,379,571) and he o he in an in e genic egion ( op SNP a bp 44,802,991) (Table 2). Milk yield In all, eigh QTL egions we e ound o milk yield (Table 2, Fig. 1, Addi ional ile 3). They we e loca ed on six di e en ch omosomes (BTA5, BTA14, BTA16, BTA19, BTA20 and BTA25). The s onges associa ion wi h milk yield was ound on BTA14 (1,448,510– 2,271,832) ha ing op a ia ion loca ed a bp 1,743,939. This a ia ion is wi hin wo o e lapping genes, CPSF1 and ADCK5.TheQTL egionis hesameaswas ound associa ed wi h a yield bu wi h di e en op a i- a ion. As o a yield, he signi icance o o he han he p e iously known causa i e DGAT1 a ia ion [35] was es ed by i ing he a ia ion K232A as a ixed e ec . None o he o he SNPs emained signi ican a e ixing he K232A e ec (Addi ional ile 15). The weakes signi ican associa ion was loca ed on BTA5 wi h op a ia ion s383553819 (posi ion 112,343,204) loca edinanin ono hegeneMKL1 (Addi ional ile 6). On BTA19 he associa ed egion con ained no anno a ed genes (61,447,138–61,449,096), and he op a ia ion s210324693 (bp 61,449,096 bp) was loca ed in an in e - genic egion (Table 2, Addi ional ile 9). Two QTL egions we e de ec ed on BTA20 (Addi ional ile 10). The QTL we e loca ed in he egions 30,531,217–32,952,019 and 37,766,226–39,183,141. The known causa i e a ia ion F279Y (bp. 31,909,478) o milk ai s in he gene GHR [36]was he opSNPinou analysis o heQTLlo- ca ed in he i s egion. The op a ia ion wi hin he second QTL was loca ed a bp 38,828,254 in he in e genic egion, bu his QTL egion also includes he PRLR gene p e iously indica ed o be linked o milk p oduc ion (e.g. [37]). BTA25 ha bo s wo QTL, he i s a bp 2,669,704 and he second in he egion om 3,494,706 bp o 3,516,671 bp, he op SNP lo- ca ed a bp 3,498,960 downs eam om he gene PAM16 (Addi ional ile 12). Fig. 1 Genome-wide Manha an plo s o he a yield (FY), milk yield (MY) and p o ein yield (PY). Red line indica es he genome-wide signi icance le el (−log 10 (p) = 8.50) Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 5 o 12 Table 2 QTL egions o each o he ai . Top SNP o each QTL a e shown including posi ion, -log 10 (p)- alues, mino allele equency (MAF), gene in o ma ion, anno a ion o he op SNP, allele subs i u ion e ec (b. alue) and s anda d e o o b. alue (SE) CHR S a (bp) End (bp) Leng h o he QTL egion (bp) Signi ican SNPs in egion No o genes wi h signi ican SNPs wi hin he QTL egion Top SNP Posi ion o he op SNP (bp) -log 10 (p) MAF Gene Anno a ion o he op SNP b. alue SE Fa yield 5 92,372,732 94,425,668 2,052,936 330 3 s209818856 93,945,694 27.49 0.38 MGST1 in on a ian −2.807 0.253 14 1,448,510 2271832 823,322 509 49 s136783505 1,807,140 42.01 0.07 DGAT1/HSF1 downs eam a ian / in on a ian −6.709 0.484 23 28,567,796 28,591,530 23,734 5 1 s381390819 28,567,796 9.36 0.45 TRIM26 in on a ian 1.148 0.184 25 8,222,347 11,507,986 3,285,639 883 16 s379546164 9,870,005 14.40 0.24 CLEC16A in on a ian −1.769 0.224 25 36,226,978 36,227,132 154 2 s109480808 36,226,978 8.72 0.16 in e genic a ian 1.635 0.272 26 22,144,777 24,793,744 2,648,967 595 32 s438420348 24,379,571 14.28 0.16 NEURL1 in on a ian −2.273 0.290 26 44,802,991 44,802,991 0 1 s135624939 44,802,991 9.72 0.28 in e genic a ian −1.474 0.231 Milk yield 5 112,343,204 11,2450,860 107,656 2 1 s383553819 112,343,204 8.70 0.36 MKL1 in on a ian 1.437 0.239 14 1,448,510 2,271,832 823,322 455 48 s133033480 1,743,939 33.01 0.09 CPSF1/ADCK5 downs eam a ian / splice egion a ian , in on a ian 6.266 0.513 16 1,322,611 1,322,611 0 1 1 s108979795 1,322,611 8.58 0.28 LAX1 ups eam a ian −1.451 0.243 19 61,447,138 61,449,096 1958 5 s210324693 61,449,096 8.93 0.32 in e genic a ian 1.490 0.244 20 30,531,217 32,952,019 2,420,802 74 5 s385640152 31,909,478 15.56 0.11 GHR missense a ian −3.877 0.472 20 37,766,226 39,183,141 1,416,915 34 3 NA 38,828,254 9.54 0.16 in e genic a ian −2.250 0.356 25 2,669,704 2,669,704 0 1 s209691835 2,669,704 9.21 0.13 in e genic a ian −2.855 0.460 25 3,494,706 3,516,671 21,965 13 4 s110749311 3,498,960 9.32 0.41 PAM16/GLIS2 downs eam a ian 1.225 0.196 P o ein yield 5 112,450,860 112,450,860 0 1 s109041054 112,450,860 8.88 0.48 in e genic a ian −1.473 0.242 14 1,802,667 1,802,667 0 1 2 NA 1,802,667 8.52 0.06 DGAT1/HSF1 in on a ian / downs eam a ian 3.354 0.564 25 1,094,996 1,257,612 162,616 12 3 s136085792 1,103,856 10.83 0.22 UNKL in on a ian 1.694 0.250 25 3,306,363 3,516,671 210,308 40 8 s110749311 3,498,960 11.70 0.41 PAM16/GLIS2 downs eam a ian 1.427 0.202 Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 6 o 12 P o ein yield In all, p o ein yield did no show as many signi ican ly associa ed SNPs as obse ed o a and milk yield (Table 2, Fig. 1, Addi ional ile 4). Ch omosomes ha ing QTL o p o ein yield we e BTA5, BTA14, and BTA25. Fo bo h BTA5 and BTA14, only one a ian o each eached he signi icance cu -o le el; on BTA5 an in e - genic a ia ion a bp 112,450,860 (Addi ional ile 6) close o he a ian ound o milk yield and on BTA14, he a ia ion (bp 1,802,667) loca ed in an in on o he DGAT1 gene (Addi ional ile 7). The wo o he QTL o p o ein yield we e bo h lo- ca ed on BTA25. The QTL on BTA25 a 3,306,363– 3,516,671 o e lapped wi h he QTL ound o milk yieldanddisplayed hesame opSNP(bp3,498,960). The o he QTL on BTA25 was unique o p o ein yield, op SNP (bp 1,103,856) loca ed in he gene UNKL. Ne wo ks o associa ed genes and en ichmen analysis The ne wo ks wi h highes sco es o each ai a e p esen ed in Addi ional iles 16, 17 and 18. The wo op ne wo ks om genes associa ed wi h a yield (Addi ional ile 16) had sco es o 49, he ne wo k associ- a ed wi h ca bohyd a e me abolism, gene exp ession and lipid me abolism is p esen ed in Addi ional ile 16a. Addi ional ile 16b shows he ne wo k gene a ed om he genes associa ed o e ili y index. The sco e alue o his ne wo k is 41 and i consis s o genes associa ed wi h in lamma o y esponse, cell- o-cell signaling and lymphoid issue s uc u e and de elopmen . E en hough common signi ican ly associa ed SNPs we e no ound be ween his s udy and ha o [17] on emale e ili y ai s, signi ican ly associa ed SNPs we e ound in i e common genes. Two o hem (GBF1 and b a-mi -146b) a e p esen in bo h a and e ili y ne wo ks (Addi ional ile 16a and b). The op ne wo k o milk yield (sco e 43, Addi ional ile 17) was associa ed wi h unc ions molecula anspo , o gan mo phology and o ganismal de elopmen . I includes he known milk candida e genes DGAT1,GHR and PRLR,as well as he op hi s om BTA5 (MLK1) and BTA16 (LAX1). The op ne wo k o p o ein yield (sco e 35, Addi ional ile 18) is connec ed o unc ions cell dea h and su i al, cance and o ganismal inju y. Al oge he 18, 18, 11 and 16 GO e ms we e signi i- can ly ( alse disco e y a e, FDR < 0.01) en iched o FY, MY, PY and FI, espec i ely. A b oad GO e m, mul i- cellula o ganismal p ocess, was he mos signi ican o all ou ai s, o ally eigh e ms we e sha ed be ween hem. Fo example all ai s we e ha ing QTL in he egions con aining signi ican en ichmen o genes ela ed o ep oduc ion and ep oduc i e p o- cesses (Table 3). Discussion A la ge numbe o a ian s we e ound signi ican ly associa ed wi h milk, a and p o ein yield. Ou indings suppo p e ious QTL indings om he No dic Red b eeds, e.g. on BTA 5, 14, 20 [37, 38] and loca e new a ia ions ha a e good candida es o be causa i e a ia ions. This is he i s ime when NRC popula ion is s udied wi h impu ed whole genome sequence a ian s in o de o e ine QTL associa ed wi h milk p oduc ion. Howe e , i is s ill di icul o pinpoin he causa i e a ian among se e al closely linked, almos equally signi ican ly associa ed a ia ions. One way o classi y he a ia ions is o look a he p edic ed unc ional con- sequence o he SNP [29]. The possibili y ha a a ia ion has an impac on he pheno ype is highe i he a ia ion causes an amino acid change (missense a ia ion) which is p edic ed (e.g. wi h SIFT analysis, [29]) o ha e an e ec o p o ein unc ion, is loca ed on splicing si e, o is loca ed downs eam o ups eam o he known gene (possible egula o y egions o he ansc ip ion). On he o he hand, genome anno a ion o ca le is s ill incom- ple e and mos egula o y elemen s emain unknown. In he sea ch o biologically ele an ma ke s he in o ma- ion o in e ac ions be ween genes in known pa hways o ne wo ks can be use ul. In his s udy, we iden i ied some in e es ing gene in e ac ion ne wo ks based on he signi i- can ly associa ed a ian s wi hin genes (e en hough he unc ional e ec s o he a ian s could no be p edic ed). The esul smaybeused oha eaclose looka alsoo he genes in he indica ed ne wo ks o unc ional a ian s. Al hough no common SNPs we e ound associa ed wi h milk p oduc ion ai s and e ili y, he i e common genes be ween a yield and e ili y gi e some indica ion o he ela ionship be ween hose ai s. The genes b a-mi -146b and GFB1 a e associa ed wi h a e ili y gene ne wo k linked wi h in lamma o y esponse and cell- o-cell signal- ing and he a yield ne wo k connec ed wi h lipid and ca bohyd a e me abolism. Fu he suppo was gained om he gene en ichmen analysis, bo h he ai s show signi ican en ichmen o he genes ela ed o example o ep oduc ion and ep oduc ion and ep oduc i e p o- cesses, al oge he ha ing 16 common GO e ms. F om he ou ch omosomes epo ed o ha bo highes numbe o QTL o milk p oduc ion [39], wo we e indi- ca ed by ou da a (BTA14 and BTA20). S ucken e al. [40] summa ized 14 genes om en di e en ch omosomes o be he majo genes in ol ed in milk p oduc ion. Among hose genes a e DGAT1 and GHR. Some commonly ound QTL (e.g. BTA6, [41]) we e no seen in ou da a; ha could be due o ixa ion o he QTL o e y low MAF in he NRC popula ion. One explana ion could be ha he EMMAX me hod chosen o associa ion analysis migh be oo conse a i e. EMMAX uses app oxima ions o con- s uc ing es o he ixed SNP e ec s o in e es in he Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 7 o 12 linea mixed model because i ing a ull linea mixed model o each SNP in u n ac oss he genome is compu a- ionally challenging [42]. This leads o sys ema ic unde - es ima ion o he mos signi ican p alues [43], bu makes EMMAX one o he as es LMM based p og ams [42]; an a gumen ha has o be conside ed when ha ing whole genome sequence le el da a om subs an ially la ge amoun o indi iduals. DGAT1 (BTA14) The s onges signal o associa ion was ound om BTA14 o a and milk yield. Also p o ein yield is signi ican ly associa ed wi h he same egion on BTA14. The op SNP a ied depending on he ai (see Table 2). Dae wyle e al. [13] analyzed associ- a ion o ea ly lac a ion milk a pe cen age wi h whole genome sequence a ia ion da a in Fleck ieh and Hols ein bulls. As in ou s udy, he p e iously epo ed causa i e a ia ion K232A (bp 1,802,266, [35]) o DGAT1 gene was no he a ian wi h he lowes p alue in Hol- s ein and Fleck ieh o milk a p oduc ion al hough K232A was among he op SNPs. Pausch e al. [44] used whole genome sequence da a o impu e Ge man Fleck ieh and Hols ein-F iesian ca le geno ypes om a la ge se o animals candida e egions and we e able o con i m he associa ion o K232A mu a ion, howe e , hei app oach was biased as only known candida e SNPs we e es ed o associa ion. In ou da a, close inspec ion o he associa ed a ia ions (−log 10 (p)≥8.50, a yield) in he DGAT1 egion e ealed ha when he e ec o he K232A mu a ion was ixed, no addi ional s a is ically signi ican SNP e ec s we e le . The causa i e mu a ion (K232A/ s109326954) o DGAT1 was epo ed al eady o e a decade ago [45] and has been unc ionally con i med [35]. The K allele inc eases milk a pe cen age [35], whe eas allele A inc eases milk p oduc ion [46]. The e a e di e en possible explana ions why K232A did no u n ou o be he mos signi ican ly associa ed a ia ion in his s udy. Impe ec impu a ion may a ec associa ion esul s. Accu acy o he impu a ion is conside - ably imp o ed by inc easing he size o e e ence panel, i.e. sequenced animals [47] and impu a ion accu acy seems o be highe when popula ions unde s udy a e combined o he impu a ion p ocesses [13]. Ou e e ence panel consis ed o a mul i-b eed popula ion wi h 1228 indi- iduals om se e al b eeds including bo h dai y and bee ca le. The DGAT1 egion would be an in e es - ing candida e o s udy wi h he in o ma ion om 1000 Bulls Genomes P ojec [13, 21]. I would gi e a chance o s udy he haplo ype s uc u e o he egion in he ca le popula ion wo ldwide and possibly ace back he e olu ion o he QTL e ec . Table 3 GO en ichmen e ms ha ing FDR < 0.01 om he genes ha ing signi ican a ia ions (−log 10 (p)≥8.50) o a yield (FY), milk yield (MY), p o ein yield (PY) and e ili y index (FI) FY MY PY FI GO e m Desc ip ion FDR FDR FDR FDR GO:0022610 Biological adhesion 6.20E-07 0.0013 2.00E-06 GO:0065007 Biological egula ion 9.40E-12 6.80E-05 0.072 2.40E-28 GO:0044085 Cellula componen biogenesis 3.60E-09 0.048 1.00E-22 GO:0016043 Cellula componen o ganiza ion 1.00E-37 2.90E-12 0.00027 1.80E-76 GO:0009987 Cellula p ocess 4.90E-14 4.60E-05 0.0074 3.80E-20 GO:0016265 Dea h 7.80E-21 2.00E-07 4.00E-25 GO:0032502 De elopmen al p ocess 8.80E-102 2.10E-41 1.40E-14 GO:0051234 Es ablishmen o localiza ion 1.80E-14 2.80E-05 0.0057 1.70E-27 GO:0040007 G ow h 5.60E-42 2.10E-17 9.40E-45 GO:0002376 Immune sys em p ocess 0.00022 0.0032 9.80E-10 GO:0051179 Localiza ion 5.90E-20 1.10E-07 0.0059 4.70E-44 GO:0008152 Me abolic p ocess 9.40E-12 0.00035 0.57 7.80E-18 GO:0032501 Mul icellula o ganismal p ocess 8.00E-126 1.10E-44 1.50E-20 1.20E-161 GO:0048519 Nega i e egula ion o biological p ocess 2.80E-50 1.30E-25 1.20E-07 GO:0048518 Posi i e egula ion o biological p ocess 1.70E-32 1.60E-08 GO:0050789 Regula ion o biological p ocess 7.10E-10 0.00011 0.32 6.90E-22 GO:0000003 Rep oduc ion 3.70E-57 3.70E-14 3.50E-09 1.20E-54 GO:0022414 Rep oduc i e p ocess 1.70E-50 7.50E-11 1.90E-07 6.30E-38 GO:0050896 Response o s imulus 6.00E-31 5.50E-06 8.80E-07 8.60E-44 Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 8 o 12 MGST1 and MKL1 (BTA5) Vii ala e al. [16] showed ha Finnish Ay shi e has a milk p oduc ion QTL a he p oximal end o BTA5. Wang e al. [48] epo ed o a QTL o milk a pe cen age in Ge man Hols ein-F iesian a he loca ion o 94,551,792 bp and sugges ed a candida e gene o be EPS8. One o he mos signi ican SNP in he s udy by Aliloo e al. [11] o milk yield (Je sey and Hols ein) was loca ed a bp 94,518,850 on BTA5. Fo a yield, we obse ed an associ- a ion peak a bp 93,945,694 (Table 2). This a ia ion is in in onic egion o he gene MGST1 ha ing a ole in oxida i e s ess eac ion [49] QTL peak o he milk yield on BTA5 is loca ed wi hin an in on o he MKL1gene (bp 112,343,204). MKL1 is ela ed o an- sc ip ion egula ion. P o ein yield associa ion peak is loca ed a bp 112,450,860 on BTA5 and is anno a ed o he non-coding egion. P e iously QTL ela ed o body weigh ha e been mapped nea by [50]. GHR and PRLR (BTA20) BTA20 is among he ch omosomes ha bo ing many QTL ela ed o milk p oduc ion [39]. GHR has been epo ed as one o he majo genes in ol ed in milk p oduc ion [40]. Fi s Blo e al. [36] ound ha a ia ion F279Y (bp 31,909,478) in he GHR gene is associa ed wi h a s ong e ec on milk yield and composi ion and o he s udies ha e con i med i (e.g. [44]). O he a ia ions han F279Y ha e also been ound o a ec milk p oduc ion nea by he GHR gene [38]. In ou s udy, F279Y has he s onges associa ion o milk yield among he NRC popula ion on BTA20. Fu he con i ma ion o he causali y o he F279Y comes om SIFT analysis p edic ing he mu a ion o be dele e ious o p o ein s uc u e, i.e. po en ially al e - ing he p o ein s uc u e hus possible leading o changes in unc ion o he p o ein. The op SNP (bp 31,909,478) in he GHR gene egion is clea ly he mos signi ican one (Addi ional ile 10) in con as o DGAT1 egion whe e se e al a ia ions a e s ongly associa ed o milk yield (Addi ional ile 7). The Y allele is p edic ed o be un a o - able o milk p oduc ion [36], bu i is s ill ai ly common in he NRC popula ion. GHR has been sugges ed o be unde balancing selec ion because o he obse ed high a ia ion in he cy oplasmic egion [51]. Ano he gene on BTA20 o special in e es is PRLR and he a ia ion S18N (posi ions 39,115,344-39,115,345) [37]. Howe e , i has been sugges ed ha S18N is a he linked o he causa i e mu a ion han being causa i e i sel [44]. We ound ha he a ia ion a bp 38,828,254 on BTA20 loca ed in he in e genic egion, was indica ed o be he mos likely candida e esponsible o he QTL e ec seen. This a ia ion is loca ed app oxima ely 245,000 base pai s downs eam o he PRLR gene and addi ional s udies a e equi ed o esol e he mechanisms how i may in luence milk p oduc ion. TRIM26 (BTA23) An in onic a ian in he gene TRIM26 a bp 28,567,796 on BTA23 has an associa ion wi h a yield. The unc ion o he TRIM26 gene, a membe o he ipa i e mo i (TRIM) gene amily, is unknown [52]. I is loca ed close o he majo his ocompa ibili y complex (MHC) class I egion. Feed in ake QTL ha e been mapped close o he associa ion peak obse ed in his s udy [53]. PAM16,UNKL and CLEC16A (BTA25) Al oge he six associa ion peaks (QTL) we e obse ed om BTA25 o di e en ai s. The same a ia ion (bp 3,498,960) in gene PAM16 is associa ed wi h bo h milk and p o ein yield. The gene has a c i ical ole in p o ein ansloca ion ac oss he inne mi ochond ial memb ane [54]. O he QTL on BTA25 we e ound o milk yield in he in e genic egion (bp 2,669,704) and o p o einyieldQTLa bp1,103,856in hegene UNKL ha has a ole on p o ein, zink ion and me al ion binding. Two dis inc QTL we e iden i ied o a yield, peak a ia ions loca ed a posi ions 9,870,005 and 36,226,978. Va ia ion a bp 9,870,005 is loca ed in he CLEC16A gene (Table 2). Va ia ions o CLEC16A in humans a e associa ed wi h inc eased ype I diabe es isk [55]. In addi ion, milk p o ein pe cen age QTL has p e i- ously been ound om he egion 9.3–10.6 Mb [56]. NEURL1 (BTA26) A e he signi icance es by ixing o he op SNPs, wo QTL we e le on BTA26 o a yield a bp 24,379,571 in NEURL1 gene and a bp 44,802,991 (Table 2). NEURL1 gene is associa ed wi h lac a ion (GO e m 0007595) hus making he a ia ion an in e es ing candi- da e o be a causa i e mu a ion. Conclusions Associa ion analyses among No dic Red Ca le using o e 15 million sequence a ia ions ac oss he whole genome impu ed o o e 4000 p ogeny es ed No dic Red Ca le bulls indica ed se e al a ia ions likely o ha e an impac o milk p oduc ion. We show ha impu a ion is obus and cos -e ec i e way o expand he in o ma ion a ailable and o inc ease knowledge o he causa i e mu a ions a ec ing ai s impo an o p o- duc ion animals. The a ailabili y o he whole genome le el sequence da a opens endless possibili ies o s udy quan i a i e ai a chi ec u e mo e closely. S ill inding he quan i a i e ai nucleo ides is challenging, wi h linkage disequilib ium and many small-e ec QTL c ea ing he puzzle ha is no easy o sol e. Fu he mo e, be e anno- a ion o he ca le genome is equi ed o be able o p edic he e ec s o a ia ions on he pheno ypes mo e accu - a ely. The knowledge om gene in e ac ions (al hough human/ oden based) may help o iden i y likely candida e Iso-Tou u e al. BMC Gene ics (2016) 17:55 Page 9 o 12