Sampling

Domestiske får (O. aries) og geder (C. hircus) blev udtaget i Iran (henholdsvis IROA- og IRCH-grupper) og Marokko (henholdsvis MOOA- og MOCH-grupper) for i alt 20 dyr pr. gruppe (Supplerende figur 6). Disse prøver blev indsamlet mellem januar 2008 og marts 2012 i den nordlige del af Marokko og mellem august 2011 og juli 2012 i det nordvestlige Iran inden for rammerne af det europæiske Nextgen-projekt (tilskudsaftale nr. 244356) i overensstemmelse med de etiske regler i Den Europæiske Unions direktiv 86/609/EØF. Øreklip blev indsamlet fra den distale del af øret på tilfældigt udvalgte dyr og blev straks opbevaret i 96 % ethanol i et døgn, inden de blev overført til silica-gelperler indtil DNA-ekstraktion.

De vilde arter asiatisk mouflon (O. orientalis) og bezoar ibex (C. aegagrus) blev udtaget som prøver i det nordvestlige Iran inden for domesticeringens vugge21,22. 13 væv fra asiatiske mufloner og 18 væv fra Bezoar ibex (henholdsvis IROO- og IRCA-grupper, supplerende figur 6) blev indsamlet enten fra dyr i fangenskab eller fra dyr, der for nylig var blevet jaget, og fra frosne prøver, der var tilgængelige i det iranske miljøministerium. Denne individbaserede prøveudtagningsmetode er designet til at minimere potentiel skævhed ved at undgå overrepræsentation af lokale effekter (f.eks, lokal indavl).

Supplerende data

Der blev desuden sammensat et verdensomspændende racepanel for får og geder (henholdsvis wpOA og wpCH). wpOA omfattede 20 prøver til gensekventering af hele genomet (WGS) med 12x dækning, der repræsenterer 20 forskellige verdensomspændende racer, som blev leveret af International Sheep Genome Consortium. wpCH bestod af 14 WGS-prøver, der blev sekventeret med 12x dækning, og som repræsenterer 9 europæiske individer, dvs, 2 franske Alpine- og 2 franske Saanen-prøver sekventeret af INRA, 5 italienske Saanen-prøver leveret af Parco Tecnologico Padano og 5 australske individer, dvs, 2 Boer-, 2 Rangeland- og 1 Cashmere-prøver leveret af CSIRO (Supplerende data 5).

Produktion af WGS-data

Genomisk DNA blev med succes ekstraheret fra alle vævsprøver ved hjælp af Macherey Nagel NucleoSpin 96 Tissue-kittet, idet producentens protokol blev tilpasset. Vævsprøverne blev udtaget i MN-blokke med firkantede brønde for at opnå 25 mg fragmenter pr. prøve. Der blev fremstillet tre og en halv MN kvadrat-96 blokke, og ekstraktionen blev udført ved hjælp af en Tecan Freedom Evo Liquid Handler i henhold til producentens protokol. Der blev udført et præ-lysetrin for at homogenisere prøverne med 180 µl T1-buffer og 25 µl proteinase K natten over ved 56 °C. For at justere bindingsbetingelserne blev der tilsat 200 µl BQ1-buffer, og prøvepladen blev inkuberet 1 time ved 70 °C; der blev efterfølgende tilsat 200 µl 100 % ethanol. Lysaterne blev overført til Nucleospin Tissue bindingsplade, og der blev sat vakuum (-0,2 bar, 5 min) for at fjerne gennemstrømningen. Der blev foretaget tre vasketrin med henholdsvis BW- og B5-buffer, og der blev igen sat vakuum på for at fjerne gennemstrømningen. Inden elueringen af genomisk DNA blev en silicamembran på Nucleospin Tissue binding plate silica-membran tørret under vakuum med mindst -0,6 bar i 10 minutter. Elueringstrinnet blev udført med 100 µl forvarmet BE-buffer (70 °C) og et centrifugeringstrin ved 3700 rcf i 5 min i 96-PCR-brønde. Genomisk DNA blev opbevaret ved 4 °C for at undgå frysetøning og testet for koncentration (som ng/μl) ved hjælp af Picogreen-metoden og ved hjælp af en Nanodrop.

