Konvergentní genomické znaky domestikace u ovcí a koz | Nature Communications

Odběr vzorků

Domácí ovce (O. aries) a kozy (C. hircus) byly odebrány v Íránu (skupiny IROA a IRCH) a v Maroku (skupiny MOOA a MOCH), celkem 20 zvířat na skupinu (doplňkový obr. 6). Tyto vzorky byly odebrány mezi lednem 2008 a březnem 2012 v severní části Maroka a mezi srpnem 2011 a červencem 2012 v severozápadním Íránu v rámci evropského projektu Nextgen (grantová dohoda č. 244356) v souladu s etickými předpisy směrnice Evropské unie 86/609/EHS. Ušní klipsy byly odebrány z distální části ucha náhodně vybraných zvířat a ihned byly uloženy na jeden den do 96% ethanolu, poté byly přeneseny do silikagelových kuliček až do extrakce DNA.

V severozápadním Íránu byly v rámci domestikační kolébky odebrány vzorky divokých druhů muflonů asijských (O. orientalis) a kozorožců bezoárových (C. aegagrus)21,22 . Třináct tkání muflonů asijských a osmnáct tkání kozorožců bezoárových (resp. skupiny IROO a IRCA, doplňkový obr. 6) bylo odebráno buď ze zvířat chovaných v zajetí, nebo nedávno ulovených, a ze zmrazených vzorků dostupných na íránském ministerstvu životního prostředí. Tento individuální přístup k odběru vzorků má minimalizovat potenciální zkreslení tím, že zamezí nadměrnému zastoupení lokálních vlivů (např, lokálnímu příbuzenskému křížení).

Další údaje

Dále byl sestaven celosvětový panel plemen ovcí a koz (wpOA a wpCH). wpOA zahrnoval 20 celogenomových resekvenovaných vzorků (WGS) s 12násobným pokrytím reprezentujících 20 různých světových plemen, které poskytlo International Sheep Genome Consortium. wpCH se skládal ze 14 WGS vzorků sekvenovaných s 12násobným pokrytím reprezentujících 9 evropských jedinců, tj, 2 francouzské alpské a 2 francouzské saanské vzorky sekvenované INRA, 5 italských saanských vzorků poskytnutých Parco Tecnologico Padano a 5 australských jedinců, tj, 2 vzorky divokých prasat, 2 vzorky rančerů a 1 vzorek kašmíru, které poskytla organizace CSIRO (doplňková data 5).

Produkce dat WGS

Genomická DNA byla úspěšně extrahována ze všech vzorků tkání pomocí soupravy Macherey Nagel NucleoSpin 96 Tissue kit, přičemž byl upraven protokol výrobce. Odběr vzorků tkáně byl proveden v blocích MN se čtvercovými jamkami, aby se získalo 25 mg fragmentů na vzorek. Bylo připraveno tři a půl bloku MN čtvercových 96 a extrakce byla provedena pomocí zařízení Tecan Freedom Evo Liquid handler podle protokolu výrobce. V předlyzačním kroku byly vzorky homogenizovány 180 µl pufru T1 a 25 µl proteinázy K přes noc při 56 °C. Pro úpravu vazebných podmínek bylo přidáno 200 µl pufru BQ1 a destička se vzorky byla inkubována 1 h při 70 °C; následně bylo přidáno 200 µl 100% ethanolu. Lyzáty byly přeneseny na vazebnou desku Nucleospin Tissue a pro odstranění průtoku bylo použito vakuum (-0,2 bar, 5 min). Byly provedeny tři promývací kroky s pufry BW a B5 a opět bylo použito vakuum k odstranění průtoku. Před elucí genomové DNA byla křemenná membrána Nucleospin Tissue binding plate vysušena ve vakuu s tlakem nejméně -0,6 baru po dobu 10 min. Eluční krok byl proveden se 100 µl předehřátého BE pufru (70 °C) a centrifugačním krokem při 3700 otáčkách za minutu po dobu 5 min v 96 jamkách. Genomická DNA byla uchovávána při 4 °C, aby se zabránilo rozmrazování, a byla testována na koncentraci (jako ng/μl) pomocí metody Picogreen a pomocí Nanodropu.

