Sampling

Domowe owce (O. aries) i kozy (C. hircus) zostały pobrane w Iranie (grupy IROA i IRCH, odpowiednio) i Maroku (grupy MOOA i MOCH, odpowiednio) w sumie 20 zwierząt na grupę (Supplementary Fig. 6). Próbki te zostały pobrane między styczniem 2008 a marcem 2012 w północnej części Maroka oraz między sierpniem 2011 a lipcem 2012 w północno-zachodnim Iranie, w ramach europejskiego projektu Nextgen (Grant Agreement no. 244356), zgodnie z przepisami etycznymi Dyrektywy Unii Europejskiej 86/609/EWG. Wycinki uszu pobrano z dystalnej części ucha losowo wybranych zwierząt i natychmiast przechowywano w 96% etanolu przez jeden dzień przed przeniesieniem do kulek żelu krzemionkowego do czasu ekstrakcji DNA.

Próbki dzikich gatunków muflona azjatyckiego (O. orientalis) i ibeksa bezoarowego (C. aegagrus) pobrano w północno-zachodnim Iranie w kolebce udomowienia21,22. Trzynaście tkanek muflonów azjatyckich i 18 tkanek ibeksa bezoarowego (odpowiednio grupy IROO i IRCA, ryc. 6 uzupełniająca) pobrano ze zwierząt trzymanych w niewoli lub niedawno upolowanych, a także z zamrożonych próbek dostępnych w irańskim Departamencie Środowiska. To indywidualne podejście do pobierania próbek ma na celu zminimalizowanie potencjalnej stronniczości poprzez unikanie nadreprezentacji efektów lokalnych (np, lokalnego chowu wsobnego).

Dodatkowe dane

Dodatkowo, ogólnoświatowy panel rasowy został zebrany dla owiec i kóz (odpowiednio wpOA i wpCH). wpOA zawierał 20 próbek ponownego sekwencjonowania całego genomu (WGS) przy 12-krotnym pokryciu reprezentujących 20 różnych światowych ras dostarczonych przez International Sheep Genome Consortium. wpCH składał się z 14 próbek WGS sekwencjonowanych przy 12-krotnym pokryciu reprezentujących 9 europejskich osobników, tj, 2 francuskie Alpine i 2 francuskie Saanen sekwencjonowane przez INRA, 5 włoskich Saanen dostarczonych przez Parco Tecnologico Padano, oraz 5 australijskich osobników, tj, 2 Boer, 2 Rangeland, i 1 Cashmere próbki dostarczone przez CSIRO (Supplementary Data 5).

Produkcja danych WGS

Genomowe DNA zostało z powodzeniem wyekstrahowane ze wszystkich próbek tkanek przy użyciu zestawu Macherey Nagel NucleoSpin 96 Tissue kit, dostosowując protokół producenta. Pobieranie próbek tkanek przeprowadzono w blokach z dołkami kwadratowymi MN, aby uzyskać 25 mg fragmentów na próbkę. Przygotowano trzy i pół MN kwadratowych 96 bloków, a ekstrakcję przeprowadzono przy użyciu urządzenia Tecan Freedom Evo Liquid handler zgodnie z protokołem producenta. Przeprowadzono etap wstępnej lizy w celu homogenizacji próbek z 180 µl buforu T1 i 25 µl proteinazy K przez noc w temperaturze 56 °C. Aby dostosować warunki wiązania, dodawano 200 µl buforu BQ1 i inkubowano płytkę z próbkami przez 1 h w 70 °C; następnie dodawano 200 µl 100% etanolu. Lizaty przenoszono na płytkę wiążącą Nucleospin Tissue i stosowano próżnię (-0,2 bar, 5 min) w celu usunięcia przesączu. Wykonano trzy etapy płukania odpowiednio buforami BW i B5, a następnie ponownie zastosowano próżnię w celu usunięcia wypływu. Przed elucją genomowego DNA, membrana krzemionkowa płytki wiążącej Nucleospin Tissue była suszona pod próżnią o ciśnieniu co najmniej -0,6 bara przez 10 min. Etap elucji przeprowadzono przy użyciu 100 µl wstępnie ogrzanego buforu BE (70 °C) i etapu wirowania przy 3700 rcf przez 5 min w studzienkach 96-PCR. Genomowe DNA przechowywano w temperaturze 4 °C, aby uniknąć zamrażania-rozmrażania i badano pod kątem stężenia (jako ng/μl) metodą Picogreen i przy użyciu Nanodrop.