Hele genomer blev resekventeret fra 500 ng genomisk DNA, der blev klippet til et område på 150-700 bp ved hjælp af Covaris® E210-instrumentet for hver prøve og anvendt til Illumina®-bibliotekspræparation ved hjælp af en halvautomatiseret protokol. Slutreparation, A-tailing og Illumina®-kompatible adaptorer (BioScientific) blev udført ved hjælp af SPRIWorks Library Preparation System og SPRI TE-instrumentet (Beckmann Coulter) i henhold til producentens protokol. Der blev anvendt en størrelsesudvælgelse på 300-600 bp for at genvinde de fleste fragmenter. DNA-fragmenterne blev amplificeret ved 12 PCR-cyklusser ved hjælp af Platinum Pfx Taq Polymerase Kit (Life® Technologies) og Illumina® adapterspecifikke primere. Bibliotekerne blev renset med 0,8x AMPure XP-perler (Beckmann Coulter) og analyseret med Agilent 2100 Bioanalyzer (Agilent® Technologies) og qPCR-kvantificering. Bibliotekerne blev sekventeret ved hjælp af 100 base-længde læsningskemi i parrede-end flowcelle på Illumina® HiSeq2000.

Illumina parrede-end læsninger for Ovis blev kortlagt til fårereferencegenomet (OAR v3.1, GenBank assemblage GCA_000298735.146), og for Capra til gedereferencegenomet (CHIR v1.0, GenBank assemblage GCA_000317765.147), ved hjælp af BWA-MEM48. Den BAM-fil, der blev produceret for hvert individ, blev sorteret ved hjælp af Picard SortSam og forbedret ved hjælp af sekventielt Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator og GATK IndelRealigner49 samt SAMtools calmd50.

Variantopdagelse blev udført ved hjælp af tre forskellige algoritmer: Samtools mpileup50, GATK UnifiedGenotyper51 og Freebayes52. Variantsteder blev identificeret uafhængigt af hinanden for hver af de seks grupper ved hjælp af algoritmernes multi-sample-tilstande for kald af algoritmer: (i) 162 prøver fra MOOA; (ii) 20 prøver fra IROA; (iii) 14 prøver fra IROO; (iv) 162 prøver fra MOCH; (v) 20 prøver fra IRCH; (vi) 19 prøver fra IRCA. For nogle grupper var WGS af flere individer tilgængelige som led i NextGen-projektet (se ovenfor). De prøver, der blev anvendt i denne undersøgelse, blev udvalgt med henblik på at opnå afbalancerede grupper på 20 individer, hvor det var muligt. For IRCA- og IROO-grupperne blev yderligere prøver tilgængelige på et senere tidspunkt og blev tilføjet til downstream-analyser. Dyr med lav tilpasning og calling-kvalitet blev fjernet for at opnå det endelige datasæt (Supplerende data 5).

Inden for hver gruppe var der to på hinanden følgende runder af filtrering af variantstedskvalitet. I filtreringsfase 1 blev opkald fra de tre algoritmer slået sammen, mens de opkald med den laveste tillid blev filtreret fra. Et variantsted blev godkendt, hvis det blev kaldt af mindst to forskellige kaldningsalgoritmer med en phred-variantkvalitet >30. En alternativ allel på et sted blev godkendt, hvis den blev kaldt af en af kalibreringsalgoritmerne, og genotypetallet var >0. Filtreringsfase 2 anvendte Variant Quality Score Recalibration by GATK. Først genererede vi et træningssæt af de variantsteder med den højeste tillid inden for gruppen, hvor (i) stedet kaldes af alle tre variant callere med phred-skaleret variantkvalitet >100, (ii) stedet er biallelisk, (iii) minor allelantallet er mindst 3, mens der kun tælles prøver med genotype phred-skaleret kvalitet >30. Træningssættet blev brugt til at opbygge en Gaussisk model ved hjælp af værktøjet GATK VariantRecalibrator ved hjælp af følgende variantannotationer fra UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Der blev anvendt en Gaussisk model på det fulde datasæt, hvorved der blev genereret en VQSLOD (log odds ratio for at være en sand variant). Steder blev filtreret, hvis VQSLOD <cutoff-værdi. Cutoff-værdien blev fastsat for hver gruppe på følgende måde: Minimum VQSLOD = {medianværdien af VQSLOD for varianter i træningssættet}-3 × {den absolutte medianafvigelse VQSLOD for varianter i træningssættet}. Forholdet mellem overgangs- og transversion-SNP’er tyder på, at det valgte cutoff-kriterium gav den bedste balance mellem selektivitet og følsomhed.