Celé genomy byly resekvenovány z 500 ng genomické DNA, které byly sestříhány na rozsah 150-700 bp pomocí přístroje Covaris® E210 pro každý vzorek a použity pro přípravu knihoven Illumina® podle poloautomatického protokolu. Oprava konců, A-tailing a ligace adaptorů kompatibilních s Illumina® (BioScientific) byly provedeny pomocí systému pro přípravu knihoven SPRIWorks a přístroje SPRI TE (Beckmann Coulter) podle protokolu výrobce. Pro získání většiny fragmentů byla použita selekce velikosti 300-600 bp. Fragmenty DNA byly amplifikovány 12 cykly PCR s použitím sady Platinum Pfx Taq Polymerase Kit (Life® Technologies) a primerů specifických pro adaptér Illumina®. Knihovny byly přečištěny pomocí 0,8x AMPure XP kuliček (Beckmann Coulter) a analyzovány pomocí bioanalyzátoru Agilent 2100 (Agilent® Technologies) a kvantifikace qPCR. Knihovny byly sekvenovány pomocí chemie čtení o délce 100 bází v párové průtokové buňce na přístroji Illumina® HiSeq2000.

Párová čtení Illumina pro Ovis byla mapována na referenční genom ovce (OAR v3.1, sestava GenBank GCA_000298735.146) a pro Capru na referenční genom kozy (CHIR v1.0, sestava GenBank GCA_000317765.147) pomocí BWA-MEM48. Soubor BAM vytvořený pro každého jedince byl roztříděn pomocí Picard SortSam a sekvenčně vylepšen pomocí Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator a GATK IndelRealigner49 a SAMtools calmd50.

Objevení variant bylo provedeno pomocí tří různých algoritmů: Samtools mpileup50, GATK UnifiedGenotyper51 a Freebayes52. Variantní místa byla identifikována nezávisle pro každou ze šesti skupin, přičemž byly použity vícevzorkové režimy vyvolávacích algoritmů: (i) 162 vzorků z MOOA; (ii) 20 vzorků z IROA; (iii) 14 vzorků z IROO; (iv) 162 vzorků z MOCH; (v) 20 vzorků z IRCH; (vi) 19 vzorků z IRCA. U některých skupin byly k dispozici WGS více jedinců v rámci projektu NextGen (viz výše). Vzorky použité v této studii byly vybrány tak, aby byly pokud možno získány vyvážené skupiny 20 jedinců. U skupin IRCA a IROO byly další vzorky k dispozici později a byly přidány pro následné analýzy. Jedinci s nízkou kvalitou zarovnání a volání byli odstraněni, aby se získal konečný soubor dat (Doplňková data 5).

V každé skupině proběhla dvě po sobě jdoucí kola filtrování kvality variantních míst. Fáze filtrování 1 sloučila volání ze tří algoritmů dohromady a zároveň odfiltrovala volání s nejnižší důvěrou. Variantní místo prošlo, pokud bylo vyvoláno alespoň dvěma různými volacími algoritmy s kvalitou variant phred >30. Alternativní alela v místě prošla, pokud byla vyvolána některým z vyvolávacích algoritmů a počet genotypů byl >0. Filtrační fáze 2 použila rekalibraci skóre kvality variant pomocí GATK. Nejprve jsme vytvořili trénovací sadu míst s nejvyšší důvěrou ve variantu v rámci skupiny, kde (i) místo je voláno všemi třemi volacími algoritmy variant s kvalitou phred-scaled variant >100, (ii) místo je bialelické, (iii) počet minoritních alel je alespoň 3, přičemž se počítají pouze vzorky s genotypem s kvalitou phred-scaled >30. Trénovací soubor byl použit k sestavení Gaussova modelu pomocí nástroje GATK VariantRecalibrator s použitím následujících anotací variant z nástroje UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Na celý soubor dat byl aplikován Gaussův model, který vygeneroval VQSLOD (logaritmický poměr šancí, že se jedná o pravou variantu). Místa byla filtrována, pokud VQSLOD <mezní hodnota. Mezní hodnota byla pro každou skupinu stanovena následujícím způsobem: Minimální hodnota VQSLOD = {medián hodnoty VQSLOD pro varianty tréninkového souboru}-3 × {medián absolutní odchylky VQSLOD variant tréninkového souboru}. Poměr přechodných/transverzních SNP naznačil, že zvolené kritérium mezní hodnoty poskytuje nejlepší rovnováhu mezi selektivitou a citlivostí.