Całe genomy poddano resekwencjonowaniu z 500 ng genomowego DNA, które ścinano do zakresu 150-700 bp przy użyciu instrumentu Covaris® E210 dla każdej próbki i wykorzystano do przygotowania biblioteki Illumina® za pomocą półautomatycznego protokołu. Naprawa końców, ogon A i ligacja adapterów kompatybilnych z Illumina® (BioScientific) zostały wykonane przy użyciu SPRIWorks Library Preparation System i urządzenia SPRI TE (Beckmann Coulter) zgodnie z protokołem producenta. Zastosowano selekcję wielkości 300-600 bp w celu odzyskania większości fragmentów. Fragmenty DNA były amplifikowane przez 12 cykli PCR przy użyciu Platinum Pfx Taq Polymerase Kit (Life® Technologies) i starterów specyficznych dla adapterów Illumina®. Biblioteki oczyszczano za pomocą kulek 0,8x AMPure XP (Beckmann Coulter) i analizowano za pomocą bioanalizatora Agilent 2100 (Agilent® Technologies) oraz kwantyzacji qPCR. Biblioteki sekwencjonowano przy użyciu chemii odczytów o długości 100 baz w paired-end flow cell na aparacie Illumina® HiSeq2000.

Illumina paired-end reads for Ovis were mapped to the sheep reference genome (OAR v3.1, GenBank assembly GCA_000298735.146), and for Capra to the goat reference genome (CHIR v1.0, GenBank assembly GCA_000317765.147), using BWA-MEM48. Plik BAM wytworzony dla każdego osobnika sortowano przy użyciu Picard SortSam i poprawiano przy użyciu sekwencyjnie Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator i GATK IndelRealigner49, oraz SAMtools calmd50.

Odkrywanie wariantów przeprowadzono przy użyciu trzech różnych algorytmów: Samtools mpileup50, GATK UnifiedGenotyper51, oraz Freebayes52. Miejsca zmienne były identyfikowane niezależnie dla każdej z sześciu grup, przy użyciu trybów wielopróbkowych algorytmów wywołujących: (i) 162 próbki z MOOA; (ii) 20 próbek z IROA; (iii) 14 próbek z IROO; (iv) 162 próbki z MOCH; (v) 20 próbek z IRCH; (vi) 19 próbek z IRCA. Dla niektórych grup, WGS większej liczby osobników były dostępne jako część projektu NextGen (patrz wyżej). Próbki użyte w obecnym badaniu zostały wybrane w celu uzyskania zrównoważonych grup 20 osobników, gdzie tylko było to możliwe. Dla grup IRCA i IROO, dodatkowe próbki stały się dostępne na późniejszym etapie i zostały dodane do dalszych analiz. Zwierzęta z niską jakością wyrównania i wywoływania zostały usunięte w celu uzyskania ostatecznego zestawu danych (Supplementary Data 5).

