Sampling

Domestic sheep (O. aries) and kecskék (C. hircus) were sampling in Iran (IROA and IRCH groups, respectively) and Morocco (MOOA and MOCH groups, respectively) for a total of 20 animals per group (Supplementary Fig. 6). Ezeket a mintákat 2008 januárja és 2012 márciusa között gyűjtötték Marokkó északi részén, illetve 2011 augusztusa és 2012 júliusa között Irán északnyugati részén, a Nextgen európai projekt (244356 számú támogatási megállapodás) keretében, az Európai Unió 86/609/EGK irányelvének etikai előírásaival összhangban. A fülcsipeszeket véletlenszerűen kiválasztott állatok fülének disztális részéből vettük, és azonnal 96%-os etanolban tároltuk egy napig, majd szilikagél gyöngyökbe helyeztük át a DNS kivonásig.

A vadon élő ázsiai muflon (O. orientalis) és a bezoari kőszálibika (C. aegagrus) fajokat Északnyugat-Iránban, a háziasítás bölcsőjében21,22 gyűjtöttük. Tizenhárom ázsiai muflon és 18 bezoari kőszálibika szövetét (IROO-, illetve IRCA-csoport, 6. kiegészítő ábra) fogságban tartott vagy nemrégiben elejtett állatokból, illetve az iráni környezetvédelmi minisztériumban rendelkezésre álló fagyasztott mintákból gyűjtötték. Ez az egyedalapú mintavételi megközelítés a helyi hatások felülreprezentálásának elkerülése révén a lehetséges torzítások minimalizálására szolgál (pl,

Kiegészítő adatok

Kiegészítésképpen összeállítottunk egy világméretű fajtapanelt a juhok és kecskék számára (wpOA és wpCH). A wpOA 20 teljes genom-újraszekvenálási (WGS) mintát tartalmazott 12x-es lefedettséggel, amelyek 20 különböző világméretű fajtát képviselnek, és amelyeket a Nemzetközi Juh Genom Konzorcium biztosított. wpCH 14 WGS-mintát tartalmazott 12x-es lefedettséggel, amelyek 9 európai egyedet képviselnek, 2 francia alpesi és 2 francia saanen mintát, amelyeket az INRA szekvenált, 5 olasz saanen mintát, amelyeket a Parco Tecnologico Padano biztosított, és 5 ausztrál egyedet, azaz, 2 Boer, 2 Rangeland és 1 kasmír minta, amelyeket a CSIRO biztosított (5. kiegészítő adat).

A WGS-adatok előállítása

A genomi DNS-t sikeresen kivontuk valamennyi szövetmintából a Macherey Nagel NucleoSpin 96 Tissue kit segítségével, a gyártó protokollját adaptálva. A szövetmintavételt MN négyszögletes lyukblokkokban végeztük, hogy mintánként 25 mg fragmentumot nyerjünk. Három és fél MN négyzet-96-os blokkot készítettünk, és az extrakciót Tecan Freedom Evo folyadékkezelővel végeztük a gyártó protokollját követve. Egy előlizálási lépést végeztünk a minták homogenizálására 180 µl T1 pufferrel és 25 µl proteináz K-val egy éjszakán át 56 °C-on. A kötési feltételek beállításához 200 µl BQ1 puffert adtunk hozzá, és a mintalemezt 1 órán át 70 °C-on inkubáltuk; ezt követően 200 µl 100%-os etanolt adtunk hozzá. A lizátumokat átvittük a Nucleospin Tissue kötőlemezre, és vákuumot (-0,2 bar, 5 perc) alkalmaztunk az átfolyás eltávolítására. Három mosási lépést végeztünk BW és B5 pufferrel, majd ismét vákuumot alkalmaztunk az átfolyás eltávolítására. A genomi DNS elúciója előtt a Nucleospin Tissue kötőlemez szilikamembránt 10 percig legalább -0,6 bar nyomáson vákuumban szárítottuk. Az elúciós lépést 100 µl előmelegített BE pufferrel (70 °C) és 5 percig 3700 rcf-en végzett centrifugálási lépéssel végeztük a 96-PCR lyukakban. A genomiális DNS-t 4 °C-on tároltuk a fagyasztás-felolvasztás elkerülése érdekében, és a koncentrációt (ng/μl-ben kifejezve) a Picogreen módszerrel és Nanodrop segítségével vizsgáltuk.

A teljes genomokat 500 ng genomiális DNS-ből reszekvenáltuk, amelyet minden egyes minta esetében a Covaris® E210 műszerrel 150-700 bp tartományba nyírtunk, és félautomatizált protokollal Illumina® könyvtárkészítésre használtunk. A végjavítást, az A-tailinget és az Illumina®-kompatibilis adaptorok (BioScientific) ligálását az SPRIWorks könyvtárkészítő rendszer és az SPRI TE műszer (Beckmann Coulter) segítségével végeztük a gyártó protokollja szerint. A legtöbb fragmentum visszanyerése érdekében 300-600 bp méretű szelekciót alkalmaztunk. A DNS-fragmentumokat 12 ciklusos PCR-rel, Platinum Pfx Taq Polymerase Kit (Life® Technologies) és Illumina® adapter-specifikus primerek használatával amplifikáltuk. A könyvtárakat 0,8x AMPure XP gyöngyökkel (Beckmann Coulter) tisztítottuk, majd Agilent 2100 Bioanalyzerrel (Agilent® Technologies) és qPCR kvantifikációval elemeztük. A könyvtárak szekvenálása 100 bázishosszúságú olvasáskémia alkalmazásával történt párosított végű áramlási cellában az Illumina® HiSeq2000-en.

Az Illumina párosított végű olvasatait az Ovis esetében a juh referencia genomra (OAR v3.1, GenBank assembly GCA_000298735.146), a Capra esetében pedig a kecske referencia genomra (CHIR v1.0, GenBank assembly GCA_000317765.147) térképeztük le a BWA-MEM48 segítségével. Az egyes egyedekre előállított BAM-fájlt a Picard SortSam segítségével rendeztük, és szekvenciálisan Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator és GATK IndelRealigner49, valamint SAMtools calmd50 segítségével javítottuk.

A változatok felfedezését három különböző algoritmus segítségével végeztük: Samtools mpileup50, GATK UnifiedGenotyper51 és Freebayes52. A variánshelyek azonosítása mind a hat csoport esetében egymástól függetlenül történt, a hívó algoritmusok többmintás módját használva: (i) 162 minta a MOOA-tól; ii) 20 minta az IROA-tól; iii) 14 minta az IROO-tól; iv) 162 minta a MOCH-tól; v) 20 minta az IRCH-tól; vi) 19 minta az IRCA-tól. Néhány csoport esetében a NextGen projekt részeként több egyed WGS-e állt rendelkezésre (lásd fent). A jelen tanulmányban felhasznált mintákat úgy választottuk ki, hogy lehetőség szerint kiegyensúlyozott, 20 egyedből álló csoportokat kapjunk. Az IRCA és az IROO csoportok esetében további minták váltak elérhetővé egy későbbi szakaszban, amelyeket a későbbi elemzésekhez adtak hozzá. Az alacsony illesztési és hívásminőségű állatokat eltávolítottuk, hogy megkapjuk a végleges adathalmazt (5. kiegészítő adat).