Byly vytvořeny soubory volání SNP pro šest skupin zvířat Ovis a Capra (tj. íránské a marocké domestikanty a divoká zvířata pro každý rod). Protože analýzy prováděné v této studii vyžadovaly meziskupinové srovnání, vytvořili jsme soubory genotypových volání na konzistentní sadě míst SNP pro všechna zvířata z kterékoli skupiny. Pro každý rod jsme sloučili místa volání variant z jeho tří skupin a ponechali jsme pouze bialelické pozice bez chybějících údajů. Genotypy byly znovu vyvolány na každém bialelickém místě SNP pro všechny zájmové jedince pomocí nástroje GATK UnifiedGenotyper s použitím volby GENOTYPE_GIVEN_ALLELES. V této fázi byl seznam jedinců rozšířen o zvířata patřící do panelů světových plemen ovcí a koz (wpOA a wpCH) a další divoké vzorky, které byly v této fázi k dispozici (4 O. orientalis a 4 C. aegagrus). Genotypy byly vylepšeny a rozfázovány v rámci skupin Beagle 453 a poté odfiltrovány tam, kde byla pravděpodobnost genotypu menší než 0,95. Nakonec jsme odfiltrovali místa, která byla monomorfní v různých podskupinách jedinců použitých v této studii (viz níže).

Pro porovnání signálů selekce zjištěných mezi Ovis a Capra jsme provedli křížové zarovnání mezi oběma referenčními genomy. Nejprve jsme použili pipeline párového zarovnání z kódové základny Ensembl verze 6954 k zarovnání referenčních genomů ovce (OARv3.1) a kozy (CHIR1.0). Tato pipeline používá LastZ55 k zarovnání na úrovni DNA, po kterém následuje následné zpracování, při kterém jsou zarovnané bloky zřetězeny podle jejich umístění v obou genomech. Pipeline párového zarovnávání LastZ je rutinně spouštěna systémem Ensembl pro všechny podporované druhy, ale koza zatím není do systému Ensembl zahrnuta. Abychom se vyhnuli zkreslení ve prospěch jednoho z druhů, vytvořili jsme dvě různá mezidruhová zarovnání. Jedno použilo ovci jako referenční genom a kozu jako nereferenční, zatímco druhé použilo kozu jako referenční genom a ovci jako nereferenční. Rozdíl spočívá v tom, že genomové oblasti referenčního druhu jsou nuceny mapovat jednoznačně k jednotlivým lokusům nereferenčního druhu, zatímco nereferenční genomové oblasti mohou být mapovány k více místům referenčního druhu. Pro segmenty chromozomů jednoho referenčního genomu jsme získali souřadnice na nereferenčním genomu. Nakonec jsme pro SNP objevené u jednoho rodu použili zarovnání celého genomu s referenčním genomem druhého rodu k určení odpovídajících pozic (Doplňková tabulka 6).

Genetická struktura

Pro popis genetické diverzity v rámci skupin jsme použili VCFtools56 k výpočtu souhrnné statistiky genetické variability na 73 jedincích pro Ovis (tj, 13 IROO, 20 IROA, 20 MOOA a 20 wpOA) a 72 jedinců pro Capra (tj. 18 IRCA, 20 IRCH, 20 MOCH a 14 wpCH). Měřenými statistikami byly celkový počet polymorfních variant (S) pro celý soubor jedinců v každém rodu a v každé skupině, průměrná nukleotidová diverzita (π) v každé skupině a koeficient inbreedingu (F) pro každého jedince. V rámci každého rodu byly rozdíly mezi divokou skupinou a každou domácí skupinou testovány pomocí jednostranného t-testu pro individuální hodnoty inbreedingu a genetické zátěže a oboustranného Mannova-Whitneyho testu pro nukleotidovou diverzitu na lokalitu.

