Semnături genomice convergente ale domesticirii la ovine și caprine | Nature Communications

Eșantionare

Oile domestice (O. aries) și caprinele (C. hircus) au fost eșantionate în Iran (grupurile IROA și, respectiv, IRCH) și Maroc (grupurile MOOA și, respectiv, MOCH) pentru un total de 20 de animale pe grup (Fig. Suplimentară 6). Aceste probe au fost colectate între ianuarie 2008 și martie 2012 în partea de nord a Marocului și între august 2011 și iulie 2012 în nord-vestul Iranului, în cadrul proiectului european Nextgen (Acordul de grant nr. 244356), în conformitate cu reglementările etice ale Directivei 86/609/CEE a Uniunii Europene. Clipurile auriculare au fost prelevate din partea distală a urechii unor animale alese la întâmplare și au fost depozitate imediat în etanol 96 % timp de o zi înainte de a fi transferate în perle de silicagel până la extragerea ADN-ului.

Speciile sălbatice mouflon asiatic (O. orientalis) și ibex Bezoar (C. aegagrus) au fost eșantionate în nord-vestul Iranului în cadrul leagănului de domesticire21,22. Treisprezece țesuturi de mufloni asiatici și 18 țesuturi de bezoar ibex (grupurile IROO și, respectiv, IRCA, Fig. suplimentară 6) au fost colectate fie de la animale captive, fie de la animale vânate recent, precum și de la eșantioane congelate disponibile la Departamentul iranian de mediu. Această abordare de eșantionare bazată pe indivizi este concepută pentru a minimiza potențialele distorsiuni prin evitarea suprareprezentării efectelor locale (de ex, consangvinizarea locală).

Date suplimentare

În mod suplimentar, a fost asamblat un panel de rase la nivel mondial pentru ovine și caprine (wpOA și, respectiv, wpCH). wpOA a inclus 20 de probe de re-secvențiere a genomului întreg (WGS) cu o acoperire de 12x, reprezentând 20 de rase diferite la nivel mondial, furnizate de International Sheep Genome Consortium. wpCH a constat din 14 probe WGS secvențiate cu o acoperire de 12x, reprezentând 9 indivizi europeni, de ex, 2 mostre Alpine franceze și 2 mostre Saanen franceze secvențiate de INRA, 5 mostre Saanen italiene furnizate de Parco Tecnologico Padano și 5 indivizi australieni, și anume 2 probe Boer, 2 probe Rangeland și 1 probă Cashmere furnizate de CSIRO (Date suplimentare 5).

Producția datelor WGS

DNA genomic a fost extras cu succes din toate probele de țesut cu ajutorul kitului Macherey Nagel NucleoSpin 96 Tissue, adaptând protocolul producătorului. Prelevarea țesuturilor a fost efectuată în blocuri MN cu godeuri pătrate pentru a obține 25 mg de fragmente per eșantion. Au fost pregătite trei blocuri pătrate MN-96 și jumătate, iar extracția a fost efectuată cu ajutorul unui manipulator de lichide Tecan Freedom Evo Liquid handler, conform protocolului producătorului. S-a efectuat o etapă de pre-liză pentru a omogeniza probele cu 180 µl de tampon T1 și 25 µl de proteinază K peste noapte la 56 °C. Pentru a ajusta condițiile de legare, s-au adăugat 200 µl de tampon BQ1 și placa de probe a fost incubată 1 h la 70 °C; ulterior, s-au adăugat 200 µl de etanol 100%. Lizații au fost transferați pe placa de legare Nucleospin Tissue și s-a aplicat un vid (-0,2 bar, 5 min) pentru a îndepărta fluxul. S-au efectuat trei etape de spălare cu tampoanele BW și, respectiv, B5, și s-a aplicat din nou un vid pentru a elimina fluxul. Înainte de eluția ADN-ului genomic, o membrană de silice Nucleospin Tissue binding plate a fost uscată în vid cu cel puțin -0,6 bar timp de 10 minute. Etapa de eluție a fost efectuată cu 100 µl de soluție tampon BE preîncălzită (70 °C) și o etapă de centrifugare la 3700 rcf timp de 5 minute în 96 de godeuri pentru PCR. ADN genomic a fost depozitat la 4 °C pentru a se evita înghețarea-dezghețarea și s-a testat concentrația (ca ng/μl) prin metoda Picogreen și cu ajutorul unui Nanodrop.