W obrębie każdej grupy, były dwie kolejne rundy filtrowania jakości miejsca wariantu. Etap filtrowania 1 łączył wywołania z trzech algorytmów, odfiltrowując jednocześnie wywołania o najniższej pewności. Miejsce wariantu przeszło, jeśli zostało wywołane przez co najmniej dwa różne algorytmy wywołujące z jakością wariantu phred >30. Alternatywny allel w danym miejscu przeszedł, jeśli został wywołany przez dowolny z algorytmów wywołujących, a liczba genotypów wynosiła >0. Etap filtrowania 2 wykorzystał rekalibrację wyników jakości wariantu przez GATK. Po pierwsze, wygenerowaliśmy zestaw treningowy miejsc o najwyższej wiarygodności w grupie, w której (i) miejsce jest wywoływane przez wszystkie trzy algorytmy wywołujące warianty z jakością wariantu >100, (ii) miejsce jest bialleliczne, (iii) liczba alleli mniejszościowych wynosi co najmniej 3 przy liczeniu tylko próbek z genotypem o jakości >30. Zbiór treningowy posłużył do zbudowania modelu gaussowskiego za pomocą narzędzia GATK VariantRecalibrator z wykorzystaniem następujących anotacji wariantów z UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Model gaussowski został zastosowany do pełnego zbioru danych, generując VQSLOD (log odds ratio of being a true variant). Miejsca były filtrowane, jeśli VQSLOD < wartość odcięcia. Wartość odcięcia została ustalona dla każdej grupy w następujący sposób: Minimalna wartość VQSLOD = {mediana wartości VQSLOD dla wariantów zestawu treningowego}-3 × {mediana bezwzględnego odchylenia VQSLOD wariantów zestawu treningowego}. Stosunek SNP przejściowych/transwersyjnych sugerował, że wybrane kryterium odcięcia dawało najlepszą równowagę między selektywnością i czułością.

Zestawy wywołań SNPs dla sześciu grup zwierząt Ovis i Capra zostały wygenerowane (tj. irańskie i marokańskie zwierzęta domowe oraz dziki dla każdego rodzaju). Ponieważ analizy przeprowadzone w tym badaniu wymagały porównań międzygrupowych, stworzyliśmy zestawy wywołań genotypów w spójnym zestawie miejsc SNP dla wszystkich zwierząt z każdej grupy. Dla każdego rodzaju, połączyliśmy miejsca wywołania wariantu z jego trzech grup i zachowaliśmy tylko pozycje bialleliczne bez brakujących danych. Genotypy zostały ponownie wywołane w każdym biallelicznym miejscu SNP dla wszystkich interesujących nas osobników przez GATK UnifiedGenotyper przy użyciu opcji GENOTYPE_GIVEN_ALLELES. Na tym etapie lista osobników została rozszerzona o zwierzęta należące do światowych paneli rasowych owiec i kóz (wpOA i wpCH) oraz dodatkowe dzikie próbki, które stały się dostępne na tym etapie (4 O. orientalis i 4 C. aegagrus). Genotypy były poprawiane i fazowane w obrębie grup przez Beagle 453, a następnie odfiltrowane, gdy prawdopodobieństwo genotypu było mniejsze niż 0,95. W końcu odfiltrowaliśmy miejsca, które były monomorficzne w różnych podzbiorach osobników użytych w tym badaniu (patrz poniżej).

W celu porównania sygnałów selekcji wykrytych między Ovis i Capra, wykonaliśmy wyrównanie krzyżowe między dwoma genomami referencyjnymi. Po pierwsze, wykorzystaliśmy potok wyrównywania parami z bazy kodowej Ensembl release 6954 do wyrównania genomów referencyjnych owcy (OARv3.1) i kozy (CHIR1.0). Potok ten wykorzystuje LastZ55 do wyrównania na poziomie DNA, po czym następuje postprocessing, w którym wyrównane bloki są łączone w łańcuchy zgodnie z ich lokalizacją w obu genomach. Potoczek wyrównywania par LastZ jest rutynowo uruchamiany przez Ensembl dla wszystkich wspieranych gatunków, ale koza nie jest jeszcze włączona do Ensembl. Aby uniknąć tendencyjności w kierunku jednego z gatunków, stworzyliśmy dwa różne wyrównania międzygatunkowe. Jedno z nich wykorzystywało owcę jako genom referencyjny, a kozę jako niereferencyjny, podczas gdy drugie wykorzystywało kozę jako genom referencyjny, a owcę jako niereferencyjny. Różnica polega na tym, że regiony genomowe gatunków referencyjnych są zmuszone do unikalnego mapowania do pojedynczych loci gatunków niereferencyjnych, podczas gdy regiony genomowe gatunków niereferencyjnych mogą mapować do wielu lokalizacji gatunków referencyjnych. Uzyskaliśmy dla segmentów chromosomów jednego genomu referencyjnego współrzędne na genomie niereferencyjnym. Wreszcie, dla SNPs odkrytych w jednym rodzaju, użyliśmy wyrównania całego genomu z genomem referencyjnym drugiego rodzaju, aby zidentyfikować odpowiadające im pozycje (Supplementary Table 6).