Minden csoporton belül két egymást követő fordulóban végeztük el a variánshelyek minőségének szűrését. Az 1. szűrési szakasz a három algoritmusból származó hívásokat egyesítette, miközben kiszűrte a legalacsonyabb konfidenciájú hívásokat. Egy variánshely akkor ment át, ha legalább két különböző hívó algoritmus hívta meg >30 phred variánsminőséggel. A hely alternatív allélja akkor ment át, ha a hívó algoritmusok bármelyike hívta, és a genotípusszám >0 volt. A 2. szűrési szakasz a GATK által végzett Variant Quality Score Recalibration-t használta. Először is létrehoztuk a csoporton belül a legmagasabb konfidenciájú variáns helyek gyakorlóhalmazát, ahol (i) a helyet mindhárom variánshívó algoritmus >100 fred-skálázott variánsminőséggel hívja, (ii) a hely biallelikus, (iii) a minor allélszám legalább 3, miközben csak a >30 fred-skálázott genotípusminőségű mintákat számoltuk. A képzési halmazból Gauss-modellt építettünk a GATK VariantRecalibrator eszközzel a következő, UnifiedGenotyperből származó variáns annotációk felhasználásával: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. A teljes adathalmazra egy Gauss-modellt alkalmaztunk, amely egy VQSLOD (valódi variánsnak lenni log odds ratio) értéket generált. A helyszíneket kiszűrtük, ha a VQSLOD <vágási érték. A cutoff-értéket az egyes csoportok esetében a következők szerint határoztuk meg: VQSLOD minimum = {az edzéskészlet változatainak VQSLOD medián értéke}-3 × {az edzéskészlet változatainak VQSLOD abszolút eltérésének mediánja}. Az átmenet/transzverzió SNP-arány azt sugallta, hogy a választott cutoff-kritérium adta a legjobb egyensúlyt a szelektivitás és az érzékenység között.