Der blev genereret SNP-kaldesæt for seks grupper af Ovis- og Capra-dyr (dvs. iranske og marokkanske husdyr og vilde dyr for hver slægt). Da de analyser, der blev udført i denne undersøgelse, krævede sammenligninger mellem grupper, skabte vi genotypekaldesæt på et konsistent sæt af SNP-steder for alle dyr fra enhver gruppe. For hver slægt fusionerede vi variantopkaldsstederne fra dens tre grupper og beholdt kun bialleliske positioner uden manglende data. Genotyper blev genkaldt på hvert biallelisk SNP-sted for alle individer af interesse af GATK UnifiedGenotyper ved hjælp af indstillingen GENOTYPE_GIVEN_ALLELES. På dette stadium blev listen over individer udvidet til at omfatte dyr, der tilhører verdensracepanelerne for får og geder (wpOA og wpCH), og yderligere vilde prøver, der blev tilgængelige på dette stadium (4 O. orientalis og 4 C. aegagrus). Genotyperne blev forbedret og faset inden for grupperne ved hjælp af Beagle 453 og derefter filtreret fra, hvis genotypesandsynligheden var mindre end 0,95. Endelig filtrerede vi steder, der var monomorfe på tværs af de forskellige delmængder af individer, der blev anvendt i denne undersøgelse (se nedenfor).

For at sammenligne de signaler af selektion, der blev påvist mellem Ovis og Capra, udførte vi en krydsaligning mellem de to referencegenomer. Først brugte vi den parvise alignment pipeline fra Ensembl release 69 code base54 til at alignere referencegenomerne af får (OARv3.1) og ged (CHIR1.0). Denne pipeline bruger LastZ55 til at tilpasse på DNA-niveau efterfulgt af en efterbehandling, hvor de tilpassede blokke kædes sammen i henhold til deres placering i begge genomer. LastZ-pipelinen til parvis tilpasning køres rutinemæssigt af Ensembl for alle understøttede arter, men geden er endnu ikke med i Ensembl. For at undgå bias i retning af en af arterne producerede vi to forskellige interspecifikke alignments. I den ene blev fåret anvendt som referencegenom og geden som ikke-reference, mens geden blev anvendt som referencegenom og fåret som ikke-reference i den anden. Forskellen er, at genomregioner fra referencearterne er tvunget til at kortlægge entydigt til enkelte loci hos ikke-referencearterne, mens genomregioner fra ikke-referencearterne har lov til at kortlægge til flere steder hos referencearterne. Vi har for kromosomsegmenter fra et referencegenom fået koordinaterne på ikke-referencegenomet. Endelig brugte vi for de SNP’er, der blev opdaget i en slægt, hele genomudligningen med referencegenomet for den anden slægt til at identificere de tilsvarende positioner (Supplerende tabel 6).

Genetisk struktur

For at beskrive den genetiske diversitet inden for grupper brugte vi VCFtools56 til at beregne sammenfattende statistikker over genetisk variation på de 73 individer for Ovis (dvs, 13 IROO, 20 IROA, 20 MOOA og 20 wpOA) og 72 individer for Capra (dvs. 18 IRCA, 20 IRCH, 20 MOCH og 14 wpCH). De statistiske data, der blev målt, var det samlede antal polymorfe varianter (S) for hele sættet af individer i hver slægt og inden for hver gruppe, den gennemsnitlige nukleotiddiversitet (π) inden for hver gruppe og indavlskoefficienten (F) for hvert individ. Inden for hver slægt blev forskellene mellem den vilde gruppe og hver domesticeret gruppe testet ved hjælp af en ensidig t-test for individuelle værdier for indavl og genetisk belastning og en tosidet Mann-Whitney-test for nukleotiddiversitet pr. lokalitet.