Celková divergence mezi čtyřmi skupinami v rámci každého rodu (tj, divokými, íránskými a marockými domestikanty a světovým panelem) byla odhadnuta pomocí všech bialelických SNP a průměrné vážené párové Fst podle Weira a Cockerhama57 implementované v nástroji VCFtools56. Genetická struktura mezi skupinami byla posouzena pomocí metody shlukování sNMF26 po ořezání datového souboru za účelem odstranění SNP s vazebnou nerovnováhou (r²) větší než 0,2 pomocí VCFtools. Vazbová nerovnováha (r²) byla vypočtena mezi páry SNP v rámci klouzavých oken po 50 SNP, přičemž jeden SNP na pár byl náhodně odstraněn, když r² bylo větší než 0,2. Pro každou analýzu sNMF bylo provedeno pět běhů se stejným počtem shluků (K) s hodnotami K od 1 do 10. K určení nejpravděpodobnějšího řešení shlukování jsme použili kritérium křížové entropie, nicméně byla zkoumána také alternativní rozdělení pro různé počty K, abychom posoudili, jak byli jedinci rozděleni mezi shluky.

Pro oddělení sdíleného původu a příměsi jsme spustili program TreeMix27 ke společnému odhadu rozdělení populace a následných příměsí pomocí ořezaného souboru dat použitého pro sNMF. Spustili jsme TreeMix s volbou -global, abychom zpřesnili naše závěry o maximální pravděpodobnosti. Strom TreeMix jsme zakořenili pomocí rozdělení na divoké a domácí jedince. Velikost bloku pro jackknifing byla -k 500 SNP, což přibližně odpovídá 150 kb, což přesahuje průměrné bloky LD zjištěné u ovcí i koz. Vytvořili jsme strom maximální věrohodnosti bez migrace a poté jsme přidali migrační události a zkoumali jsme přírůstkovou změnu rozptylu vysvětleného modelem a hodnoty reziduí mezi jedinci. Cílem bylo odhalit případné vysoké reziduální hodnoty nebo migrační hrany mezi divokými a domácími jedinci. Abychom dále prozkoumali statistickou relevanci možných vektorů příměsi identifikovaných programem TreeMix (doplňková tabulka 3), vypočítali jsme třípopulační test f328 jako formální test genetické introgrese pomocí programu qp3Pop ze sady ADMIXTOOLS58 pro každou kombinaci skupin. U plemene Capra byla skupina wpCH rozdělena mezi australská plemena, francouzská plemena a italská plemena. Výsledky jsou uvedeny v doplňkových údajích 2.

Demografická inference

Pro každý rod jsme provedli analýzy demografické inference předků pomocí modelu MSMC implementovaného v softwaru MSMC225. MSMC je založen na párově sekvenční markovské koalescenci59; jako vstup však používá haplotypy fázových sekvenčních dat genomu. Pro každou analýzu jsme použili dva jedince z jedné skupiny, tedy 4 haplotypy. Každá analýza byla opakována pro další náhodný soubor dvou jedinců, tj. opakování analýzy pro každou skupinu. Vstupní a výstupní soubory byly generovány a analyzovány pomocí skriptů python dodaných se softwarem MSMC a naleznete je na adrese https://github.com/stschiff/msmc-tools. Parametry analýz byly ponechány ve výchozím nastavení s výjimkou mutační rychlosti, která byla nastavena na 2,5×10-8, a délky generace na 2 roky. Abychom odhadli nejistotu časových odhadů, měnili jsme tyto parametry (mutační rychlost 2,5×10-8 a 1,0×10-8 v kombinaci s délkou generace 2 a 4 roky) a získali hrubý odhad doby domestikace (viz doplňkový obr. 2).

Genetická zátěž