Genomii întregi au fost resecvențializați din 500 ng de ADN genomic care au fost secționați la un interval de 150-700 bp cu ajutorul instrumentului Covaris® E210 pentru fiecare probă și au fost utilizați pentru pregătirea bibliotecilor Illumina® printr-un protocol semi-automatizat. Repararea capetelor, ligatura cozii A și a adaptorilor compatibili Illumina® (BioScientific) au fost efectuate cu ajutorul sistemului de pregătire a bibliotecilor SPRIWorks și al instrumentului SPRI TE (Beckmann Coulter), conform protocolului producătorului. S-a aplicat o selecție de dimensiuni de 300-600 bp pentru a recupera majoritatea fragmentelor. Fragmentele de ADN au fost amplificate prin 12 cicluri de PCR cu ajutorul kitului de polimerază Platinum Pfx Taq Polymerase (Life® Technologies) și a primerilor specifici adaptorului Illumina®. Bibliotecile au fost purificate cu margele AMPure XP 0,8x (Beckmann Coulter) și analizate cu bioanalizatorul Agilent 2100 (Agilent® Technologies) și cuantificarea qPCR. Bibliotecile au fost secvențiate utilizând o chimie de citire cu lungimea de 100 de baze în celulă de flux împerecheată pe Illumina® HiSeq2000.

Lecturile împerecheate de la Illumina pentru Ovis au fost cartografiate la genomul de referință al oilor (OAR v3.1, ansamblul GenBank GCA_000298735.146), iar pentru Capra la genomul de referință al caprelor (CHIR v1.0, ansamblul GenBank GCA_000317765.147), utilizând BWA-MEM48. Fișierul BAM produs pentru fiecare individ a fost sortat folosind Picard SortSam și îmbunătățit folosind secvențial Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator și GATK IndelRealigner49, și SAMtools calmd50.

Descoperirea variantelor a fost efectuată folosind trei algoritmi diferiți: Samtools mpileup50, GATK UnifiedGenotyper51, și Freebayes52. Locurile variantelor au fost identificate în mod independent pentru fiecare dintre cele șase grupuri, utilizând modurile multi-eșantion ale algoritmilor de apelare: (i) 162 de eșantioane de la MOOA; (ii) 20 de eșantioane de la IROA; (iii) 14 eșantioane de la IROO; (iv) 162 de eșantioane de la MOCH; (v) 20 de eșantioane de la IRCH; (vi) 19 eșantioane de la IRCA. Pentru unele grupuri, WGS ale mai multor indivizi au fost disponibile ca parte a proiectului NextGen (a se vedea mai sus). Eșantioanele utilizate în prezentul studiu au fost selectate pentru a obține grupuri echilibrate de 20 de indivizi ori de câte ori a fost posibil. Pentru grupurile IRCA și IROO, eșantioane suplimentare au devenit disponibile într-o etapă ulterioară și au fost adăugate pentru analizele din aval. Animalele cu o calitate scăzută a alinierii și a apelării au fost eliminate pentru a obține setul final de date (Date suplimentare 5).

