Sampling
Domesticerade får (O. aries) och getter (C. hircus) provtogs i Iran (IROA- respektive IRCH-grupperna) och Marocko (MOOA- respektive MOCH-grupperna) med sammanlagt 20 djur per grupp (kompletterande fig. 6). Dessa prover samlades in mellan januari 2008 och mars 2012 i den norra delen av Marocko och mellan augusti 2011 och juli 2012 i nordvästra Iran inom ramen för det europeiska Nextgen-projektet (bidragsöverenskommelse nr 244356) i enlighet med de etiska bestämmelserna i Europeiska unionens direktiv 86/609/EEG. Öronklipp samlades in från den distala delen av örat på slumpmässigt utvalda djur och förvarades omedelbart i 96 % etanol i ett dygn innan de överfördes till kiselgelpärlor fram till DNA-extraktion.
De vilda arterna asiatisk mufflon (O. orientalis) och bezoar ibex (C. aegagrus) provtogs i nordvästra Iran inom ramen för domesticeringens vagga21,22. Tretton vävnader från asiatiska mufflondjur och 18 vävnader från Bezoar ibex (IROO- respektive IRCA-grupperna, kompletterande figur 6) samlades in antingen från djur i fångenskap eller från nyligen jagade djur och från frysta prover som fanns tillgängliga på det iranska miljödepartementet. Denna individbaserade provtagningsmetod är utformad för att minimera potentiell bias genom att undvika överrepresentation av lokala effekter (t.ex, lokal inavel).
Tilläggsdata
Tilläggsvis sammanställdes en världsomspännande raspanel för får och getter (wpOA respektive wpCH). wpOA omfattade 20 prover för omsekvensering av hela arvsmassan (WGS) med 12x täckning som representerar 20 olika världsomspännande raser och som tillhandahölls av International Sheep Genome Consortium. wpCH bestod av 14 WGS-prover som sekvenserades med 12x täckning och som representerade 9 europeiska individer, dvs, 2 franska Alpine- och 2 franska Saanen-prover som sekvenserats av INRA, 5 italienska Saanen-prover som tillhandahållits av Parco Tecnologico Padano, och 5 australiensiska individer, dvs, 2 Boer-, 2 Rangeland- och 1 Cashmere-prov som tillhandahölls av CSIRO (Supplementary Data 5).
Produktion av WGS-data
Genomiskt DNA extraherades framgångsrikt från alla vävnadsprov med hjälp av Macherey Nagel NucleoSpin 96 Tissue kit, med anpassning av tillverkarens protokoll. Vävnadsprovtagningen utfördes i MN-kvadratbrunnsblock för att erhålla 25 mg fragment per prov. Tre och ett halvt MN kvadrat-96 block förbereddes och extraktionen utfördes med hjälp av en Tecan Freedom Evo Liquid handler enligt tillverkarens protokoll. Ett förlysningssteg utfördes för att homogenisera proverna med 180 µl T1-buffert och 25 µl proteinas K över natten vid 56 °C. För att justera bindningsförhållandena tillsattes 200 µl BQ1-buffert och provplattan inkuberades 1 timme vid 70 °C. Därefter tillsattes 200 µl 100 % etanol. Lysaten överfördes till Nucleospin Tissue-bindningsplatta och vakuum (-0,2 bar, 5 min) användes för att avlägsna genomströmningen. Tre tvättsteg gjordes med BW- respektive B5-buffer, och vakuum applicerades igen för att avlägsna flödet. Före eluering av genomiskt DNA torkades ett silikamembran från Nucleospin Tissue binding plate under vakuum med minst -0,6 bar i 10 minuter. Elueringssteget utfördes med 100 µl förvärmd BE-buffert (70 °C) och ett centrifugeringssteg vid 3700 rcf i 5 minuter i 96-PCR-brunnar. Genomiskt DNA förvarades vid 4 °C för att undvika frysning och upptining och testades för koncentration (som ng/μl) med hjälp av Picogreen-metoden och med hjälp av en Nanodrop.
Helgenomerna omsekvenserades från 500 ng genomiskt DNA som klipptes till 150-700 bp med hjälp av Covaris® E210-instrumentet för varje prov och användes för Illumina®-biblioteksframställning genom ett halvautomatiserat protokoll. Slutreparation, A-tailing och Illumina®-kompatibla adaptrar (BioScientific) ligering utfördes med hjälp av SPRIWorks Library Preparation System och SPRI TE-instrumentet (Beckmann Coulter) enligt tillverkarens protokoll. Ett storleksurval på 300-600 bp tillämpades för att återvinna de flesta fragmenten. DNA-fragmenten amplifierades genom 12 PCR-cykler med hjälp av Platinum Pfx Taq Polymerase Kit (Life® Technologies) och Illumina® adapterspecifika primers. Biblioteken renades med 0,8x AMPure XP-beads (Beckmann Coulter) och analyserades med Agilent 2100 Bioanalyzer (Agilent® Technologies) och qPCR-kvantifiering. Biblioteken sekvenserades med hjälp av 100 baslängders läsningskemi i parvisa flödesceller på Illumina® HiSeq2000.
Illuminas parvisa läsning för Ovis kartlades till fårreferensgenomet (OAR v3.1, GenBank assembly GCA_000298735.146) och för Capra till getreferensgenomet (CHIR v1.0, GenBank assembly GCA_000317765.147), med hjälp av BWA-MEM48. Den BAM-fil som producerades för varje individ sorterades med hjälp av Picard SortSam och förbättrades med hjälp av sekventiellt Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator och GATK IndelRealigner49 samt SAMtools calmd50.
Variantupptäckt gjordes med hjälp av tre olika algoritmer: Samtools mpileup50, GATK UnifiedGenotyper51 och Freebayes52. Variantplatser identifierades oberoende av varandra för var och en av de sex grupperna, med hjälp av algoritmernas multiproverläge: (i) 162 prover från MOOA, (ii) 20 prover från IROA, (iii) 14 prover från IROO, (iv) 162 prover från MOCH, (v) 20 prover från IRCH, (vi) 19 prover från IRCA. För vissa grupper fanns WGS från fler individer tillgängliga som en del av NextGen-projektet (se ovan). De prover som användes i den aktuella studien valdes ut för att i möjligaste mån få balanserade grupper om 20 individer. För IRCA- och IROO-grupperna blev ytterligare prover tillgängliga i ett senare skede och lades till för nedströmsanalyser. Djur med låg anpassning och kallelsekvalitet togs bort för att erhålla den slutliga datamängden (Supplementary Data 5).
Inom varje grupp genomfördes två på varandra följande omgångar av filtrering av variantplatsernas kvalitet. I filtreringssteg 1 slogs anrop samman från de tre algoritmerna, samtidigt som man filtrerade bort de anrop som hade lägst konfidens. En variantplats godkändes om den anropades av minst två olika anropningsalgoritmer med phred variantkvalitet >30. En alternativ allel på en plats godkändes om den anropades av någon av anropningsalgoritmerna och genotypantalet var >0. I filtreringssteg 2 användes Variant Quality Score Recalibration by GATK. Först genererade vi en träningsuppsättning av de variantplatser med högsta konfidens inom gruppen där (i) platsen kallas av alla tre variantkallare med phred-skalad variantkvalitet >100, (ii) platsen är biallelisk, (iii) minor allelantalet är minst 3 samtidigt som endast prover med genotyp phred-skalad kvalitet >30 räknas med. Träningsuppsättningen användes för att bygga en Gauss-modell med hjälp av verktyget GATK VariantRecalibrator med följande variantannotationer från UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. En gaussisk modell tillämpades på hela datamängden och genererade en VQSLOD (log oddskvot för att vara en sann variant). Platser filtrerades om VQSLOD <gränsvärde. Cutoff-värdet fastställdes för varje grupp enligt följande: Minsta VQSLOD = {medianvärdet av VQSLOD för varianter i träningsuppsättningen}-3 × {medianvärdet av den absoluta avvikelsen VQSLOD för varianter i träningsuppsättningen}. Förhållandet mellan SNP-övergångar och SNP-transversioner tyder på att det valda avgränsningskriteriet ger den bästa balansen mellan selektivitet och känslighet.
SNP-utlysningsuppsättningar för sex grupper av Ovis- och Capra-djur genererades (dvs. iranska och marockanska husdjur och vilda djur för varje släkte). Eftersom de analyser som utfördes i den här studien krävde jämförelser mellan grupperna skapade vi genotypsuppsättningar vid en konsekvent uppsättning SNP-platser för alla djur från varje grupp. För varje släkte sammanfogade vi de alternativa utlysningsplatserna från dess tre grupper och behöll endast bialleliska positioner utan saknade data. Genotyperna återkallades vid varje biallelisk SNP-plats för alla individer av intresse av GATK UnifiedGenotyper med hjälp av alternativet GENOTYPE_GIVEN_ALLELES. I detta skede utökades listan över individer till att omfatta djur som tillhör de globala raspanelerna för får och getter (wpOA och wpCH) och ytterligare vilda prover som blev tillgängliga i detta skede (4 O. orientalis och 4 C. aegagrus). Genotyperna förbättrades och fasades in i grupperna med Beagle 453 och filtrerades sedan bort där genotypsannolikheten var mindre än 0,95. Slutligen filtrerade vi bort platser som var monomorfa i de olika undergrupperna av individer som användes i denna studie (se nedan).
För att jämföra de signaler av selektion som upptäcktes mellan Ovis och Capra utförde vi en korsjustering mellan de två referensgenomerna. Först använde vi pipeline för parvis anpassning från kodbasen Ensembl release 6954 för att anpassa referensgenomerna för får (OARv3.1) och get (CHIR1.0). Denna pipeline använder LastZ55 för att anpassa på DNA-nivå, följt av en efterbearbetning där anpassade block kedjas ihop enligt deras placering i de båda genomerna. LastZ-pipeline för parvis anpassning körs rutinmässigt av Ensembl för alla arter som stöds, men geten ingår ännu inte i Ensembl. För att undvika bias till förmån för någon av arterna producerade vi två olika interspecifika anpassningar. I den ena användes får som referensgenom och get som icke-referens, i den andra användes get som referensgenom och får som icke-referens. Skillnaden är att genomiska regioner hos referensarten tvingas kartlägga unikt till en enda loci hos icke-referensarten, medan genomiska regioner hos icke-referensarten tillåts kartlägga till flera platser hos referensarten. För segment av kromosomer i ett referensgenom erhöll vi koordinaterna på icke-referensgenomet. Slutligen, för de SNPs som upptäcktes i ett släkte använde vi hela genomutjämningen med referensgenomet i det andra släktet för att identifiera motsvarande positioner (kompletterande tabell 6).
Genisk struktur
För att beskriva den genetiska mångfalden inom grupperna använde vi VCFtools56 för att beräkna sammanfattande statistik över den genetiska variationen på de 73 individerna för Ovis (dvs, 13 IROO, 20 IROA, 20 MOOA och 20 wpOA) och 72 individer för Capra (dvs. 18 IRCA, 20 IRCH, 20 MOCH och 14 wpCH). Den statistik som mättes var det totala antalet polymorfa varianter (S) för hela uppsättningen individer i varje släkte och inom varje grupp, den genomsnittliga nukleotiddiversiteten (π) inom varje grupp och inavelskoefficienten (F) för varje individ. Inom varje släkte testades skillnaderna mellan den vilda gruppen och varje domesticerad grupp med hjälp av ett ensidigt t-test för individuella värden för inavel och genetisk belastning, och ett tvåsidigt Mann-Whitney-test för nukleotiddiversitet per plats.
Den totala divergensen mellan de fyra grupperna inom varje släkte (dvs, vilda, iranska och marockanska husdjur och världspanelen) uppskattades med hjälp av alla bialleliska SNPs och den genomsnittliga viktade parvisa Fst enligt Weir och Cockerham57 som implementerades i VCFtools56. Den genetiska strukturen mellan grupperna bedömdes med klustermetoden sNMF26 , efter beskärning av datamängden för att ta bort SNP:er med kopplingsobalans (r²) större än 0,2 med hjälp av VCFtools. Kopplingsobalansen (r²) beräknades mellan par av SNP:er inom glidande fönster med 50 SNP:er, med en SNP per par som slumpmässigt togs bort när r² var större än 0,2. För varje sNMF-analys utfördes fem körningar med samma antal kluster (K) med värden för K från 1 till 10. Vi använde korsentropikriteriet för att identifiera den mest sannolika klusterlösningen, men alternativa partitioner för olika antal K undersöktes också för att bedöma hur individer delades upp mellan kluster.
För att skilja mellan delad härstamning och blandning körde vi TreeMix27 för att gemensamt skatta befolkningsuppdelningar och efterföljande blandningshändelser med hjälp av den beskurna datamängd som användes för sNMF. Vi körde TreeMix med alternativet -global för att förfina våra slutsatser om maximal sannolikhet. Vi rotade TreeMix-trädet med uppdelningen mellan vilda och tama individer. Blockstorleken för jackknifing var -k 500 SNPs, vilket ungefär motsvarar 150 kb, vilket överstiger de genomsnittliga block av LD som finns hos både får och getter. Vi genererade ett Maximum Likelihood-träd utan migration och lade sedan till migrationshändelser och undersökte den inkrementella förändringen i den varians som förklaras av modellen och restvärdena mellan individer. Målet var att upptäcka eventuella höga restvärden eller migrationskanter mellan vilda och tama individer. För att ytterligare utforska den statistiska relevansen av möjliga blandningsvektorer som identifierats av TreeMix (kompletterande tabell 3) beräknade vi trepopulationstestet f328 som ett formellt test av genetisk introgression, med hjälp av programmet qp3Pop i ADMIXTOOLS-sviten58 för varje kombination av grupper. För Capra delades wpCH-gruppen upp mellan australiska raser, franska raser och italienska raser. Resultaten redovisas i Supplementary Data 2.
Demografisk inferens
För varje släkte utförde vi analyser av demografisk inferens från förfäderna med hjälp av MSMC-modellen som implementerats i programvaran MSMC225. MSMC bygger på den parvisa sekventiella Markovian coalescent59; den använder dock haplotyper av fasade genomsekvensdata som indata. För varje analys använde vi två individer från en grupp, alltså 4 haplotyper. Varje analys upprepades för en annan slumpmässig uppsättning av två individer, dvs. en replik av analysen per grupp. Ingångs- och utdatafiler genererades och analyserades med de pythonskript som tillhandahålls med MSMC-programvaran och som finns på https://github.com/stschiff/msmc-tools. Analysparametrarna behölls som standardparametrar, utom mutationshastigheten som sattes till 2,5×10-8 och generationslängden sattes till 2 år. För att uppskatta osäkerheten i tidsuppskattningarna varierade vi dessa parametrar (mutationshastighet på 2,5×10-8 och 1,0×10-8 i kombination med generationslängd på 2 och 4 år) och gav en grov uppskattning av domesticeringsperioden (se kompletterande figur 2).
Genetisk belastning
Den genetiska belastningen uppskattades på två sätt. För det första genom att beräkna den genetiska belastningen för varje individ som summan av skadliga fitnesseffekter över alla proteinkodande genomiska positioner enligt Librados et al.60 metod. I korthet använde vi PhyloP-poängen från den 46-vägiga däggdjursanpassningen (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/) som en proxy för evolutionär begränsning. Från denna anpassning identifierade vi proteinkodande platser som utvecklas under funktionella begränsningar (PhyloP-poäng ≥1,5). För varje Ovis- eller Capra-genom undersökte vi sedan om dessa platser var muterade. Om så var fallet summerade vi phyloP-poängen för alla muterade platser, så att mutationer på platser med stora begränsningar bidrar proportionellt sett mer till uppskattningen av den totala belastningen. Detta gav en skattning av belastningen för varje får- och getgenom. För att få fram en genomsnittlig belastning per webbplats dividerade vi slutligen med det totala antalet analyserade positioner. Det är värt att notera att vi betingade oss på homozygota platser för att undvika att modellera dominanskoefficienten för mutationer på heterozygota platser (t.ex. recessiv, intermediär, dominant). För det andra jämförde vi gen för gen den genetiska skadliga belastningen i vilda och domesticerade Ovis-grupper genom att utföra ett Wilcoxon-test, med den alternativa hypotesen att de domesticerade djuren har större belastning än vilda släktingar. p-värdena korrigerades för multipel testning61 och vi tillämpade ett tröskelvärde för justerade p-värden < 0,05. Vi utförde en genontologisk anrikningsanalys på den uppsättning gener som visade en signifikant ökning av den genetiska belastningen med hjälp av WebGestalt62,63. Eftersom referensgenomerna är dåligt annoterade för gener förlitade vi oss på ortologer med en enda kopia mellan vår art och människa och mus. Gener från X-kromosomen uteslöts från bakgrundsmängden. Vi utförde inte denna analys på Capra på grund av den högre inavel som observerades i de vilda proverna.
Detektering av urvalssignaturer
För att upptäcka signaturer av urval relaterat till domesticering använde vi alla bialleliska SNP:er som uppvisade en mindre allelfrekvens som var större än 0,10 i minst en av de tre grupperna som testades (dvs. de iranska och marockanska domesticeringsgrupperna, samt den vilda gruppen för varje släkte). Eftersom vi förväntade oss att signaturer av selektion i samband med domesticeringsprocessen skulle finnas hos alla husdjur antog vi följande allmänna strategi: vi testade med hapFLK29 (se kompletterande anmärkning 5 och kompletterande figurer 9, 10 och 11) för varje släkte den vilda gruppen mot var och en av de traditionellt hanterade husdjursgrupperna (dvs. iranska och marockanska) och inriktade oss på de gemensamma regioner som antas vara föremål för selektion och som upptäcktes i båda fallen. Gruppernas provstorlekar (n = 13-20) var förenliga med metodens krav29. Vi kontrollerade visuellt om de konsekventa signaturer av selektion som hittades med hapFLK också fanns i motsvarande världspaneluppsättning för varje släkte, men inkluderade inte dessa grupper i det statistiska testet på grund av deras sammansättning med flera raser. Slutligen letade vi efter gemensamma signaler av selektion mellan Ovis och Capra med hjälp av en stratifierad FDR-metod. Strategin visas i kompletterande figur 4.
Vi utförde hapFLK-tester för att kontrastera den vilda gruppen mot var och en av de iranska och marockanska grupperna i varje släkte. Släktskapsmatrisen beräknades från Reynolds genetiska avstånd64 mellan par av grupper, med hjälp av en slumpmässig delmängd av en procent av varianterna. Det härledda populationsträdet byggdes upp med hjälp av algoritmen neighbor-joining. För varje SNP utförde vi hapFLK-testet som innehåller haplotypisk information för att öka kraften för att upptäcka selektiva svepningar. För varje testad SNP beräknade hapFLK-statistiken avvikelsen av haplotypiska frekvenser i förhållande till den neutrala modell som uppskattas av släktskapsmatrisen65. För att utnyttja information om kopplingsdifferens använder hapFLK Scheet och Stephens’66 multipunktsmodell för genotyper på flera platser som kan anpassas till ofasade data. En av de viktigaste tillämpningarna av denna modell är att utföra fasuppskattning (programvaran fastPHASE66). I vår analys har modellen tränats på ofasade data, och därför tar vår analys hänsyn till fasosäkerhet. Metoden användes för att omgruppera lokala haplotyper längs kromosomer i ett specificerat antal kluster K som sattes till 25, med hjälp av en Hidden Markov-modell.
För att identifiera de gemensamma regioner som antas vara föremål för selektion i de två traditionellt skötta inhemska grupperna för varje släkte kombinerade vi de två tidigare hapFLK-analyserna. För varje analys anpassades hapFLK-poängen till en χ2-fördelning för att få fram p-värden (skriptet finns tillgängligt på https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Resultaten av de två kontrasterna mellan den vilda gruppen och var och en av husdjursgrupperna kombinerades med Stouffers metod67 för att få fram enskilda p-värden för jämförelsen mellan vilda och tama djur. Slutligen tillämpades FDR-ramen68 på hela uppsättningen SNPs för att omvandla de kombinerade p-värdena till q-värden. SNPs som visade q-värden < 10-2 behölls och grupperades i genomiska regioner när de låg mindre än 50 kb från varandra.
För att undersöka om selektionssignalen delades mellan Ovis och Capra använde vi först korsanpassningen av de två referensgenomerna för att identifiera homologa segment. Vi tillämpade sedan en stratifierad FDR-ram69. Detta tillvägagångssätt bygger på det faktum att det finns en inneboende stratifiering i testerna med tanke på förhandsinformationen i de genetiska uppgifterna69 , eftersom den underliggande fördelningen av sanna alternativa hypoteser kan vara olika beroende på olika dynamik i olika genomiska regioner, vilket leder till olika fördelningar av p-värden. Detta kräver att man får FDR-justerade p-värden (dvs. q-värden) separat för de olika strata. Vi sökte efter konvergenser i varje släkte genom att separera de regioner som är homologa med dem som upptäckts i det andra släktet (kallat det delade stratumet) och resten av genomet (kallat det allmänna stratumet). Vi extraherade p-värdena separat för vart och ett av de två definierade strata och beräknade sedan q-värdena med hjälp av FDR-ramen. Dessa stratifierade q-värden var de slutliga kvantiteter som beaktades för statistisk signifikans (<10-2) för att upptäcka SNP:er som är föremål för urval och slå samman dem i motsvarande genomiska regioner.
För att testa konvergerande signaturer av selektion som skiljer vilda från tama djur i de båda släktena undersökte vi förhållandet mellan den signifikanströskel som tillämpas på q-värden (som vi lät variera från 0,2 till 0,002) i ett släkte och den uppskattade sannolikheten för att en SNP är selekterad i det delade stratumet i det andra släktet med hjälp av Storey et al.70. tillvägagångssätt. En ökning av den härledda sannolikheten med en minskning av det tröskelvärde som tillämpas på q-värdet (ökad stringens) indikerar att ju mer signifikant regionen är i ett släkte, desto mer sannolikt är det att vi skulle hitta signifikanta SNP:er i det andra släktet.
Vi filtrerade bort de selektionssignaler som inte var konsekventa mellan de tre inhemska grupperna. För varje upptäckt region använde vi de fasade haplotyperna hos individer som klustrades med hjälp av Neighbor-Joining-träd baserat på den procentuella identiteten mellan sekvenserna. Endast regioner som visade konsekventa signaler behölls (kompletterande figur 5).
För att kunna härleda om de selektionssignaler som upptäcktes med hapFLK indikerade avspänning av selektion eller positiv selektion hos husdjuren uppskattade vi skillnaden i nukleotiddiversitet (π) på varje förmodad region under selektion mellan de vilda och husdjursgrupperna. Vi uttryckte denna skillnad som Δπ-indexet, som beräknades för varje genomisk region som skillnaden mellan π beräknat för den vilda gruppen och genomsnittet av π för de iranska och marockanska domesticerade grupperna, minus skillnaden i π mellan dessa två grupper beräknat över hela genomet:
Ett negativt värde skulle indikera att nukleotiddiversiteten är lägre i den vilda gruppen jämfört med genomsnittet av de två inhemska grupperna, och skulle anses visa på en avmattning av urvalet i de sistnämnda grupperna, diversifierande urval i de inhemska eller positivt urval i de vilda. Omvänt skulle ett positivt värde indikera ett riktat positivt eller stabiliserande urval som skett i de inhemska grupperna. Vi använde också haplotypklustering för att manuellt verifiera i varje region om den upptäckta selektiva svepningen bekräftade de indikationer som gavs av Δπ-indexet.
Vi genomförde funktionella tolkningar på följande sätt. För varje region under urval betraktade vi regionen plus 50 kb på varje sida för att identifiera funktionella roller och 5 kb uppströms och nedströms från generna och vi bedömde överlappningen mellan dessa koordinater för att behålla de intressanta generna. Slutligen ansåg vi att en gen var relaterad till en given upptäckt region när regionens och genens positioner överlappade varandra. Vi bedömde sedan vilken gen som sannolikt var den mest sannolika måltavlan för urvalet genom att ta hänsyn till den gen som ligger närmast den högsta signalen, dvs. positionen för det lägsta q-värdet inom regionen. Generna annoterades funktionellt med hjälp av Uniprot (http://www.uniprot.org/), genom att beakta deras inblandning i 30 underordnade termer (dvs. termernas direkta nedstigningar) i kategorin ”Biologisk process” (dvs. GO:0008150). Vi hämtade alla GO-termer som motsvarar varje gen (Supplementary Data 4) för 30 av de 33 kategorierna, eftersom vi inte tog hänsyn till tre termer som inte var involverade i däggdjursfunktioner (dvs. GO:0006791 svavelutnyttjande, GO:0006794 fosforutnyttjande, GO:0015976 kolutnyttjande). Vi utförde två χ2-tester för att jämföra fördelningen av gener i GO-kategorierna, dvs. (i) gener under urval från genusspecifika regioner jämfört med gener från homologa regioner och (ii) alla gener under urval jämfört med de 18 689 mänskliga gener som är associerade med GO-termer i Swiss-Prot. För att tolka genernas funktioner i ett djursammanhang hämtade vi också den information som fanns tillgänglig i litteraturen om deras fenotypiska effekter.
För att slutligen hitta de SNPs inom de tidigare upptäckta regionerna som var de som skiljde sig mest mellan vilda och tama grupper använde vi FLK-statistiken. När det gäller hapFLK representerar den avvikelsen av allelfrekvenser för enskilda markörer med avseende på den neutrala modellen som uppskattas av släktskapsmatrisen65. Samma förfarande användes för att anpassa poängen från de två analyserna till en χ2-fördelning och kombinera de erhållna p-värdena som användes för hapFLK-testet. Den ojämna fördelningen av p-värdena omöjliggjorde dock en tillämpning av FDR-ramen och vi valde ut SNP:er inom de regioner som upptäcktes med hapFLK och som uppvisade p-värden <10-4. För dessa SNP:er använde vi VEP-annotationerna71 (Variant Effect Predictor) som genererades från Ensembl v74 får OARv3.1-genomannotationen för Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) och från get CHIR1.0-genomannotationen för Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/), som producerats av NCBI:s pipeline för eukaryotisk genomannotation. SNPs klassificerades som intergeniska, uppströms och nedströms (inklusive UTR) samt introniska och exoniska positioner. Skillnaderna mellan fördelningen av SNPs med FLK p-värden <10-4 och alla SNPs som användes för att upptäcka urvalssignaturer undersöktes med ett χ2-test.
Datatillgänglighet
Sekvenser och metadata som genererats för de 73 Ovis- och 72 Capra-proverna som användes i dessa analyser är offentligt tillgängliga. Allmän information och alla vcf-filer finns på webbplatsen Ensembl (http://projects.ensembl.org/nextgen/). Alla Fastq-filer, Bam-filer och de novo-sammansättningar av O. orientalis och C. aegagrus finns på European Nucleotide Archive (https://www.ebi.ac.uk/ena) under accessionskoden för Nextgen-projektet (PRJEB7436).