Genetická zátěž byla odhadována dvěma způsoby. Zaprvé výpočtem genetické zátěže pro každého jedince jako součtu škodlivých fitness efektů na všech genomických pozicích kódujících proteiny podle metody Librado et al.60. Dále byl proveden výpočet genetické zátěže pro každého jedince jako součtu škodlivých fitness efektů na všech genomických pozicích kódujících proteiny podle metody Librado et al. Stručně řečeno, jako zástupný ukazatel evolučního omezení jsme použili skóre PhyloP ze zarovnání 46 savců (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Z tohoto zarovnání jsme identifikovali místa kódující proteiny vyvíjející se pod funkčním omezením (skóre phyloP ≥1,5). U každého genomu Ovis nebo Capra jsme pak zjišťovali, zda tato místa byla mutována. Pokud ano, sečetli jsme skóre phyloP za všechna mutovaná místa, takže mutace ve vysoce omezených místech přispívají úměrně více k odhadu celkové zátěže. Tím jsme získali odhad zátěže pro každý genom ovce/kozy. Nakonec jsme pro získání průměrného zatížení na jedno místo vydělili celkovým počtem analyzovaných pozic. Stojí za zmínku, že jsme se podmínili homozygotními místy, abychom se vyhnuli modelování koeficientu dominance mutací na heterozygotních místech (např. recesivní, intermediární, dominantní). Zadruhé jsme gen po genu porovnali genetickou škodlivou zátěž u divokých a domestikovaných skupin Ovis provedením Wilcoxonova testu, přičemž alternativní hypotézou bylo, že domácí zvířata mají větší zátěž než divocí příbuzní. p-hodnoty byly korigovány na vícenásobné testování61 a použili jsme prahovou hodnotu upravené p-hodnoty < 0,05. Na souboru genů vykazujících významné zvýšení genetické zátěže jsme provedli analýzu obohacení genové ontologie pomocí programu WebGestalt62,63 . Vzhledem k tomu, že referenční genomy jsou u genů nedostatečně anotovány, vycházeli jsme z ortologů s jednou kopií mezi naším druhem a člověkem a myší. Geny z chromozomu X byly z podkladového souboru vyloučeny. Tuto analýzu jsme neprovedli u Capry kvůli vyššímu příbuzenskému křížení pozorovanému u divokých vzorků.

Detekce selekčních znaků

Pro detekci znaků selekce související s domestikací jsme použili všechny bialelické SNP vykazující menší frekvenci alel než 0,10 alespoň v jedné ze tří testovaných skupin (tj. íránské a marocké domestikované skupiny a divoké skupiny pro každý rod). Protože jsme očekávali, že znaky selekce související s procesem domestikace budou přítomny u všech domácích zvířat, přijali jsme následující obecnou strategii: testovali jsme pomocí hapFLK29 (viz doplňková poznámka 5 a doplňkové obrázky 9, 10 a 11) pro každý rod divokou skupinu proti každé z tradičně obhospodařovaných domácích skupin (tj. íránské a marocké) a zaměřili jsme se na ty společné oblasti domněle podléhající selekci, které byly zjištěny v obou případech. Velikosti vzorků skupin (n = 13-20) byly v souladu s požadavky metody29. Vizuálně jsme zkontrolovali, zda konzistentní znaky selekce zjištěné pomocí hapFLK byly přítomny také v odpovídajícím souboru světových panelů každého rodu, ale tyto skupiny jsme do statistického testu nezahrnuli vzhledem k jejich vícekmenovému složení. Nakonec jsme hledali společné signály selekce mezi druhy Ovis a Capra pomocí stratifikovaného přístupu FDR. Strategie je znázorněna na doplňkovém obr. 4.

Provedli jsme hapFLK testy pro kontrast divoké skupiny s každou z íránských a marockých skupin v každém rodu. Matice příbuznosti byla vypočtena z Reynoldových genetických vzdáleností64 mezi dvojicemi skupin s použitím náhodné podmnožiny jednoho procenta variant. Odvozený populační strom byl sestaven pomocí algoritmu spojování sousedů. Pro každý SNP jsme provedli test hapFLK, který zahrnuje haplotypovou informaci, aby se zvýšila síla pro detekci selekčních zásahů. Pro každý testovaný SNP statistika hapFLK vypočítala odchylku haplotypových frekvencí vzhledem k neutrálnímu modelu odhadnutému maticí příbuznosti65. K využití informací o vazebné nerovnováze používá hapFLK Scheetův a Stephensův66 vícebodový model pro multilokusové genotypy, který lze přizpůsobit nefázovaným datům. Jednou z hlavních aplikací tohoto modelu je provádění fázového odhadu (software fastPHASE66). V naší analýze byl model vycvičen na nefázovaných datech, a proto naše analýza zohledňuje nejistotu fáze. Metoda byla použita k přeskupení lokálních haplotypů podél chromozomů do určeného počtu shluků K nastaveného na 25 pomocí skrytého Markovova modelu.

Pro určení společných oblastí údajně podléhajících selekci ve dvou tradičně spravovaných domácích skupinách pro každý rod jsme zkombinovali dvě předchozí analýzy hapFLK. Pro každou analýzu bylo hapFLK skóre přizpůsobeno χ2 rozdělení pro získání p-hodnot (skript je k dispozici na https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Výsledky obou kontrastů mezi divokou skupinou a každou z domácích skupin byly zkombinovány pomocí Stoufferovy metody67 , aby se získaly jediné p-hodnoty pro srovnání divokých a domácích zvířat. Nakonec byl na celý soubor SNP aplikován rámec FDR68 , aby se kombinované p-hodnoty převedly na q-hodnoty. SNP vykazující q-hodnoty < 10-2 byly ponechány a seskupeny do genomových oblastí, pokud byly od sebe vzdáleny méně než 50 kb.

Pro zjištění, zda byl signál selekce sdílen mezi Ovis a Capra, jsme nejprve použili křížové zarovnání obou referenčních genomů k identifikaci homologických segmentů. Poté jsme použili stratifikovaný rámec FDR69. Tento přístup vychází ze skutečnosti, že vzhledem k předběžným informacím v genetických datech69 existuje přirozená stratifikace testů, protože základní rozdělení pravdivých alternativních hypotéz se může lišit podle různé dynamiky různých genomických oblastí, což vede k různým rozdělením p-hodnot. To vyžaduje získat p-hodnoty upravené o FDR (tj. q-hodnoty) odděleně pro různé vrstvy. Konvergence jsme hledali v každém rodu tak, že jsme oddělili oblasti homologické s oblastmi zjištěnými u druhého rodu (označované jako sdílená vrstva) a zbytek genomu (označované jako obecná vrstva). Extrahovali jsme p-hodnoty zvlášť pro každou z obou definovaných vrstev a poté jsme vypočítali q-hodnoty pomocí rámce FDR. Tyto stratifikované q-hodnoty byly konečnými veličinami uvažovanými pro statistickou významnost (<10-2) pro detekci SNP podléhajících selekci a jejich sloučení do příslušných genomových oblastí.

Pro testování konvergentních znaků selekce odlišujících divoká a domácí zvířata v obou rodech jsme zkoumali vztah mezi prahem významnosti použitým pro q-hodnoty (které jsme nechali kolísat od 0,2 do 0,002) v jednom rodě a odhadovanou pravděpodobností, že SNP je selektován ve sdílené vrstvě druhého rodu, pomocí Storeyho a kol.70 přístup. Zvýšení odvozené pravděpodobnosti se snížením prahové hodnoty aplikované na hodnotu q (zvýšení přísnosti) naznačuje, že čím významnější je oblast v jednom rodu, tím pravděpodobněji bychom našli významné SNP v druhém rodu.

Odfiltrovali jsme selekční signály, které nebyly konzistentní mezi třemi domácími skupinami. Pro každou zjištěnou oblast jsme použili fázové haplotypy jedinců, které jsme shlukovali pomocí stromů Neighbor-Joining na základě procenta identity mezi sekvencemi. Ponechány byly pouze oblasti vykazující konzistentní signály (doplňkový obr. 5).

Abychom mohli odvodit, zda signály selekce zjištěné pomocí hapFLK naznačují uvolnění selekce nebo pozitivní selekci u domestiků, odhadli jsme rozdíl v nukleotidové diverzitě (π) na každé domnělé oblasti podléhající selekci mezi divokými a domácími skupinami. Tento rozdíl jsme vyjádřili jako index Δπ, který jsme pro každou genomovou oblast vypočítali jako rozdíl mezi π vypočteným pro divokou skupinu a průměrem π pro íránskou a marockou domácí skupinu minus rozdíl v π mezi těmito dvěma skupinami vypočtený pro celý genom:

$$\Delta \pi = \left( {\pi _{{\rm wilds}} – \pi _{{\mathrm{\operatorname{{iran-morocco}}}}}} \right)_{{\mathrm{\operatorname{genomická-oblast}}}} – \left( {\pi _{{\rm wilds}} – \pi _{{\mathrm{\operatorname{iran-morocco}}}}} \right)_{{\mathrm{\operatorname{celý-genom}}}}$$

Záporná hodnota by znamenala, že nukleotidová diverzita je nižší ve skupině wild ve srovnání s průměrem obou skupin domestiků, a byla by považována za projev uvolnění selekce v těchto posledních skupinách, diverzifikační selekce u domestiků nebo pozitivní selekce ve wilds. Naopak kladná hodnota by ukazovala na směrovou pozitivní nebo stabilizující selekci, která proběhla v domestikovaných skupinách. Pomocí haplotypového shlukování jsme také v každé oblasti ručně ověřili, zda zjištěný selekční zásah potvrdil indicie dané indexem Δπ.

Funkční interpretace jsme provedli následovně. Pro každou oblast podléhající selekci jsme vzali v úvahu oblast plus 50 kb na každé straně pro určení funkčních rolí a 5 kb před a za geny a vyhodnotili jsme překrytí těchto souřadnic pro zachování zájmových genů. Nakonec jsme měli za to, že gen souvisí s danou zjištěnou oblastí, pokud se pozice oblasti a genu překrývaly. Poté jsme posoudili, který gen byl s největší pravděpodobností cílem selekce, a to tak, že jsme vzali v úvahu nejbližší gen k hornímu signálu, tj. pozici s nejnižší hodnotou q v rámci regionu. Geny byly funkčně anotovány pomocí Uniprot (http://www.uniprot.org/) tak, že se zvážilo jejich zapojení do 30 podřízených termínů (tj. přímých sestupů termínů) kategorie „Biologický proces“ (tj. GO:0008150). Vyhledali jsme všechny termíny GO odpovídající každému genu (doplňková data 4) pro 30 z 33 kategorií, protože jsme nebrali v úvahu tři termíny, které nebyly zapojeny do funkcí savců (tj. využití síry GO:0006791, využití fosforu GO:0006794, využití uhlíku GO:0015976). Provedli jsme dva χ2-testy, abychom porovnali rozložení genů v kategoriích GO, tj. i) geny pod selekcí z rodově specifických oblastí oproti genům z homologických oblastí a ii) všechny geny pod selekcí oproti 18 689 lidským genům přiřazeným k termínům GO ve Swiss-Prot. Abychom mohli interpretovat funkce genů v kontextu hospodářských zvířat, vyhledali jsme také dostupné informace z literatury o jejich fenotypových účincích.

Nakonec jsme k nalezení SNP v rámci dříve zjištěných oblastí, které se nejvíce lišily mezi divokými a domácími skupinami, použili statistiku FLK. Stejně jako hapFLK představuje odchylku alelických frekvencí jednotlivých markerů vzhledem k neutrálnímu modelu odhadnutému pomocí matice příbuznosti65. Pro dosazení výsledků z obou analýz do χ2 rozdělení a kombinaci získaných p-hodnot byl použit stejný postup, jaký byl použit pro test hapFLK. Nerovnoměrné rozdělení p-hodnot však vylučovalo použití rámce FDR a vybrali jsme SNP v rámci oblastí zjištěných pomocí hapFLK, které vykazovaly p-hodnoty <10-4. Pro tyto SNP jsme použili anotace VEP (Variant Effect Predictor)71 , které byly vygenerovány z anotace genomu ovcí OARv3.1 v databázi Ensembl v74 pro Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) a z anotace genomu koz CHIR1.0 vytvořené pipeline pro anotaci eukaryotických genomů NCBI pro Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). SNP byly klasifikovány jako intergenické, upstream a downstream (včetně UTR) a intronické a exonické pozice. Rozdíly mezi distribucí SNP s p-hodnotami FLK <10-4 a všemi SNP použitými pro detekci selekčních znaků byly zkoumány pomocí χ2-testu.

Dostupnost dat

Sekvence a metadata generovaná pro 73 vzorků Ovis a 72 vzorků Capra použitých v těchto analýzách jsou veřejně dostupné. Obecné informace a všechny soubory vcf lze nalézt na webových stránkách Ensembl (http://projects.ensembl.org/nextgen/). Všechny soubory Fastq, soubory Bam a sestavy de novo pro O. orientalis a C. aegagrus lze nalézt na stránkách European Nucleotide Archive (https://www.ebi.ac.uk/ena) pod přístupovým kódem projektu Nextgen (PRJEB7436).