În cadrul fiecărui grup, au existat două runde succesive de filtrare a calității situsurilor variantelor. Etapa 1 de filtrare a fuzionat apelurile împreună de la cei trei algoritmi, filtrând în același timp apelurile cu cea mai scăzută încredere. Un sit de variantă a trecut dacă a fost apelat de cel puțin doi algoritmi de apelare diferiți cu o calitate a variantei phred >30. O alelă alternativă la un sit a fost acceptată dacă a fost apelată de oricare dintre algoritmii de apelare, iar numărul de genotipuri a fost >0. În etapa de filtrare 2 s-a utilizat varianta Variant Quality Score Recalibration de către GATK. În primul rând, am generat un set de antrenament al siturilor de variante cu cea mai mare încredere din cadrul grupului în care (i) situl este apelat de toți cei trei apelatori de variante cu o calitate a variantei cu scala phred >100, (ii) situl este bialelic, (iii) numărul de alele minore este de cel puțin 3, numărând în același timp numai eșantioanele cu o calitate a genotipului cu scala phred >30. Setul de antrenament a fost utilizat pentru a construi un model gaussian cu ajutorul instrumentului GATK VariantRecalibrator folosind următoarele adnotări de variante din UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Un model gaussian a fost aplicat la întregul set de date, generând un VQSLOD (log odds ratio de a fi o variantă adevărată). Site-urile au fost filtrate în cazul în care VQSLOD <valoarea limită. Valoarea limită a fost stabilită pentru fiecare grup după cum urmează: VQSLOD minim = {valoarea mediană a VQSLOD pentru variantele din setul de formare}-3 × {deviația absolută mediană VQSLOD a variantelor din setul de formare}. Raportul SNP de tranziție/transversie SNP a sugerat că criteriul de cutoff ales a oferit cel mai bun echilibru între selectivitate și sensibilitate.

S-au generat seturi de apeluri SNP pentru șase grupuri de animale Ovis și Capra (și anume, animale domestice iraniene și marocane, și animale sălbatice pentru fiecare gen). Deoarece analizele efectuate în acest studiu au necesitat comparații între grupuri, am creat seturi de apeluri de genotipuri la un set consistent de situri SNP pentru toate animalele din orice grup. Pentru fiecare gen, am fuzionat site-urile de apelare a variantelor din cele trei grupuri ale sale și am păstrat doar pozițiile bialelice fără date lipsă. Genotipurile au fost rechemate la fiecare situs SNP bialelic pentru toți indivizii de interes de către GATK UnifiedGenotyper, utilizând opțiunea GENOTYPE_GIVEN_ALLELES. În această etapă, lista indivizilor a fost extinsă pentru a include animalele aparținând panelurilor mondiale de rase de ovine și caprine (wpOA și wpCH) și eșantioane sălbatice suplimentare care au devenit disponibile în această etapă (4 O. orientalis și 4 C. aegagrus). Genotipurile au fost îmbunătățite și eșalonate în cadrul grupurilor de către Beagle 453, iar apoi au fost filtrate în cazul în care probabilitatea genotipului a fost mai mică de 0,95. În cele din urmă, am filtrat site-urile care au fost monomorfe în diferitele subseturi de indivizi utilizate în acest studiu (a se vedea mai jos).

Pentru a compara semnalele de selecție detectate între Ovis și Capra, am efectuat o aliniere încrucișată între cele două genomuri de referință. În primul rând, am utilizat pipeline-ul de aliniere pe perechi din baza de coduri Ensembl versiunea 6954 pentru a alinia genomurile de referință ale oilor (OARv3.1) și ale caprelor (CHIR1.0). Această conductă utilizează LastZ55 pentru a se alinia la nivel de ADN, urmată de o postprocesare în care blocurile aliniate sunt înlănțuite în funcție de locația lor în ambele genomuri. Conducta de aliniere pe perechi LastZ este rulată în mod obișnuit de Ensembl pentru toate speciile acceptate, dar capra nu este încă inclusă în Ensembl. Pentru a evita prejudecata în favoarea uneia dintre specii, am produs două alinieri interspecifice diferite. Unul a folosit oaia ca genom de referință și capra ca non-referință, în timp ce celălalt a folosit capra ca genom de referință și oaia ca non-referință. Diferența constă în faptul că regiunile genomice ale speciei de referință sunt forțate să se coreleze în mod unic cu un singur loci al speciei de nereferință, în timp ce regiunile genomice de nereferință pot fi corelate cu mai multe locații ale speciei de referință. Am obținut pentru segmentele de cromozomi dintr-un genom de referință coordonatele pe genomul de nereferință. În cele din urmă, pentru SNP-urile descoperite într-un gen, am utilizat alinierea întregului genom cu genomul de referință al celuilalt gen pentru a identifica pozițiile corespunzătoare (tabelul suplimentar 6).