Den samlede divergens mellem de fire grupper inden for hver slægt (dvs, vilde, iranske og marokkanske husdyr og verdenspanel) blev estimeret ved hjælp af alle bialleliske SNP’er og den gennemsnitlige vægtede parvise Fst i henhold til Weir og Cockerham57 som implementeret i VCFtools56. Den genetiske struktur mellem grupperne blev vurderet med clusteringmetoden sNMF26 , efter at datasættet var blevet beskåret for at fjerne SNP’er med linkage disquilibrium (r²) større end 0,2 ved hjælp af VCFtools. Linkage disquilibrium (r²) blev beregnet mellem par af SNP’er inden for glidende vinduer med 50 SNP’er, idet en SNP pr. par blev fjernet tilfældigt, når r² var større end 0,2. For hver sNMF-analyse blev der udført fem kørsler med det samme antal klynger (K) med værdier for K fra 1 til 10. Vi brugte krydsentropikriteriet til at identificere den mest sandsynlige klyngeløsning, men alternative partitioner for forskellige antal K blev også undersøgt for at vurdere, hvordan individer blev fordelt mellem klynger.

For at adskille mellem fælles forfædre og blanding kørte vi TreeMix27 for i fællesskab at estimere populationsopdelinger og efterfølgende blandingsbegivenheder ved hjælp af det beskårede datasæt, der blev brugt til sNMF. Vi kørte TreeMix med -global-indstillingen for at forfine vores maximum likelihood-slutninger. Vi rodede TreeMix-træet med opdelingen mellem vilde og tamme individer. Blokstørrelsen for jackknifing var -k 500 SNP’er, hvilket omtrent svarer til 150 kb, hvilket overstiger de gennemsnitlige blokke af LD, der er fundet hos både får og geder. Vi genererede et Maximum Likelihood-træ uden migration og tilføjede derefter migrationsbegivenheder og undersøgte den inkrementelle ændring i den varians, der forklares af modellen, og restværdierne mellem individer. Målet var at opdage enhver potentiel høj restværdi eller migrationskant mellem vilde og tamme individer. For yderligere at undersøge den statistiske relevans af mulige blandingsvektorer, der er identificeret af TreeMix (Supplerende tabel 3), beregnede vi trepopulationstesten f328 som en formel test af genetisk introgression ved hjælp af qp3Pop-programmet i ADMIXTOOLS-pakken58 for hver kombination af grupper. For Capra blev wpCH-gruppen opdelt mellem australske racer, franske racer og italienske racer. Resultaterne er rapporteret i Supplerende data 2.

Demografisk inferens

For hver slægt udførte vi analyser af forfædres demografiske inferens ved hjælp af MSMC-modellen, der er implementeret i MSMC2-softwaren25. MSMC er baseret på den parvise sekventielt markoviske koalescens59; den bruger dog haplotyper af fasede genomsekvensdata som input. Til hver analyse anvendte vi to individer fra en gruppe, dvs. 4 haplotyper. Hver analyse blev gentaget for et andet tilfældigt sæt af to individer, dvs. en gentagelse af analysen pr. gruppe. Input- og outputfiler blev genereret og analyseret med de python-scripts, der leveres med MSMC-softwaren, og som findes på https://github.com/stschiff/msmc-tools. Analyseparametrene blev holdt som standardparametre, bortset fra mutationsraten, der blev sat til 2,5×10-8, og generationslængden blev sat til 2 år. For at estimere usikkerheden på tidsestimaterne varierede vi disse parametre (mutationshastighed på 2,5×10-8 og 1,0×10-8 i kombination med generationslængde på 2 og 4 år) og gav et groft estimat af domesticeringsperioden (se supplerende figur 2).

Genetisk belastning