SNPs híváskészleteket generáltunk az Ovis és Capra állatok hat csoportjára (azaz az iráni és marokkói háziállatokra, valamint a vadakra minden nemzetség esetében). Mivel az ebben a tanulmányban végzett elemzésekhez csoportközi összehasonlításokra volt szükség, minden csoportból minden állat esetében konzisztens SNP-helyekre vonatkozó genotípus-híváskészleteket hoztunk létre. Minden egyes nemzetség esetében egyesítettük a három csoportból származó variáns hívóhelyeket, és csak a hiányzó adatok nélküli biallelikus pozíciókat tartottuk meg. A genotípusokat minden egyes biallelikus SNP-helyen újrahívtuk a GATK UnifiedGenotyperrel az összes érdekes egyedre vonatkozóan, a GENOTYPE_GIVEN_ALLELES opciót használva. Ebben a szakaszban az egyedek listáját kibővítettük a juh- és kecskefajták világtenyésztési paneljeihez tartozó állatokkal (wpOA és wpCH), valamint az ebben a szakaszban elérhetővé vált további vadon élő mintákkal (4 O. orientalis és 4 C. aegagrus). A genotípusokat a Beagle 453 segítségével javítottuk és csoporton belül fázisosítottuk, majd kiszűrtük azokat, ahol a genotípus valószínűsége 0,95-nél kisebb volt. Végül kiszűrtük azokat a helyeket, amelyek monomorfikusak voltak a vizsgálatban használt egyedek különböző alcsoportjaiban (lásd alább).

Az Ovis és Capra között észlelt szelekciós jelek összehasonlítása érdekében keresztalignálást végeztünk a két referencia genom között. Először az Ensembl 69-es kiadású kódbázis54 páronkénti igazítási csővezetékét használtuk a juh (OARv3.1) és a kecske (CHIR1.0) referencia genomjainak igazításához. Ez a csővezeték a LastZ55 -t használja a DNS-szintű igazításhoz, amelyet utófeldolgozás követ, amelyben az igazított blokkokat a két genomban elfoglalt helyük szerint láncoljuk össze. A LastZ páros igazítási csővezetéket az Ensembl rutinszerűen futtatja minden támogatott faj esetében, de a kecske még nem szerepel az Ensemblben. Annak érdekében, hogy elkerüljük a torzítást bármelyik faj irányába, két különböző fajközi összehangolást készítettünk. Az egyikben a juhot használtuk referencia genomként és a kecskét nem referenciaként, míg a másikban a kecskét használtuk referencia genomként és a juhot nem referenciaként. A különbség abban áll, hogy a referenciafaj genomi régióit arra kényszerítettük, hogy egyértelműen a nem referenciafaj egyetlen lókuszához kapcsolódjanak, míg a nem referenciafaj genomi régiói a referenciafaj több lókuszához is kapcsolódhatnak. Az egyik referencia genom kromoszómáinak szegmenseire megkaptuk a koordinátákat a nem referencia genomon. Végül az egyik nemzetségben felfedezett SNP-k esetében a teljes genomnak a másik nemzetség referencia genomjával való összehangolását használtuk a megfelelő pozíciók azonosítására (6. kiegészítő táblázat).