Structura genetică

Pentru a descrie diversitatea genetică în cadrul grupurilor, am utilizat VCFtools56 pentru a calcula statisticile rezumative ale variației genetice pe cei 73 de indivizi pentru Ovis (i.e, 13 IROO, 20 IROA, 20 MOOA, și 20 wpOA) și 72 de indivizi pentru Capra (adică 18 IRCA, 20 IRCH, 20 MOCH, și 14 wpCH). Statisticile măsurate au fost numărul total de variante polimorfe (S) pentru întregul set de indivizi din fiecare gen și din cadrul fiecărui grup, diversitatea nucleotidică medie (π) în cadrul fiecărui grup și coeficientul de consangvinizare (F) pentru fiecare individ. În cadrul fiecărui gen, diferențele dintre grupul sălbatic și fiecare grup domestic au fost testate cu ajutorul unui test t unilateral pentru valorile de consangvinizare și de încărcătură genetică individuale și a unui test Mann-Whitney bilateral pentru diversitatea nucleotidelor per situs.

Diferența globală între cele patru grupuri din cadrul fiecărui gen (adică, sălbatic, domestici iranieni și marocani, și panelul mondial) a fost estimată folosind toate SNP-urile bialelice și media Fst ponderată pe perechi după Weir și Cockerham57 , așa cum a fost implementată în VCFtools56. Structura genetică între grupuri a fost evaluată cu ajutorul metodei de grupare sNMF26, după ce a fost curățat setul de date pentru a elimina SNP-urile cu dezechilibru de legătură (r²) mai mare de 0,2, utilizând VCFtools. Dezechilibrul de legătură (r²) a fost calculat între perechile de SNP în cadrul unor ferestre glisante de 50 de SNP, cu un SNP pe pereche eliminat în mod aleatoriu atunci când r² a fost mai mare de 0,2. Pentru fiecare analiză sNMF, au fost efectuate cinci rulări ale aceluiași număr de clustere (K) cu valori ale lui K de la 1 la 10. Am utilizat criteriul de entropie încrucișată pentru a identifica soluția de grupare cea mai probabilă, însă au fost explorate și partiții alternative pentru diferite numere de K pentru a evalua modul în care indivizii au fost împărțiți între grupe.

Pentru a face distincția între strămoșii comuni și amestecul, am rulat TreeMix27 pentru a estima în comun diviziunile populației și evenimentele ulterioare de amestec folosind setul de date curățat utilizat pentru sNMF. Am rulat TreeMix cu opțiunea -global pentru a rafina inferențele noastre de maximă verosimilitate. Am înrădăcinat arborele TreeMix cu diviziunea dintre indivizii sălbatici și cei domestici. Dimensiunea blocului pentru jackknifing a fost de -k 500 SNP, ceea ce corespunde aproximativ la 150 kb, depășind blocurile medii de LD găsite atât la ovine, cât și la caprine. Am generat un arbore de maximă verosimilitate fără migrație și apoi am adăugat evenimente de migrație și am examinat modificarea incrementală a varianței explicate de model și a valorilor reziduale între indivizi. Scopul a fost de a detecta orice potențială valoare reziduală ridicată sau margine de migrație între indivizii sălbatici și cei domestici. Pentru a explora în continuare relevanța statistică a posibililor vectori de amestec identificați de TreeMix (tabelul suplimentar 3), am calculat testul celor trei populații f328 ca un test formal de introgresie genetică, utilizând programul qp3Pop din suita ADMIXTOOLS58 pentru fiecare combinație de grupuri. Pentru Capra, grupul wpCH a fost împărțit între rasele australiene, rasele franceze și rasele italiene. Rezultatele sunt raportate în datele suplimentare 2.