Genetisk belastning blev estimeret på to måder. For det første ved at beregne den genetiske belastning for hvert individ som summen af skadelige fitness-effekter over alle proteinkodende genomiske positioner efter metoden af Librado et al.60. Kort sagt, som en proxy for evolutionær begrænsning brugte vi PhyloP-scorerne fra 46-vejs pattedyrs alignment (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/) som en proxy for evolutionær begrænsning. Ud fra denne tilpasning identificerede vi protein-kodende steder, der udviklede sig under funktionelle begrænsninger (PhyloP-score ≥1,5). For hvert Ovis- eller Capra-genom undersøgte vi derefter, om disse steder var muteret. Hvis det var tilfældet, summerede vi phyloP-scorerne over alle muterede steder, således at mutationer på steder med store begrænsninger bidrager forholdsmæssigt mere til det samlede belastningsestimat. Dette gav et belastningsestimat for hvert får/gedegenom. For at opnå en gennemsnitlig belastning pr. websted dividerede vi endelig med det samlede antal analyserede positioner. Det er værd at bemærke, at vi betingede os homozygote steder for at undgå at modellere dominanskoefficienten for mutationer på heterozygote steder (f.eks. recessiv, intermediær, dominerende). For det andet sammenlignede vi gen-for-gen den genetiske skadelige belastning i vilde og domesticerede Ovis-grupper ved at udføre en Wilcoxon-test, hvor den alternative hypotese er, at husdyrene har større belastning end vilde slægtninge. p-værdierne blev korrigeret for multiple test61 , og vi anvendte en tærskelværdi for justerede p-værdier < 0,05. Vi udførte en genontologiberigelsesanalyse på det sæt gener, der viste en signifikant stigning i genetisk belastning, ved hjælp af WebGestalt62,63. Da referencegenomerne er dårligt annoteret for gener, var vi afhængige af enkeltkopi-ortologer mellem vores art og menneske og mus. Gener fra X-kromosomet blev udelukket fra baggrundssættet. Vi udførte ikke denne analyse på Capra på grund af den højere indavl, der blev observeret i de vilde prøver.

Detektion af udvælgelsessignaturer

For at detektere signaturer af udvælgelse relateret til domesticering anvendte vi alle bialleliske SNP’er, der viste en mindre allelfrekvens på over 0,10 i mindst en af de tre testede grupper (dvs. de iranske og marokkanske domesticeringsgrupper og den vilde gruppe for hver slægt). Da vi forventede, at signaturer af selektion i forbindelse med domesticeringsprocessen ville være til stede i alle husdyr, vedtog vi følgende generelle strategi: Vi testede med hapFLK29 (se supplerende note 5 og supplerende figurer 9, 10 og 11) for hver slægt den vilde gruppe mod hver af de traditionelt forvaltede husdyrgrupper (dvs. iranske og marokkanske) og fokuserede på de fælles regioner, der formodentlig var under udvælgelse, og som blev påvist i begge tilfælde. Gruppens prøvestørrelser (n = 13-20) var forenelige med kravene til metoden29. Vi kontrollerede visuelt, om de konsistente signaturer af udvælgelse, der blev fundet med hapFLK, også var til stede i det tilsvarende verdenspanelsæt for hver slægt, men inkluderede ikke disse grupper i den statistiske test på grund af deres sammensætning af flere racer. Endelig søgte vi efter fælles signaler af selektion mellem Ovis og Capra ved hjælp af en stratificeret FDR-tilgang. Strategien er afbildet i supplerende fig. 4.

Vi udførte hapFLK-tests for at kontrastere den vilde gruppe til hver af de iranske og marokkanske grupper i hver slægt. Slægtskabsmatrixen blev beregnet ud fra Reynolds genetiske afstande64 mellem par af grupper ved hjælp af en tilfældig delmængde af en procent af varianterne. Det udledte populationstræ blev opbygget ved hjælp af naboforbindelsesalgoritmen. For hver SNP udførte vi hapFLK-testen, der inkorporerer haplotypiske oplysninger for at øge kraften til at påvise selektive fejninger. For hver testet SNP beregnede hapFLK-statistikken afvigelsen af haplotypiske frekvenser i forhold til den neutrale model, der er estimeret ved hjælp af slægtskabsmatrixen65 . For at udnytte information om linkage disquilibrium bruger hapFLK Scheet og Stephens’66 multipoint-modellen for multilocus-genotyper, som kan tilpasses til ikke-fasede data. En af de vigtigste anvendelser af denne model er at udføre faseskøn (fastPHASE-software66). I vores analyse blev modellen trænet på ufasede data, og derfor tager vores analyse højde for faseusikkerhed. Metoden blev anvendt til at omgruppere lokale haplotyper langs kromosomer i et bestemt antal klynger K sat til 25 ved hjælp af en skjult Markov-model.