Struktura genetyczna

W celu opisania różnorodności genetycznej w obrębie grup, użyliśmy VCFtools56 do obliczenia statystyk zbiorczych zmienności genetycznej na 73 osobnikach dla Ovis (tj., 13 IROO, 20 IROA, 20 MOOA, i 20 wpOA) i 72 osobników dla Capra (tj. 18 IRCA, 20 IRCH, 20 MOCH, i 14 wpCH). Mierzono całkowitą liczbę wariantów polimorficznych (S) dla całego zbioru osobników w każdym rodzaju i w obrębie każdej grupy, uśrednioną różnorodność nukleotydową (π) w obrębie każdej grupy oraz współczynnik inbredu (F) dla każdego osobnika. W obrębie każdego rodzaju, różnice między grupą dziką a każdą grupą domową testowano za pomocą jednostronnego testu t dla indywidualnych wartości inbredu i obciążenia genetycznego oraz dwustronnego testu Manna-Whitneya dla różnorodności nukleotydów na miejsce.

Ogólna dywergencja między czterema grupami w obrębie każdego rodzaju (tj, dzikimi, udomowionymi irańskimi i marokańskimi oraz panelem światowym) została oszacowana przy użyciu wszystkich biallelicznych SNP i średniej ważonej parami Fst zgodnie z Weir i Cockerham57, jak zaimplementowano w VCFtools56. Strukturę genetyczną między grupami oceniono za pomocą metody grupowania sNMF26, po przycięciu zbioru danych w celu usunięcia SNPs z nierównowagą sprzężeń (r²) większą niż 0,2 za pomocą VCFtools. Nierównowaga pokrewieństwa (r²) została obliczona między parami SNPs w przesuwnych oknach 50 SNPs, z jednym SNP na parę losowo usuniętym, gdy r² było większe niż 0,2. Dla każdej analizy sNMF wykonano pięć przebiegów z taką samą liczbą klastrów (K) przy wartościach K od 1 do 10. Użyliśmy kryterium entropii krzyżowej, aby zidentyfikować najbardziej prawdopodobne rozwiązanie klastrowania, jednak alternatywne partycje dla różnych liczb K były również badane, aby ocenić, jak jednostki zostały podzielone między klastry.

Aby rozdzielić wspólne pochodzenie i domieszkowanie, uruchomiliśmy TreeMix27 , aby wspólnie oszacować podziały populacji i późniejsze zdarzenia domieszkowania przy użyciu przyciętego zestawu danych użytego do sNMF. Uruchomiliśmy TreeMix z opcją -global, aby udoskonalić nasze wnioskowanie na podstawie maksymalnego prawdopodobieństwa. Zakorzeniliśmy drzewo TreeMix z podziałem na osobniki dzikie i domowe. Rozmiar bloku dla jackknifingu wynosił -k 500 SNPs, co w przybliżeniu odpowiada 150 kb, przekraczając średnie bloki LD znalezione zarówno u owiec jak i kóz. Wygenerowaliśmy drzewo Maximum Likelihood bez migracji, a następnie dodaliśmy zdarzenia migracyjne i zbadaliśmy przyrostową zmianę wariancji wyjaśnionej przez model oraz wartości rezydualne między osobnikami. Celem było wykrycie wszelkich potencjalnych wysokich wartości rezydualnych lub krawędzi migracji pomiędzy dzikimi i domowymi osobnikami. Aby dokładniej zbadać znaczenie statystyczne możliwych wektorów domieszek zidentyfikowanych przez TreeMix (Tabela uzupełniająca 3), obliczyliśmy test trzech populacji f328 jako formalny test introgresji genetycznej, używając programu qp3Pop z pakietu ADMIXTOOLS58 dla każdej kombinacji grup. W przypadku Capry grupę wpCH podzielono między rasy australijskie, rasy francuskie i rasy włoskie. Wyniki przedstawiono w Supplementary Data 2.