Inferență demografică

Pentru fiecare gen, am efectuat analize de inferență demografică ancestrală folosind modelul MSMC implementat în software-ul MSMC225. MSMC se bazează pe coalescența markoviană secvențială pe perechi59; cu toate acestea, utilizează ca intrare haplotipuri ale datelor de secvență genomică fazată. Pentru fiecare analiză am utilizat doi indivizi dintr-un grup, deci 4 haplotipuri. Fiecare analiză a fost repetată pentru un alt set aleatoriu de doi indivizi, adică o replică a analizei pentru fiecare grup. Fișierele de intrare și de ieșire au fost generate și analizate cu ajutorul scripturilor python furnizate împreună cu software-ul MSMC și care se găsesc la https://github.com/stschiff/msmc-tools. Parametrii de analiză au fost păstrați ca fiind cei impliciți, cu excepția ratei de mutație, care a fost stabilită la 2,5×10-8, iar durata generației a fost stabilită la 2 ani. Pentru a estima incertitudinea privind estimările de timp, am variat acești parametri (rata de mutație de 2,5×10-8 și 1,0×10-8 în combinație cu lungimea generației de 2 și 4 ani) și am furnizat o estimare aproximativă a perioadei de domesticire (a se vedea figura suplimentară 2).

Carcina genetică