Genetikai struktúra

A csoportokon belüli genetikai diverzitás leírásához a VCFtools56 segítségével kiszámítottuk a genetikai variáció összefoglaló statisztikáit az Ovis 73 egyedére (pl, 13 IROO, 20 IROA, 20 MOOA és 20 wpOA) és 72 egyedet a Capra esetében (azaz 18 IRCA, 20 IRCH, 20 MOCH és 14 wpCH). A mért statisztikák a következők voltak: a polimorf variánsok teljes száma (S) a teljes egyedkészletre vonatkozóan minden nemzetségben és minden csoporton belül, az átlagos nukleotiddiverzitás (π) minden csoporton belül és a beltenyésztési együttható (F) minden egyedre vonatkozóan. Az egyes nemzetségeken belül a vadon élő csoport és az egyes hazai csoportok közötti különbségeket egyoldalú t-próbával teszteltük az egyéni beltenyésztési és genetikai terhelés értékek esetében, valamint kétoldalú Mann-Whitney-teszttel a nukleotiddiverzitás helyenkénti értékére.

A négy csoport közötti teljes divergencia az egyes nemzetségeken belül (azaz, vad, iráni és marokkói háziállatok, valamint a világpanel) közötti különbséget az összes biallelikus SNP és a Weir és Cockerham57 szerinti átlagos súlyozott páronkénti Fst segítségével becsültük meg, a VCFtools56-ban implementált módon. A csoportok közötti genetikai struktúrát az sNMF26 klaszterezési módszerrel értékeltük, miután az adathalmazt a VCFtools segítségével a 0,2-nél nagyobb kapcsolati egyenlőtlenséggel (r²) rendelkező SNP-k eltávolítása céljából megkurtítottuk. A kapcsolódási egyenlőtlenséget (r²) az SNP-párok között 50 SNP-ből álló csúszóablakon belül számoltuk ki, és páronként egy SNP-t véletlenszerűen eltávolítottunk, ha az r² nagyobb volt 0,2-nél. Minden sNMF-elemzéshez öt futtatást végeztünk ugyanannyi klaszterrel (K), a K értéke 1 és 10 között volt. A kereszt-entrópia kritériumot használtuk a legvalószínűbb klaszterezési megoldás azonosítására, azonban alternatív partíciókat is megvizsgáltunk különböző K számok esetén, hogy felmérjük, hogyan oszlanak meg az egyének a klaszterek között.

A közös származás és a keveredés szétválasztása érdekében lefuttattuk a TreeMix27 programot a populációs felosztások és a későbbi keveredési események együttes becslésére az sNMF-hez használt metszett adathalmaz felhasználásával. A TreeMixet a -global opcióval futtattuk, hogy finomítsuk a maximális valószínűségű következtetéseinket. A TreeMix-fát a vadon élő és a házi egyedek közötti felosztással gyökereztettük. A jackknifinghez használt blokkméret -k 500 SNP volt, ami körülbelül 150 kb-nak felel meg, ami meghaladja a juhoknál és a kecskéknél talált átlagos LD blokkokat. Létrehoztunk egy Maximum Likelihood fát migráció nélkül, majd migrációs eseményeket adtunk hozzá, és megvizsgáltuk a modell által magyarázott variancia és az egyedek közötti reziduális értékek növekményes változását. A cél az volt, hogy a vadon élő és a háziasított egyedek közötti esetleges magas reziduális értékeket vagy migrációs éleket felderítsük. A TreeMix által azonosított lehetséges admixtúra-vektorok statisztikai relevanciájának további feltárása érdekében (3. kiegészítő táblázat) a genetikai introgresszió formális tesztjeként kiszámítottuk az f328 hárompopulációs tesztet az ADMIXTOOLS csomag qp3Pop programjának58 segítségével minden egyes csoportkombinációra. A Capra esetében a wpCH csoportot az ausztrál fajták, a francia fajták és az olasz fajták között osztottuk fel. Az eredményeket a 2. kiegészítő adatokban közöljük.