For at identificere de fælles regioner, der formodentlig er under udvælgelse i de to traditionelt forvaltede husholdningsgrupper for hver slægt, kombinerede vi de to tidligere hapFLK-analyser. For hver analyse blev hapFLK-scorerne tilpasset til en χ2-fordeling for at opnå p-værdier (script findes på https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Resultaterne af de to kontraster mellem den vilde gruppe og hver af husdyrgrupperne blev kombineret ved hjælp af Stouffers metode67 for at opnå enkelt p-værdier for sammenligningen af vilde vs. husdyr. Endelig blev FDR-rammen68 anvendt på hele sættet af SNP’er for at konvertere de kombinerede p-værdier til q-værdier. SNP’er, der viste q-værdier < 10-2, blev bevaret og grupperet i genomiske regioner, når de var mindre end 50 kb fra hinanden.

For at undersøge, om selektionssignalet blev delt mellem Ovis og Capra, brugte vi først krydsjusteringen af de to referencegenomer til at identificere homologe segmenter. Vi anvendte derefter en stratificeret FDR-ramme69. Denne tilgang er baseret på det faktum, at der er en iboende stratificering i testene i betragtning af de forudgående oplysninger i de genetiske data69 , fordi den underliggende fordeling af sande alternative hypoteser kan være forskellig i henhold til den forskellige dynamik i forskellige genomiske regioner, hvilket fører til forskellige fordelinger af p-værdierne. Dette kræver, at man får FDR-justerede p-værdier (dvs. q-værdier) særskilt for de forskellige strata. Vi søgte efter konvergenser i hver slægt ved at adskille de regioner, der er homologe med dem, der er påvist i den anden slægt (kaldet det fælles stratum), og resten af genomet (kaldet det generelle stratum). Vi udvandt p-værdierne separat for hvert af de to definerede strata og beregnede derefter q-værdierne ved hjælp af FDR-rammen. Disse stratificerede q-værdier var de endelige størrelser, der blev taget i betragtning for statistisk signifikans (<10-2) med henblik på at påvise SNP’er under udvælgelse og sammenlægge dem i de tilsvarende genomiske regioner.

For at teste for konvergerende signaturer af udvælgelse, der differentierer vilde fra tamdyr i begge slægter, undersøgte vi forholdet mellem den signifikanstærskel, der anvendes på q-værdier (som vi lod variere fra 0,2 til 0,002) i den ene slægt, og den estimerede sandsynlighed for, at en SNP er udvalgt i det fælles stratum i den anden slægt ved hjælp af Storey et al.70 tilgang. En stigning i den udledte sandsynlighed med et fald i den tærskel, der anvendes på q-værdien (stigning i stringens), indikerer, at jo mere signifikant regionen er i én slægt, jo mere sandsynligt er det, at vi ville finde signifikante SNP’er i den anden slægt.

Vi filtrerede de selektionssignaler fra, som ikke var konsistente blandt de tre husholdningsgrupper. For hver detekteret region brugte vi de fasede haplotyper af individer, som blev grupperet ved hjælp af Neighbor-Joining-træer baseret på den procentvise identitet mellem sekvenserne. Kun regioner, der viste konsistente signaler, blev beholdt (Supplerende figur 5).

For at udlede, om de signaler for selektion, der blev påvist med hapFLK, indikerede afslapning af selektion eller positiv selektion hos husdyrene, estimerede vi forskellen i nukleotiddiversitet (π) på hver formodet region under selektion mellem de vilde og husdyrgrupperne. Vi udtrykte denne forskel som Δπ-indekset, der blev beregnet for hver genomisk region som forskellen mellem π beregnet for den vilde gruppe og gennemsnittet af π for de iranske og marokkanske tamme grupper, minus forskellen i π mellem disse to grupper beregnet over hele genomet:

En negativ værdi vil indikere, at nukleotiddiversiteten er lavere i den vilde gruppe sammenlignet med gennemsnittet af de to tamme grupper, og vil blive anset for at vise en afslapning af selektion i disse sidste grupper, diversificerende selektion i de tamme eller positiv selektion i de vilde grupper. Omvendt ville en positiv værdi være et tegn på retningsbestemt positiv eller stabiliserende selektion, der har fundet sted i de tamme grupper. Vi brugte også haplotypeklyngen til manuelt at verificere i hver region, om den påviste selektive fejning bekræftede de indikationer, der blev givet af Δπ-indekset.

Vi foretog funktionelle fortolkninger som følger. For hver region under udvælgelse betragtede vi regionen plus 50 kb på hver side for at identificere funktionelle roller og 5 kb opstrøms og nedstrøms for generne, og vi vurderede overlapningen mellem disse koordinater for at bevare de interessante gener. Endelig mente vi, at et gen var relateret til en given detekteret region, når regionens og genets positioner var overlappende. Vi vurderede derefter, hvilket gen der var det mest sandsynlige mål for udvælgelse ved at tage hensyn til det nærmeste gen til topsignalet, dvs. den position med den laveste q-værdi inden for regionen. Generne blev funktionelt annoteret ved hjælp af Uniprot (http://www.uniprot.org/) ved at tage hensyn til deres involvering i 30 underordnede termer (dvs. termernes direkte nedstammer) i kategorien “Biologisk proces” (dvs. GO:0008150). Vi hentede alle GO-termer, der svarer til hvert gen (Supplerende data 4) for 30 af de 33 kategorier, fordi vi ikke tog hensyn til tre termer, der ikke var involveret i pattedyrsfunktioner (dvs. GO:0006791 svovludnyttelse, GO:0006794 fosforudnyttelse, GO:0015976 kulstofudnyttelse). Vi udførte to χ2-tests for at sammenligne fordelingen af gener i GO-kategorierne, dvs. (i) gener under udvælgelse fra slægtsspecifikke regioner versus den fra homologe regioner, og (ii) alle gener under udvælgelse versus de 18.689 menneskelige gener, der er associeret til GO-termer i Swiss-Prot. For at fortolke genernes funktioner i en husdyrkontekst hentede vi også de tilgængelige oplysninger fra litteraturen om deres fænotypiske virkninger.

Slutteligt brugte vi FLK-statistikken for at finde de SNP’er inden for de tidligere påviste regioner, der var mest differentieret mellem vilde og tamme grupper, for at finde de SNP’er inden for de tidligere påviste regioner, der var mest differentieret mellem vilde og tamme grupper. Som for hapFLK repræsenterer den afvigelsen af allelfrekvenserne for enkeltmarkørerne i forhold til den neutrale model, der er estimeret ved hjælp af slægtskabsmatricen65. Den samme procedure blev anvendt til at tilpasse scorerne fra de to analyser til en χ2-fordeling og kombinere de opnåede p-værdier, som blev anvendt til hapFLK-testen. Den uensartede fordeling af p-værdierne udelukkede imidlertid anvendelsen af FDR-rammen, og vi udvalgte SNP’er inden for de regioner, der blev påvist med hapFLK, og som viste p-værdier <10-4. For disse SNP’er anvendte vi Variant Effect Predictor (VEP)-annotationerne71 , der blev genereret fra Ensembl v74-fårenes OARv3.1-genomannotation for Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) og fra gedens CHIR1.0-genomannotation produceret af NCBI’s eukaryotiske genomannotationspipeline for Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). SNP’er blev klassificeret som intergeniske, opstrøms og nedstrøms (herunder UTR’er) og introniske og exoniske positioner. Forskellene mellem fordelingen af SNP’er med FLK p-værdier <10-4 og alle de SNP’er, der blev brugt til at detektere selektionssignaturer, blev undersøgt med en χ2-test.

Datatilgængelighed

Sekvenser og metadata, der er genereret for de 73 Ovis- og 72 Capra-prøver, der blev brugt i disse analyser, er offentligt tilgængelige. Generelle oplysninger og alle vcf-filer kan findes på Ensembl-webstedet (http://projects.ensembl.org/nextgen/). Alle Fastq-filer, Bam-filer og de novo-samlinger af O. orientalis og C. aegagrus kan findes på European Nucleotide Archive (https://www.ebi.ac.uk/ena) under Nextgen-projektets accessionskode (PRJEB7436).