Carcina genetică a fost estimată în două moduri. În primul rând, prin calcularea încărcăturii genetice pentru fiecare individ ca sumă a efectelor de fitness dăunătoare pe toate pozițiile genomice codificatoare de proteine, urmând metoda lui Librado et al.60. Pe scurt, ca proxy pentru constrângerea evolutivă, am folosit scorurile PhyloP din alinierea mamiferelor în 46 de direcții (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Din această aliniere, am identificat locurile de codificare a proteinelor care evoluează sub constrângeri funcționale (scor phyloP ≥1,5). Pentru fiecare genom Ovis sau Capra, am investigat apoi dacă aceste situri au suferit mutații. În caz afirmativ, am însumat scorurile phyloP pentru toate siturile mutate, astfel încât mutațiile din siturile cu constrângeri ridicate să contribuie proporțional mai mult la estimarea sarcinii totale. Acest lucru a furnizat o estimare a încărcăturii pentru fiecare genom de oaie/capră. În cele din urmă, pentru a obține o sarcină medie pe site, am împărțit-o la numărul total de poziții analizate. Este demn de remarcat faptul că am condiționat pe site-urile homozigote pentru a evita modelarea coeficientului de dominanță al mutațiilor la site-urile heterozigote (de exemplu, recesiv, intermediar, dominant). În al doilea rând, am comparat, genă cu genă, încărcătura genetică deleterioasă în grupurile de Ovis sălbatic și domestic, efectuând un test Wilcoxon, ipoteza alternativă fiind că animalele domestice au o încărcătură mai mare decât rudele sălbatice. valorile p au fost corectate pentru teste multiple61 și am aplicat un prag de valori p ajustate < 0,05. Am efectuat o analiză de îmbogățire ontologică a genelor pe setul de gene care prezintă o creștere semnificativă a încărcăturii genetice folosind WebGestalt62,63. Deoarece genomurile de referință sunt slab adnotate pentru gene, ne-am bazat pe ortologii cu o singură copie între specia noastră și cea umană și de șoarece. Genele din cromozomul X au fost excluse din setul de referință. Nu am efectuat această analiză pe Capra din cauza consangvinizării mai mari observate în eșantioanele sălbatice.

Detecția semnăturilor de selecție

Pentru detectarea semnăturilor de selecție legate de domesticire, am folosit toate SNP-urile bialelice care prezintă o frecvență a alelei minore mai mare de 0,10 în cel puțin unul dintre cele trei grupuri testate (adică grupul iranian și grupurile domestice marocane, precum și grupul sălbatic pentru fiecare gen). Deoarece ne așteptam ca semnăturile de selecție legate de procesul de domesticire să fie prezente la toate animalele domestice, am adoptat următoarea strategie generală: am testat cu hapFLK29 (a se vedea nota suplimentară 5 și figurile suplimentare 9, 10 și 11) pentru fiecare gen grupul sălbatic față de fiecare dintre grupurile domestice gestionate în mod tradițional (adică grupul iranian și cel marocan) și ne-am concentrat asupra acelor regiuni comune presupuse a fi supuse selecției care au fost detectate în ambele cazuri. Dimensiunile eșantioanelor de grup (n = 13-20) au fost compatibile cu cerințele metodei29. Am verificat vizual dacă semnăturile consecvente ale selecției găsite cu hapFLK erau, de asemenea, prezente în setul de paneluri mondiale corespunzător fiecărui gen, dar nu am inclus aceste grupuri în testul statistic din cauza compoziției lor multiraselor. În cele din urmă, am căutat semnale comune de selecție între Ovis și Capra folosind o abordare FDR stratificată. Strategia este descrisă în Fig. suplimentară 4.

Am efectuat teste hapFLK pentru a contrasta grupul sălbatic cu fiecare dintre grupurile iraniene și marocane din fiecare gen. Matricea de rudenie a fost calculată din distanțele genetice ale lui Reynold64 între perechile de grupuri, folosind un subset aleatoriu de un procent din variante. Arborele de populație dedus a fost construit cu ajutorul algoritmului neighbor-joining. Pentru fiecare SNP, am efectuat testul hapFLK, care încorporează informații haplotipice pentru a crește puterea de detectare a măturilor selective. Pentru fiecare SNP testat, statistica hapFLK a calculat abaterea frecvențelor haplotipice în raport cu modelul neutru estimat prin matricea de rudenie65. Pentru a exploata informațiile privind dezechilibrul de legătură, hapFLK utilizează modelul multipunct al lui Scheet și Stephens’66 pentru genotipurile multilocus, care poate fi ajustat la datele nefazate. Una dintre principalele aplicații ale acestui model este de a efectua estimarea fazei (software-ul fastPHASE66). În analiza noastră, modelul a fost antrenat pe date nefazate și, prin urmare, analiza noastră ține cont de incertitudinea fazei. Metoda a fost utilizată pentru a regrupa haplotipurile locale de-a lungul cromozomilor într-un număr specificat de clustere K stabilit la 25, utilizând un model Markov ascuns.

Pentru a identifica regiunile comune presupuse a fi supuse selecției în cele două grupuri domestice gestionate în mod tradițional pentru fiecare gen, am combinat cele două analize hapFLK anterioare. Pentru fiecare analiză, scorurile hapFLK au fost ajustate la o distribuție χ2 pentru a obține valori p (script disponibil la https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Rezultatele celor două contraste dintre grupul sălbatic și fiecare dintre grupurile domestice au fost combinate folosind metoda lui Stouffer67 pentru a obține valori p unice pentru comparația dintre animalele sălbatice și cele domestice. În cele din urmă, cadrul FDR68 a fost aplicat la întregul set de SNP-uri pentru a converti valorile p combinate în valori q. SNP-urile care prezentau valori q < 10-2 au fost reținute și grupate în regiuni genomice atunci când se aflau la o distanță mai mică de 50 kb unul față de celălalt.

Pentru a investiga dacă semnalul de selecție a fost împărtășit între Ovis și Capra, am folosit mai întâi alinierea încrucișată a celor două genomuri de referință pentru a identifica segmentele omoloage. Apoi am aplicat un cadru FDR stratificat69. Această abordare se bazează pe faptul că există o stratificare inerentă a testelor, având în vedere informațiile anterioare din datele genetice69, deoarece distribuția subiacentă a ipotezelor alternative adevărate ar putea fi diferită în funcție de dinamica diferită a diferitelor regiuni genomice, ceea ce duce la distribuții diferite ale valorilor p. Acest lucru necesită obținerea valorilor p ajustate cu FDR (adică valorile q) separat pentru diferitele straturi. Am căutat convergențe în fiecare gen prin separarea regiunilor omoloage cu cele detectate în celălalt gen (denumite stratul comun) și restul genomului (denumit stratul general). Am extras valorile p separat pentru fiecare dintre cele două straturi definite și apoi am calculat valorile q prin intermediul cadrului FDR. Aceste valori q stratificate au fost cantitățile finale luate în considerare pentru semnificația statistică (<10-2) pentru a detecta SNP-urile supuse selecției și pentru a le unifica în regiunile genomice corespunzătoare.

Pentru a testa dacă există semnături convergente de selecție care diferențiază animalele sălbatice de cele domestice în ambele genuri, am examinat relația dintre pragul de semnificație aplicat valorilor q (pe care l-am făcut să varieze de la 0,2 la 0,002) într-un gen și probabilitatea estimată ca un SNP să fie selectat în stratul comun al celuilalt gen folosind Storey et al.70 abordare. O creștere a probabilității deduse odată cu scăderea pragului aplicat la valoarea q (creșterea rigurozității) indică faptul că, cu cât regiunea este mai semnificativă într-un gen, cu atât este mai probabil să găsim SNP-uri semnificative în celălalt gen.

Am filtrat semnalele de selecție care nu au fost consecvente între cele trei grupuri domestice. Pentru fiecare regiune detectată, am folosit haplotipurile fazate ale indivizilor care au fost grupate cu ajutorul arborilor Neighbor-Joining pe baza procentului de identitate dintre secvențe. Au fost păstrate doar regiunile care prezentau semnale consistente (Fig. Suplimentară 5).

Pentru a deduce dacă semnalele de selecție detectate cu hapFLK indicau relaxarea selecției sau o selecție pozitivă la domestici, am estimat diferența de diversitate nucleotidică (π) pe fiecare regiune putativă supusă selecției între grupurile sălbatice și domestice. Am exprimat această diferență sub forma indicelui Δπ, care a fost calculat pentru fiecare regiune genomică ca diferență între π calculat pentru grupul sălbatic și media lui π pentru grupurile domestice iranian și marocan, minus diferența de π dintre aceste două grupuri calculată pe întregul genom:

$$\Delta \pi = \left( {\pi _{{{\rm wilds}}} – \pi _{{\mathrm{\operatorname{{{iran-morocco}}}}}} \right)_{{{mathrm{{operatorname{genomic-region}}}} – \left( {\pi _{{{\rm wilds}} – \pi _{{\mathrm{\operatorname{iran-morocco}}}}} \right)_{{{mathrm{\operatorname{whole-genome}}}}$$$

O valoare negativă ar indica faptul că diversitatea nucleotidelor este mai mică în grupul sălbatic în comparație cu media celor două grupuri domestice și ar fi considerată ca arătând o relaxare a selecției în aceste ultime grupuri, o selecție diversificatoare în cele domestice sau o selecție pozitivă în cele sălbatice. Dimpotrivă, o valoare pozitivă ar indica o selecție direcțională pozitivă sau stabilizatoare care a avut loc în grupurile domestice. Am folosit, de asemenea, gruparea haplotipurilor pentru a verifica manual în fiecare regiune dacă măturarea selectivă detectată a confirmat indicațiile date de indicele Δπ.

Am efectuat interpretări funcționale după cum urmează. Pentru fiecare regiune supusă selecției, am considerat regiunea plus 50 kb de fiecare parte pentru a identifica rolurile funcționale și 5 kb în amonte și în aval de gene și am evaluat suprapunerea dintre aceste coordonate pentru a reține genele de interes. În cele din urmă, am considerat că o genă era legată de o anumită regiune detectată atunci când pozițiile regiunii și ale genei se suprapun. Am evaluat apoi ce genă a fost cel mai probabil vizată de selecție, luând în considerare cea mai apropiată genă de semnalul de vârf, adică poziția cu cea mai mică valoare q în cadrul regiunii. Genele au fost adnotate din punct de vedere funcțional cu ajutorul Uniprot (http://www.uniprot.org/), luând în considerare implicarea lor în 30 de termeni copii (adică descendenții direcți ai termenilor) din categoria „Biological Process” (adică GO:0008150). Am recuperat toți termenii GO care corespund fiecărei gene (Date suplimentare 4) pentru 30 din cele 33 de categorii, deoarece nu am luat în considerare trei termeni care nu erau implicați în funcțiile mamiferelor (de exemplu, GO:0006791 utilizarea sulfului, GO:0006794 utilizarea fosforului, GO:0015976 utilizarea carbonului). Am efectuat două teste χ2 pentru a compara distribuțiile genelor în categoriile GO, și anume: (i) genele supuse selecției din regiunile specifice genurilor față de cele din regiunile omoloage și (ii) toate genele supuse selecției față de cele 18.689 de gene umane asociate termenilor GO din Swiss-Prot. Pentru a interpreta funcțiile genelor într-un context zootehnic, am recuperat, de asemenea, informațiile disponibile în literatura de specialitate cu privire la efectele lor fenotipice.

În cele din urmă, pentru a găsi SNP-urile din cadrul regiunilor detectate anterior care au fost cele mai diferențiate între grupurile sălbatice și cele domestice, am folosit statistica FLK. În ceea ce privește hapFLK, aceasta reprezintă abaterea frecvențelor alelice ale unui singur marker în raport cu modelul neutru estimat prin matricea de rudenie65. Aceeași procedură a fost utilizată pentru a adapta scorurile din cele două analize la o distribuție χ2 și pentru a combina valorile p obținute, așa cum a fost utilizată pentru testul hapFLK. Cu toate acestea, distribuția neuniformă a valorilor p a împiedicat aplicarea cadrului FDR și am selectat SNP-urile din cadrul regiunilor detectate cu hapFLK care prezentau valori p <10-4. Pentru aceste SNP-uri am utilizat adnotările Variant Effect Predictor (VEP)71 care au fost generate din adnotarea Ensembl v74 a genomului ovin OARv3.1 pentru Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) și din adnotarea genomului caprin CHIR1.0 produsă de NCBI eukaryotic genome annotation pipeline pentru Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). SNP-urile au fost clasificate în poziții intergenice, în amonte și în aval (inclusiv UTR-uri) și în poziții intronice și exonice. Diferențele dintre distribuțiile SNP-urilor cu valori FLK p <10-4 și toate SNP-urile utilizate pentru detectarea semnăturilor de selecție au fost examinate cu un test χ2.

Disponibilitatea datelor

Secvențele și datele de metadate generate pentru cele 73 de probe Ovis și 72 de probe Capra utilizate în aceste analize sunt disponibile public. Informații generale și toate fișierele vcf pot fi găsite pe site-ul Ensembl (http://projects.ensembl.org/nextgen/). Toate fișierele Fastq, fișierele Bam și ansamblurile de novo ale O. orientalis și C. aegagrus pot fi găsite pe European Nucleotide Archive (https://www.ebi.ac.uk/ena) sub codul de acces al proiectului Nextgen (PRJEB7436).