Demográfiai következtetés

Minden nemzetség esetében elvégeztük az ősi demográfiai következtetés elemzését az MSMC2 szoftverben25 implementált MSMC modell segítségével. Az MSMC a páronkénti szekvenciális Markov-koaleszcencián59 alapul, azonban bemenetként a fázisos genomszekvencia-adatok haplotípusait használja. Minden elemzéshez egy csoportból két egyedet, tehát 4 haplotípust használtunk. Minden elemzést megismételtünk egy másik véletlenszerű, két egyedből álló csoportra, azaz csoportonként egy ismétléssel. A bemeneti és kimeneti fájlokat az MSMC szoftverhez mellékelt python szkriptekkel generáltuk és elemeztük, amelyek a https://github.com/stschiff/msmc-tools címen találhatók. Az elemzési paraméterek alapértelmezettek maradtak, kivéve a mutációs rátát, amelyet 2,5×10-8-ra, a generáció hosszát pedig 2 évre állítottuk be. Az időbecslések bizonytalanságának becslése érdekében ezeket a paramétereket variáltuk (2,5×10-8 és 1,0×10-8 mutációs ráta a 2 és 4 éves generációs hosszal kombinálva), és így kaptunk egy durva becslést a domesztikációs időszakra (lásd a 2. kiegészítő ábrát).

Genetikai terhelés