Wnioskowanie demograficzne

Dla każdego rodzaju przeprowadziliśmy analizy wnioskowania demograficznego o przodkach przy użyciu modelu MSMC zaimplementowanego w oprogramowaniu MSMC225. MSMC opiera się na sekwencyjnej koalescencji parami59; jednak jako dane wejściowe wykorzystuje haplotypy z etapowych danych sekwencji genomu. Do każdej analizy używaliśmy dwóch osobników z jednej grupy, a więc 4 haplotypy. Każda analiza była powtarzana dla innego losowego zestawu dwóch osobników, czyli replikacja analizy na grupę. Pliki wejściowe i wyjściowe były generowane i analizowane za pomocą skryptów pythonowych dostarczonych z oprogramowaniem MSMC i znajdujących się na stronie https://github.com/stschiff/msmc-tools. Parametry analizy zostały zachowane jako domyślne, z wyjątkiem współczynnika mutacji, który został ustawiony na 2.5×10-8, a długość generacji została ustawiona na 2 lata. W celu oszacowania niepewności oszacowania czasu, zmieniliśmy te parametry (tempo mutacji 2.5×10-8 i 1.0×10-8 w połączeniu z długością pokolenia 2 i 4 lata) i dostarczyliśmy przybliżone oszacowanie okresu udomowienia (patrz Suplementary Fig. 2).

Obciążenie genetyczne

Obciążenie genetyczne zostało oszacowane na dwa sposoby. Po pierwsze, poprzez obliczenie obciążenia genetycznego dla każdego osobnika jako sumy szkodliwych efektów kondycyjnych na wszystkich pozycjach genomowych kodujących białka, zgodnie z metodą Librado i wsp.60. W skrócie, jako przybliżenie dla ograniczeń ewolucyjnych, wykorzystaliśmy wyniki PhyloP z 46-drożnego wyrównania ssaków (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Na podstawie tego wyrównania zidentyfikowaliśmy miejsca kodujące białka ewoluujące pod wpływem ograniczeń funkcjonalnych (PhyloP score ≥1,5). Dla każdego genomu Ovis lub Capra sprawdziliśmy, czy miejsca te były zmutowane. Jeśli tak, sumowaliśmy wyniki phyloP score dla wszystkich zmutowanych miejsc, tak aby mutacje w miejscach o wysokim ograniczeniu przyczyniały się proporcjonalnie bardziej do całkowitego oszacowania obciążenia. W ten sposób uzyskaliśmy oszacowanie obciążenia dla każdego genomu owcy/kóz. Ostatecznie, aby uzyskać średni ładunek na miejsce, podzieliliśmy go przez całkowitą liczbę analizowanych pozycji. Warto zauważyć, że uwarunkowaliśmy się na miejsca homozygotyczne, aby uniknąć modelowania współczynnika dominacji mutacji w miejscach heterozygotycznych (np. recesywne, pośrednie, dominujące). Po drugie, porównaliśmy gen po genie genetyczne obciążenie deleteryjne w dzikich i udomowionych grupach Ovis, wykonując test Wilcoxona, z alternatywną hipotezą, że zwierzęta domowe mają większe obciążenie niż dzicy krewni. wartości p zostały skorygowane o wielokrotne testowanie61 i zastosowaliśmy próg skorygowanych wartości p < 0,05. Przeprowadziliśmy analizę wzbogacenia ontologii genów na zestawie genów wykazujących znaczący wzrost obciążenia genetycznego przy użyciu WebGestalt62,63. Ponieważ genomy referencyjne są słabo anotowane dla genów, oparliśmy się na ortologach z pojedynczą kopią między naszym gatunkiem a człowiekiem i myszą. Geny z chromosomu X zostały wyłączone z zestawu tła. Nie przeprowadziliśmy tej analizy na Capra z powodu wyższego chowu wsobnego obserwowanego w dzikich próbkach.

Detekcja sygnatur selekcji

Do wykrywania sygnatur selekcji związanych z udomowieniem, użyliśmy wszystkich biallelicznych SNPs wykazujących mniejszą częstotliwość alleli większą niż 0,10 w co najmniej jednej z trzech badanych grup (tj. Irańskich i marokańskich grup domowych oraz dzikiej grupy dla każdego rodzaju). Ponieważ spodziewaliśmy się, że sygnatury selekcji związanej z procesem udomowienia będą obecne u wszystkich zwierząt domowych, przyjęliśmy następującą ogólną strategię: testowaliśmy za pomocą hapFLK29 (patrz Uzupełniająca uwaga 5 i Uzupełniające ryciny 9, 10 i 11) dla każdego rodzaju grupę dziką w stosunku do każdej z tradycyjnie zarządzanych grup domowych (tj. irańskiej i marokańskiej) i skupiliśmy się na tych wspólnych regionach rzekomo podlegających selekcji, które zostały wykryte w obu przypadkach. Wielkości próbek grupowych (n = 13-20) były zgodne z wymaganiami metody29. Wizualnie sprawdziliśmy, czy spójne sygnatury selekcji znalezione za pomocą hapFLK były również obecne w odpowiednim zestawie panelu światowego dla każdego rodzaju, ale nie włączyliśmy tych grup do testu statystycznego ze względu na ich wielorasowy skład. Wreszcie, szukaliśmy wspólnych sygnałów selekcji między Ovis i Capra przy użyciu podejścia stratyfikowanego FDR. Strategia ta jest przedstawiona w Supplementary Fig. 4.

Wykonaliśmy testy hapFLK dla kontrastu dzikiej grupy z każdą z irańskich i marokańskich grup w każdym rodzaju. Macierz pokrewieństwa obliczono na podstawie genetycznych odległości Reynolda64 między parami grup, używając losowego podzbioru jednego procenta wariantów. Drzewo populacji zostało zbudowane przy użyciu algorytmu neighbor-joining. Dla każdego SNP wykonaliśmy test hapFLK, który uwzględnia informacje haplotypowe, aby zwiększyć moc wykrywania selektywnych przemiatań. Dla każdego testowanego SNP statystyka hapFLK obliczała odchylenie częstotliwości haplotypowych w odniesieniu do neutralnego modelu oszacowanego przez macierz pokrewieństwa65. Aby wykorzystać informacje o nierównowadze sprzężeń, hapFLK wykorzystuje wielopunktowy model Scheeta i Stephensa66 dla genotypów wieloogniskowych, który może być dopasowany do danych niefazowanych. Jednym z głównych zastosowań tego modelu jest estymacja fazy (oprogramowanie fastPHASE66). W naszej analizie model był trenowany na niefazowanych danych, dlatego też nasza analiza uwzględnia niepewność fazy. Metoda została wykorzystana do przegrupowania lokalnych haplotypów wzdłuż chromosomów w określoną liczbę klastrów K ustaloną na 25, przy użyciu Ukrytego Modelu Markowa.

Aby zidentyfikować wspólne regiony przypuszczalnie podlegające selekcji w dwóch tradycyjnie zarządzanych grupach domowych dla każdego rodzaju, połączyliśmy dwie poprzednie analizy hapFLK. Dla każdej analizy wyniki hapFLK zostały dopasowane do rozkładu χ2 w celu uzyskania wartości p (skrypt dostępny na stronie https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Wyniki dwóch kontrastów między grupą dziką a każdą z grup domowych połączono przy użyciu metody Stouffera67 , aby uzyskać pojedyncze wartości p dla porównania zwierząt dzikich i domowych. W końcu, ramy FDR68 zostały zastosowane do całego zestawu SNPs, aby przekształcić połączone wartości p w wartości q. SNPs wykazujące wartości q < 10-2 zostały zachowane i pogrupowane w regiony genomowe, gdy były mniej niż 50 kb odległe od siebie.

Aby zbadać, czy sygnał selekcji był wspólny między Ovis i Capra, najpierw użyliśmy wyrównania krzyżowego dwóch genomów referencyjnych, aby zidentyfikować homologiczne segmenty. Następnie zastosowaliśmy warstwowe ramy FDR69. Podejście to opiera się na fakcie, że istnieje nieodłączne rozwarstwienie w testach, biorąc pod uwagę wcześniejsze informacje w danych genetycznych69, ponieważ podstawowy rozkład prawdziwych hipotez alternatywnych może być różny w zależności od różnej dynamiki różnych regionów genomowych, co prowadzi do różnych rozkładów wartości p. Wymaga to uzyskania wartości p skorygowanych o FDR (tj. wartości q) oddzielnie dla różnych warstw. Poszukiwaliśmy zbieżności w każdym rodzaju, oddzielając regiony homologiczne do tych wykrytych w innym rodzaju (określane jako warstwa wspólna) od reszty genomu (określanej jako warstwa ogólna). Wyodrębniliśmy wartości p oddzielnie dla każdej z dwóch zdefiniowanych warstw, a następnie obliczyliśmy wartości q za pomocą metody FDR. Te stratyfikowane wartości q były ostatecznymi wielkościami branymi pod uwagę dla istotności statystycznej (<10-2) w celu wykrycia SNP pod selekcją i połączenia ich z odpowiednimi regionami genomowymi.

Aby przetestować konwergentne sygnatury selekcji różnicujące dzikie zwierzęta od zwierząt domowych w obu rodzajach, zbadaliśmy związek między progiem istotności zastosowanym do wartości q (które uczyniliśmy zmiennymi od 0,2 do 0,002) w jednym rodzaju i szacowanym prawdopodobieństwem, że SNP jest wybrany we wspólnej warstwie drugiego rodzaju przy użyciu Storey i wsp.70 podejście. Wzrost wnioskowanego prawdopodobieństwa z obniżeniem progu zastosowanego do wartości q (wzrost rygorystyczności) wskazuje, że im bardziej znaczący region jest w jednym rodzaju, tym bardziej prawdopodobne jest, że znajdziemy znaczące SNP w drugim rodzaju.

Odfiltrowaliśmy sygnały selekcji, które nie były spójne wśród trzech grup domowych. Dla każdego wykrytego regionu użyliśmy fazowanych haplotypów osobników, które zostały zgrupowane przy użyciu drzew Neighbor-Joining opartych na procencie identyczności między sekwencjami. Tylko regiony wykazujące spójne sygnały zostały zachowane (Supplementary Fig. 5).

W celu wywnioskowania, czy sygnały selekcji wykryte za pomocą hapFLK wskazywały na rozluźnienie selekcji lub pozytywną selekcję w domestics, oszacowaliśmy różnicę w różnorodności nukleotydów (π) na każdym przypuszczalnym regionie podlegającym selekcji między grupami dzikimi i domowymi. Wyraziliśmy tę różnicę jako indeks Δπ, który został obliczony dla każdego regionu genomowego jako różnica między π obliczoną dla grupy dzikiej i średnią π dla irańskiej i marokańskiej grupy domowej, minus różnica π między tymi dwiema grupami obliczona dla całego genomu:

Wartość ujemna wskazywałaby, że różnorodność nukleotydów jest niższa w grupie dzikiej w porównaniu ze średnią dwóch grup krajowych i byłaby uważana za wskazującą na rozluźnienie selekcji w tych ostatnich grupach, selekcję różnicującą w grupach domowych lub selekcję pozytywną w grupach dzikich. I odwrotnie, wartość dodatnia wskazywałaby na kierunkową selekcję pozytywną lub stabilizującą, która wystąpiła w grupach krajowych. Wykorzystaliśmy również klasteryzację haplotypów do ręcznej weryfikacji w każdym regionie, czy wykryty przemiat selekcyjny potwierdził wskazania podane przez indeks Δπ.

Interpretacje funkcjonalne przeprowadziliśmy w następujący sposób. Dla każdego regionu poddawanego selekcji rozważaliśmy region plus 50 kb z każdej strony, aby zidentyfikować role funkcjonalne i 5 kb upstream i downstream genów i ocenialiśmy nakładanie się tych współrzędnych, aby zachować geny zainteresowania. Ostatecznie uznaliśmy, że gen był związany z danym wykrytym regionem, gdy pozycje regionu i genu pokrywały się. Następnie ocenialiśmy, który gen był najbardziej prawdopodobnym celem selekcji, biorąc pod uwagę gen najbliższy najwyższemu sygnałowi, tj. pozycję najniższej wartości q w regionie. Geny zostały zanotowane funkcjonalnie przy użyciu Uniprot (http://www.uniprot.org/), biorąc pod uwagę ich zaangażowanie w 30 terminów potomnych (tj. terminów bezpośrednio zstępujących) kategorii „Proces biologiczny” (tj. GO:0008150). Wyszukaliśmy wszystkie terminy GO odpowiadające każdemu genowi (Supplementary Data 4) dla 30 z 33 kategorii, ponieważ nie wzięliśmy pod uwagę trzech terminów, które nie były zaangażowane w funkcje ssaków (np. GO:0006791 wykorzystanie siarki, GO:0006794 wykorzystanie fosforu, GO:0015976 wykorzystanie węgla). Przeprowadziliśmy dwa testy χ2 w celu porównania rozkładu genów w kategoriach GO, tj. (i) genów podlegających selekcji z regionów specyficznych dla rodzaju w porównaniu do genów z regionów homologicznych oraz (ii) wszystkich genów podlegających selekcji w porównaniu do 18 689 ludzkich genów powiązanych z terminami GO w Swiss-Prot. W celu zinterpretowania funkcji genów w kontekście zwierząt gospodarskich, pobraliśmy również informacje dostępne w literaturze na temat ich efektów fenotypowych.

Wreszcie, aby znaleźć SNPs w obrębie wcześniej wykrytych regionów, które były najbardziej zróżnicowane pomiędzy grupami dzikimi i domowymi, użyliśmy statystyki FLK. Podobnie jak w przypadku hapFLK, reprezentuje ona odchylenie częstotliwości alleli pojedynczych markerów w odniesieniu do modelu neutralnego oszacowanego przez macierz pokrewieństwa65. Zastosowano tę samą procedurę, aby dopasować wyniki z dwóch analiz do rozkładu χ2 i połączyć uzyskane wartości p, jak to miało miejsce w przypadku testu hapFLK. Jednakże, niejednolity rozkład wartości p wykluczył zastosowanie ram FDR i wybraliśmy SNP w obrębie regionów wykrytych za pomocą hapFLK, wykazujących wartości p <10-4. Dla tych SNP użyliśmy anotacji Variant Effect Predictor (VEP)71 , które zostały wygenerowane z anotacji genomu owcy Ensembl v74 OARv3.1 dla Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) i z anotacji genomu kozy CHIR1.0 wyprodukowanej przez NCBI eukaryotyczny potok anotacji genomu dla Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). SNP zostały sklasyfikowane jako pozycje intergeniczne, upstreamowe i downstreamowe (w tym UTR) oraz intronowe i egzoniczne. Różnice między rozkładami SNP z FLK p-values <10-4 i wszystkimi SNP użytymi do wykrywania sygnatur selekcji zbadano testem χ2.

Dostępność danych

Sekwencje i dane metadanych wygenerowane dla 73 próbek Ovis i 72 próbek Capra użytych w tych analizach są publicznie dostępne. Informacje ogólne i wszystkie pliki vcf można znaleźć na stronie Ensembl (http://projects.ensembl.org/nextgen/). Wszystkie pliki Fastq, Bam i asemblacje de novo O. orientalis i C. aegagrus można znaleźć w European Nucleotide Archive (https://www.ebi.ac.uk/ena) pod kodem akcesyjnym projektu Nextgen (PRJEB7436).

.