A genetikai terhelés becslése kétféle módon történt. Először is úgy, hogy a genetikai terhelést minden egyes egyedre a káros fitneszhatások összegeként számoltuk ki az összes fehérjekódoló genomi pozícióra Librado et al.60 módszerét követve. Röviden, az evolúciós kényszer helyettesítőjeként a PhyloP-pontszámokat használtuk a 46-os emlős-illesztésből származó PhyloP-pontszámokat (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Ebből az összehangolásból azonosítottuk a funkcionális korlátok között fejlődő fehérjekódoló helyeket (phyloP score ≥1,5). Ezután minden egyes Ovis vagy Capra genom esetében megvizsgáltuk, hogy ezek a helyek mutálódtak-e. Ha igen, akkor a phyloP-pontszámokat összesítettük az összes mutálódott helyen, így az erősen korlátozott helyeken bekövetkezett mutációk arányosan nagyobb mértékben járulnak hozzá a teljes terhelésbecsléshez. Így kaptunk egy terhelésbecslést minden egyes juh/kecske genomra. Végül, hogy megkapjuk az oldalankénti átlagos terhelést, elosztottuk az elemzett pozíciók teljes számával. Érdemes megjegyezni, hogy a homozigóta helyeket kondicionáltuk, hogy elkerüljük a heterozigóta helyeken lévő mutációk dominancia együtthatójának modellezését (pl. recesszív, intermedier, domináns). Másodszor, génről génre összehasonlítottuk a genetikai károsító terhelést a vadon élő és a háziasított Ovis-csoportokban Wilcoxon-teszt elvégzésével, azzal az alternatív hipotézissel, hogy a háziasított állatokban nagyobb a terhelés, mint a vadon élő rokonokban. p-értékeket korrigáltuk a többszörös tesztelésre61 , és a korrigált p-értékek < 0,05 küszöbértékét alkalmaztuk. A genetikai terhelés szignifikáns növekedését mutató génkészleten génontológiai gazdagodási elemzést végeztünk a WebGestalt62,63 segítségével. Mivel a referencia genomok gyengén annotáltak a gének tekintetében, a fajunk, valamint az ember és az egér közötti egykópiás ortológokra támaszkodtunk. Az X-kromoszómáról származó géneket kizártuk a háttérkészletből. A Capra esetében nem végeztük el ezt az elemzést a vad mintákban megfigyelt nagyobb mértékű beltenyésztés miatt.

Szelekciós szignatúrák kimutatása

A domesztikációval kapcsolatos szelekciós szignatúrák kimutatásához minden olyan biallelikus SNP-t felhasználtunk, amely a három vizsgált csoport (azaz az iráni és a marokkói házi csoport, valamint a vad csoport minden nemzetség esetében) közül legalább az egyikben 0,10-nél nagyobb minor allélfrekvenciát mutatott. Mivel azt vártuk, hogy a háziasítási folyamathoz kapcsolódó szelekciós jelek minden háziállatban jelen lesznek, a következő általános stratégiát alkalmaztuk: a hapFLK29 segítségével (lásd az 5. kiegészítő megjegyzést és a 9., 10. és 11. kiegészítő ábrát) minden egyes nemzetség esetében a vad csoportot a hagyományosan kezelt házi csoportok (azaz az iráni és a marokkói) mindegyikével szemben teszteltük, és azokra a feltételezhetően szelekció alatt álló közös régiókra összpontosítottunk, amelyeket mindkét esetben kimutattunk. A csoportminták mérete (n = 13-20) megfelelt a módszer követelményeinek29. Vizuálisan ellenőriztük, hogy a hapFLK-val talált konzisztens szelekciós jelek jelen vannak-e az egyes nemzetségek megfelelő világpanelkészletében is, de ezeket a csoportokat nem vontuk be a statisztikai tesztbe a több fajtából álló összetételük miatt. Végül az Ovis és a Capra között a szelekció közös jeleit kerestük rétegzett FDR megközelítéssel. A stratégiát a 4. kiegészítő ábra mutatja be.

HapFLK-teszteket végeztünk a vad csoport és az egyes nemzetségek iráni és marokkói csoportjainak szembeállítására. A rokonsági mátrixot a csoportpárok közötti Reynold-féle genetikai távolságokból64 számoltuk ki, a variánsok egy százalékának véletlenszerű részhalmazát felhasználva. A levezetett populációfát a neighbor-joining algoritmus segítségével építettük fel. Minden egyes SNP esetében elvégeztük a hapFLK-tesztet, amely a haplotípusos információkat is beépíti, hogy növelje a szelektív seprések kimutatásának erejét. A hapFLK-statisztika minden egyes vizsgált SNP esetében kiszámította a haplotípusos gyakoriságok eltérését a rokonsági mátrix65 által becsült semleges modellhez képest. A kapcsolási kiegyensúlyozatlansági információk kiaknázásához a hapFLK a Scheet és Stephens66 többpontos modelljét használja a többfókuszú genotípusokra, amely a fázis nélküli adatokra illeszthető. E modell egyik fő alkalmazása a fázisbecslés elvégzése (fastPHASE szoftver66). Elemzésünkben a modellt fázis nélküli adatokon képeztük ki, ezért elemzésünk figyelembe veszi a fázis bizonytalanságát. A módszert a helyi haplotípusok kromoszómák mentén történő átcsoportosítására használtuk egy Rejtett Markov-modell segítségével egy meghatározott számú, K 25-re beállított klaszterbe.

A két hagyományosan kezelt hazai csoportban minden nemzetség esetében a feltételezhetően szelekció alatt álló közös régiók azonosításához egyesítettük a két korábbi hapFLK-elemzést. Mindegyik elemzéshez a hapFLK-pontszámokat χ2-eloszlásra illesztettük, hogy p-értékeket kapjunk (a szkript elérhető https://forge-dga.jouy.inra.fr/projects/hapflk/documents). A vadon élő csoport és az egyes háziállat-csoportok közötti két kontraszt eredményeit a Stouffer-módszer67 segítségével kombináltuk, hogy egyetlen p-értéket kapjunk a vadon élő és a háziállatok összehasonlítására. Végül az FDR-keretet68 alkalmaztuk az SNP-k teljes halmazára, hogy a kombinált p-értékeket q-értékekké alakítsuk át. A < 10-2 q-értéket mutató SNP-ket megtartottuk és genomi régiókba csoportosítottuk, ha azok egymástól kevesebb mint 50 kb távolságra voltak egymástól.

Az Ovis és a Capra közötti szelekciós jel megosztottságának vizsgálatához először a két referencia genom keresztillesztését használtuk a homológ szakaszok azonosítására. Ezután rétegzett FDR keretrendszert alkalmaztunk69. Ez a megközelítés azon alapul, hogy a genetikai adatok előzetes információinak ismeretében a tesztekben eredendően rétegződés van69 , mivel az igaz alternatív hipotézisek mögöttes eloszlása a különböző genomi régiók eltérő dinamikája szerint eltérő lehet, ami a p-értékek eltérő eloszlásához vezet. Ez megköveteli, hogy az FDR-rel korrigált p-értékeket (azaz q-értékeket) külön-külön kapjuk meg a különböző rétegekre. Az egyes nemzetségekben úgy kerestük a konvergenciákat, hogy szétválasztottuk a másik nemzetségben kimutatottakkal homológ régiókat (a továbbiakban: közös réteg) és a genom többi részét (a továbbiakban: általános réteg). A két meghatározott rétegre külön-külön kivontuk a p-értékeket, majd az FDR keretrendszer segítségével q-értékeket számoltunk. Ezek a rétegzett q-értékek voltak a statisztikai szignifikancia szempontjából figyelembe vett végső mennyiségek (<10-2) a szelekció alatt álló SNP-k kimutatásához és a megfelelő genomi régiókba való beillesztéséhez.

A vadon élő és háziállatokat a két nemzetségben megkülönböztető szelekció konvergens szignatúráinak tesztelésére megvizsgáltuk a kapcsolatot a q-értékekre alkalmazott szignifikancia küszöbérték (amelyet 0,2 és 0,002 között változtattunk) között az egyik nemzetségben és annak becsült valószínűsége között, hogy egy SNP a másik nemzetség közös rétegében szelektálódik, a Storey et al.70 segítségével. megközelítéssel. A q-értékre alkalmazott küszöbérték csökkenésével (a szigor növekedésével) a levezetett valószínűség növekedése azt jelzi, hogy minél szignifikánsabb a régió az egyik nemzetségben, annál valószínűbb, hogy a másik nemzetségben is találunk szignifikáns SNP-ket.

A három hazai csoport között nem konzisztens szelekciós jeleket kiszűrtük. Minden egyes kimutatott régióhoz az egyedek fázisos haplotípusait használtuk, amelyeket a szekvenciák közötti azonosság százalékos aránya alapján Neighbor-Joining fák segítségével klasztereztünk. Csak a konzisztens jeleket mutató régiókat tartottuk meg (5. kiegészítő ábra).

Azért, hogy következtetni tudjunk arra, hogy a hapFLK-val észlelt szelekciós jelek a szelekció lazulását vagy pozitív szelekciót jeleznek-e a háziasított egyedeknél, megbecsültük a nukleotiddiverzitás különbségét (π) minden egyes feltételezett szelekció alatt álló régióban a vadon élő és a háziasított csoportok között. Ezt a különbséget a Δπ indexként fejeztük ki, amelyet minden egyes genomi régióra a vadon élő csoportra számított π és az iráni és marokkói házi csoportokra számított π átlagának különbségeként számoltunk ki, mínusz a két csoport közötti π különbség, amelyet a teljes genomra számítottunk:

Egy negatív érték azt jelezné, hogy a nukleotiddiverzitás alacsonyabb a vad csoportban a két hazai csoport átlagához képest, és úgy tekinthetnénk, hogy ez utóbbi csoportokban a szelekció lazulását, a hazaiaknál diverzifikáló szelekciót, a vadaknál pozitív szelekciót mutat. Ezzel szemben a pozitív érték a háziasított csoportokban bekövetkezett irányított pozitív vagy stabilizáló szelekcióra utalna. A haplotípus-klaszterezést arra is használtuk, hogy manuálisan ellenőrizzük az egyes régiókban, hogy az észlelt szelekciós söprés megerősíti-e a Δπ-index által adott jelzéseket.

Funkcionális értelmezéseket végeztünk az alábbiak szerint. Minden egyes szelekció alatt álló régió esetében a funkcionális szerepek azonosításához figyelembe vettük a régiót plusz 50 kb-t mindkét oldalon, valamint 5 kb-t a gének előtt és után, és értékeltük az átfedést e koordináták között az érdeklődésre számot tartó gének megtartása érdekében. Végül úgy tekintettük, hogy egy gén akkor kapcsolódik egy adott detektált régióhoz, ha a régió és a gén pozíciói átfedésben vannak. Ezután felmértük, hogy melyik gén a legvalószínűbb célpontja a szelekciónak, figyelembe véve a legfelső jelhez legközelebbi gént, azaz a régióban a legalacsonyabb q-értékkel rendelkező pozíciót. A géneket funkcionálisan annotáltuk az Uniprot (http://www.uniprot.org/) segítségével, figyelembe véve a “Biológiai folyamat” kategória (azaz GO:0008150) 30 gyermekterminusában való részvételüket (azaz a terminusok közvetlen leszármazottait). A 33 kategóriából 30 esetében az egyes géneknek megfelelő összes GO-terminust (Kiegészítő adatok 4) visszakértük, mivel három olyan terminust nem vettünk figyelembe, amelyek nem érintettek emlősök funkcióiban (pl. GO:0006791 kénhasznosítás, GO:0006794 foszforhasznosítás, GO:0015976 szénhasznosítás). Két χ2-tesztet végeztünk a GO-kategóriákba tartozó gének eloszlásainak összehasonlítására, azaz (i) a nemzetségspecifikus régiókból származó szelekció alatt álló gének és a homológ régiókból származó gének, valamint (ii) az összes szelekció alatt álló gén és a Swiss-Protban található 18 689, GO-terminusokhoz társított emberi gén összehasonlítására. Annak érdekében, hogy a gének funkcióit állattartási kontextusban értelmezzük, a fenotípusos hatásaikról az irodalomból elérhető információkat is lekérdeztük.

Végezetül, hogy megtaláljuk azokat az SNP-ket a korábban detektált régiókon belül, amelyek a leginkább megkülönböztetik a vadon élő és a háziasított csoportokat, az FLK-statisztikát használtuk. A hapFLK-hoz hasonlóan az egymarker allélfrekvenciák eltérését jelenti a rokonsági mátrix65 által becsült semleges modellhez képest. Ugyanazzal az eljárással illesztettük a két elemzésből származó pontszámokat egy χ2 eloszláshoz, és kombináltuk a kapott p-értékeket, mint a hapFLK-teszt esetében. A p-értékek nem egyenletes eloszlása azonban kizárta az FDR-keret alkalmazását, és a hapFLK-val kimutatott régiókon belül olyan SNP-ket választottunk ki, amelyek p-értéke <10-4 volt. Ezekhez az SNP-khez a Variant Effect Predictor (VEP) annotációkat71 használtuk, amelyeket az Ensembl v74 juh OARv3.1 genom annotációjából generáltunk az Ovis esetében (http://www.ensembl.org/Ovis_aries/Tools/VEP) és a kecske CHIR1.0 genom annotációjából, amelyet az NCBI eukarióta genom annotációs csővezetéke készített a Capra esetében (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). Az SNP-ket intergenikus, upstream és downstream (beleértve az UTR-eket is), valamint intronikus és exonikus pozíciókba soroltuk. A <10-4 FLK p-értékkel rendelkező SNP-k és a szelekciós szignatúrák kimutatására használt összes SNP eloszlása közötti különbségeket χ2-próbával vizsgáltuk.

Adatok elérhetősége

Az elemzésekben felhasznált 73 Ovis- és 72 Capra-mintára generált szekvenciák és metaadatok nyilvánosan elérhetők. Általános információk és az összes vcf-fájl megtalálható az Ensembl weboldalán (http://projects.ensembl.org/nextgen/). Az O. orientalis és a C. aegagrus összes Fastq fájlja, Bam fájlja és de novo összeállítása megtalálható az Európai Nukleotid Archívumban (https://www.ebi.ac.uk/ena) a Nextgen projekt csatlakozási kódja alatt (PRJEB7436).

.