Whole-genome sequencing of 128 camels across Asia reveals origin and migration of domestic Bactrian camels

Sample collection and whole-genome sequencing

W sumie 105 domowych wielbłądów baktryjskich z całej Azji, 19 dzikich wielbłądów baktryjskich z regionu Gobi-Altai w MG, jak również 4 dromadery z IRANU zostały zebrane do tego badania (Supplementary Fig. 1 i uzupełniająca tabela 1). Krajowe wielbłądy baktryjskie zostały wybrane tak, aby objąć jak najwięcej głównych regionów geograficznych, w tym 55 z Inner MG (IMG), Xinjiang (XJ) i Qinghai w Chinach, 28 z MG, 6 z KAZA, 10 z Rosji (RUS) i 6 z IRAN. Ponieważ w związku z szerokim wykorzystaniem domowych wielbłądów baktryjskich w Chinach i MG ukształtowała się różnorodność ras lokalnych, z regionów tych wybrano osiem różnych ras reprezentatywnych. Pozostałe domowe wielbłądy baktryjskie z Azji Środkowej żyły wokół Morza Kaspijskiego.

Po ekstrakcji DNA, poszczególne genomy zostały zsekwencjonowane do średniego 13-krotnego pokrycia (Supplementary Fig. 2 i Supplementary Table 2). Odczyty sekwencji zostały wyrównane do naszego poprzedniego złożenia genomu wielbłąda baktryjskiego29 w celu wywołania wariantów. Po rygorystycznym filtrowaniu (Supplementary Fig. 3), całkowicie zidentyfikowaliśmy 13,83 milionów polimorfizmów pojedynczego nukleotydu (SNPs) i 1,41 milionów małych indeli. Warto zauważyć, że stosunek przejścia do transwersji surowych SNP (2,29) był niższy niż ten opisany u dromaderów (2,31-2,34)31, ale został zwiększony do 2,44 po procedurach filtrowania, co sugeruje poprawę jakości zidentyfikowanych SNP. Funkcjonalna anotacja wariantów wykazała, że około 63,10% z nich było intergenicznych, 33,62% intronicznych, a 0,94% egzonicznych (Tabela 3). Zidentyfikowano 13,73 mln, 6,39 mln i 10,55 mln wariantów odpowiednio u domowych wielbłądów baktryjskich, dzikich wielbłądów baktryjskich i dromaderów. Chociaż dromadery były bardziej rozbieżne filogenetycznie od obu gatunków wielbłądów baktryjskich, domowe wielbłądy baktryjskie dzieliły więcej wariantów z dromaderami (66,73%) niż z dzikimi wielbłądami baktryjskimi (39,31%) (ryc. 4) ze względu na ogromną redukcję wariantów genetycznych zaobserwowaną u istniejących dzikich wielbłądów baktryjskich oraz przepływ genów między dromaderami i domowymi wielbłądami baktryjskimi. Wśród domowych wielbłądów baktryjskich było 12,68 milionów i 11,61 milionów wariantów zidentyfikowanych odpowiednio w populacjach wschodnioazjatyckich i środkowoazjatyckich (uzupełniające ryc. 4). Chociaż wielbłądy domowe pobrane z Azji Wschodniej były liczniejsze niż te z Azji Środkowej, liczba wariantów w każdej populacji nie wykazywała znaczącego zróżnicowania pomiędzy tymi dwoma obszarami (P-value = 0,77, dwuogonowy test t; Tabela uzupełniająca 4).

Różnorodność genetyczna i zróżnicowanie

Dla bardziej szczegółowego porównania różnorodności genetycznej pomiędzy różnymi populacjami, najpierw usunęliśmy 14 osobników wykazujących bliskie pokrewieństwo genetyczne z pozostałymi (Tabela uzupełniająca 5). Różnorodność nukleotydów π (Ryc. 1a) dromaderów (1,54 × 10-3) była znacząco wyższa niż wielbłądów baktryjskich ze wszystkich regionów geograficznych (0,88 × 10-3-1,11 × 10-3; Tabela uzupełniająca 6), co było sprzeczne z wcześniejszymi szacunkami heterozygotyczności opartymi na indywidualnych genomach10. Jedną z ważnych przyczyn może być praktyka hybrydyzacji między dromaderami a wielbłądami baktryjskimi w Azji Środkowej30. Wśród wielbłądów baktryjskich, dzika populacja wykazała najniższą π (0,88 × 10-3) w porównaniu ze wszystkimi populacjami krajowymi (Rys. 1a i Tabela uzupełniająca 6). Chociaż wynik ten narusza wiele przypadków, że dzikie zwierzęta zwykle mają wyższą różnorodność genetyczną niż ich domowe odpowiedniki, takie jak psy24, świnie27, i króliki32, to zdarzyłoby się to w przypadku zagrożonych dzikich zwierząt o wyjątkowo małej liczebności populacji, takich jak konie33. Ponadto populacje domowe żyjące w Azji Środkowej wykazywały generalnie większą różnorodność (1,03 × 10-3-1,11 × 10-3) niż te żyjące w Azji Wschodniej (0,95 × 10-3-1,02 × 10-3; ryc. 1a). Tendencję tę potwierdzał również współczynnik θ Wattersona (ryc. 5). Również w tym przypadku hybrydyzacja z dromaderami w Azji Centralnej może odpowiadać za większą różnorodność w tym regionie.

Potem zmierzyliśmy parami odległość genetyczną między populacjami wielbłądów za pomocą Fst Weira (ryc. 1b). Wynik był dobrze zgodny ze znaną filogenezą, która wskazywała, że dromadery miały najwyższy Fst z wielbłądami baktryjskimi (0,54-0,64), a dzikie wielbłądy baktryjskie miały drugi najwyższy Fst z wielbłądami domowymi (0,27-0,31). Zróżnicowanie wśród domowych wielbłądów baktryjskich było znacznie niższe, co wskazuje na ich niedawne pojedyncze pochodzenie. Co ciekawe, wśród domowych wielbłądów baktryjskich, te z IRANU wykazywały największą rozbieżność z pozostałymi (0.05-0.06). Aby potwierdzić zróżnicowanie populacji, skonstruowaliśmy drzewo NJ (ang. neighbor-joining) dla wszystkich osobników w oparciu o ich macierz identyczności parami (ang. identity-by-state – IBS) (Supplementary Fig. 6). Drzewo NJ również wspierało monofiletyczny klad wszystkich domowych wielbłądów baktryjskich, w ramach którego IRAN tworzył najgłębszą gałąź.

Potencjalnym problemem z oszacowaniami genomu populacji był błąd referencyjny, gdzie użycie pojedynczego genomu referencyjnego prowadziłoby do niskiej wydajności wywoływania wariantów dla osobników, które bardzo się od niego różniły34. Aby zbadać to uprzedzenie, porównaliśmy brakującą liczbę wariantów wśród trzech gatunków, biorąc głębokość sekwencjonowania jako kowariant (Tabela uzupełniająca 7). Analiza wariancji (ANOVA) wykazała, że chociaż domowe wielbłądy baktryjskie nie miały znaczącej różnicy z dzikimi (P-value = 0,50), rzeczywiście miały niższą liczbę brakujących wariantów niż dromadery (P-value = 4,38 × 10-3). Aby ocenić wpływ błędu systematycznego na nasze szacunki, ponownie obliczyliśmy różnorodność genetyczną i Fst przy użyciu tylko synonimicznych SNP (Supplementary Fig. 7), ponieważ sekwencje kodujące były bardziej prawdopodobne, aby być niezmienne dla różnych gatunków. W rezultacie, oszacowania oparte na synonimicznych SNP były w dobrej zgodzie z całym genomem dla wszystkich gatunków, co sugeruje, że błąd odniesienia miał tylko niewielki wpływ na nasze oszacowania genomu populacji.

Struktura populacji z domieszką

Aby ujawnić ogólną strukturę populacji z potencjalną domieszką, przycięliśmy SNP usuwając te z wysokim stopniem niezgodności powiązań i potencjalnymi efektami funkcjonalnymi. Analiza skalowania wielowymiarowego (MDS) oparta na przyciętym podzbiorze odtworzyła podobny wynik jak w przypadku pełnego zestawu (Rys. 2a i Supplementary Fig. 8). Zgodnie z oczekiwaniami, dromadery i dzikie wielbłądy baktryjskie mogły być oddzielone odpowiednio przez pierwszą i drugą współrzędną. Kiedy trzecia współrzędna została włączona do MDS, IRAN został oddzielony od wszystkich innych domowych wielbłądów baktryjskich (Ryc. 2a).

Aby oszacować różne proporcje przodków, przeprowadziliśmy analizę struktury populacji z Admixture35 zakładając K populacji przodków (Ryc. 2b). Procedura walidacji krzyżowej poparła, że K = 3 było optymalne (Supplementary Fig. 9), pokazując wyraźny podział między dromaderami, dzikimi wielbłądami baktryjskimi i domowymi wielbłądami baktryjskimi. Zaobserwowano ewidentną introgresję wielbłądów domowych Bactrian do irańskich dromaderów, przynajmniej u jednego dromadera. W związku z tym, rodowód dromadera był powszechny w środkowoazjatyckich populacjach wielbłądów baktryjskich, w tym w IRANIE, KAZIE i RUSIE, z udziałem szacowanym na 1-10%. Ponadto, zaobserwowaliśmy rodowód wielbłąda domowego u kilku dzikich osobników z udziałem 7-15%. Może to wynikać z polimorfizmu przodków, ale może też być spowodowane hybrydyzacją introgresywną, którą zaobserwowano w przypadku mtDNA36 i chromosomów Y15, i którą uznano za zagrożenie dla odrębności genetycznej dzikich gatunków. Co zaskakujące, dzikie wielbłądy nie wniosły prawie nic do rodowodu populacji domowych, nawet do populacji mongolskich, które dzielą bliskie siedliska z dzikimi wielbłądami. Chociaż większość domowych wielbłądów baktryjskich nie wykazywała zróżnicowania, gdy K rosło, IRAN był pierwszą populacją, która oddzieliła się z unikalnym rodowodem (K = 5; ryc. 2b).

Jako inną metodę badania struktury populacji z domieszkami, wywnioskowaliśmy drzewo populacji dla wielbłądów używając TreeMix37 (ryc. 2c). Gdy nie dodano żadnego śladu migracji, topologia drzewa ponownie wskazywała, że IRAN był pierwszą populacją wyodrębnioną wśród wszystkich krajowych wielbłądów baktryjskich (Supplementary Fig. 10). Zwiększenie liczby ścieżek migracji (m = 1-3) mogło znacznie poprawić dopasowanie modelu (Supplementary Fig. 11), który zidentyfikował przepływ genów z dromaderów do domowych wielbłądów baktryjskich w Azji Środkowej, w tym KAZA, RUS i IRAN z wagami migracji od 4% do 9% (Supplementary Table 8). Warto zwrócić uwagę, że choć szlak migracyjny wskazywał na koniec gałęzi dromaderów do IRANU, to w połowie gałęzi dromaderów do KAZA i RUS (ryc. 2c). Może to sugerować istnienie populacji-ducha związanej z dromaderami irańskimi, która wniosła wkład w rodowód KAZA i RUS. Dodatkowa ścieżka migracji (m = 4) mogła w dalszym ciągu poprawiać dopasowanie modelu, co wskazywało na migrację dromadera do rasy XJ (ryc. 2c). Chociaż TreeMix nie wykrył silnego sygnału migracji między dzikimi i domowymi wielbłądami baktryjskimi, pozostałości wykazały umiarkowaną domieszkę między dzikimi i wschodnioazjatyckimi rasami (Supplementary Fig. 11). Następnie użyliśmy mniej sparametryzowanego testu trzech i czterech populacji (F3/F4)38 do oceny statystycznej istotności tych zdarzeń domieszkowych. Ponownie, test F3 silnie potwierdził domieszkę dromaderów i wielbłądów baktryjskich w KAZA, RUS i IRAN (Tabela uzupełniająca 9). Bardziej czuły test F4 potwierdził istotnie większy zakres domieszki między dromaderami i wielbłądami baktryjskimi w Azji Środkowej w porównaniu z wielbłądami z Azji Wschodniej (tab. 10). Wśród tych ostatnich, wyższy stopień domieszki z dromaderami został wykryty w XJ niż w MG/IMG.

Dowód na pochodzenie z Azji Środkowej poprzez usunięcie introgresji

Wschodnia i Środkowa Azja były dwoma alternatywnymi regionami udomowienia wielbłądów baktryjskich w oparciu o dowody archeologiczne1,12,17, ale najbardziej prawdopodobny pozostał nierozstrzygnięty. Chociaż zaobserwowaliśmy największe zróżnicowanie genetyczne między populacją irańską a wszystkimi innymi krajowymi, istnienie domieszki między dromaderami i wielbłądami baktryjskimi w Azji Środkowej osłabiłoby poparcie dla wnioskowania o pochodzeniu. Aby zmniejszyć ten efekt, próbowaliśmy usunąć introgresowane segmenty dromaderów z genomu wielbłąda baktryjskiego za pomocą testu „BABA/ABBA „39. Pogrupowaliśmy populacje wschodnio- i środkowoazjatyckie, a następnie porównaliśmy wymianę alleli między tymi dwiema grupami z dromaderami (ryc. 3a). Ponieważ rodowód wielbłądów baktryjskich u jednego dromadera, jak również rodowód wielbłądów domowych u trzech dzikich (ryc. 2b), byłyby czynnikami zakłócającymi, usunęliśmy te cztery osobniki z analizy. Użyliśmy statystyki fd, odpornej wersji D Pattersona, do zlokalizowania segmentów introgresowanych40 i zastosowaliśmy ścisły poziom istotności Z-score = 2 przy użyciu procedury Jackknife (Supplementary Fig. 12). W sumie 21 153 nienakładających się 100 kb segmentów w całym genomie, było znacznie więcej segmentów z przypuszczalnymi sygnałami introgresji w populacjach Azji Środkowej (11 711, Z-score > 2) niż w populacjach Azji Wschodniej (3891, Z-score < -2), jak oczekiwano. Następnie przeprowadziliśmy analizę Admixture w oparciu o pozostałe segmenty i potwierdziliśmy, że introgresja dromaderów została skutecznie ograniczona (Supplementary Fig. 13). Ponowne obliczenie parami Fst po usunięciu introgresji nadal wskazywało, że IRAN był najbardziej zróżnicowany (0,04-0,06) wśród wszystkich populacji krajowych (ryc. 3b). Aby uzyskać lepszy wgląd w filogenezę populacji, zrekonstruowaliśmy drzewo NJ w oparciu o Fst i wykonaliśmy test bootstrapowy (Ryc. 3b). Potwierdził on, że IRAN był pierwszą populacją, która wyodrębniła się spośród wszystkich krajowych populacji baktrianów, a następnie KAZA i RUS. Bayesowska binarna analiza Markov Chain Monte Carlo (MCMC) oparta na filogenezie silnie wspierała Azję Środkową jako obszar przodków domowych wielbłądów baktryjskich (prawdopodobieństwo = 99,78%) i późniejszą drogę dyspersji z Azji Środkowej do Wschodniej (Supplementary Fig. 14).

Jako niezależny dowód, zrekonstruowaliśmy również drzewo maksymalnego prawdopodobieństwa pełnej długości mtDNA w oparciu o 128 próbek, które zsekwencjonowaliśmy w tym badaniu, a także 39 dodatkowych próbek dostępnych z Genbank (ryc. 3c i tabela uzupełniająca 11). Introgresja mtDNA mogła być łatwo zidentyfikowana i wykluczona z drzewa. Na przykład, dwa nowo zsekwencjonowane mtDNA z KAZA i RUS zostały połączone w klastry z dromaderami. W obrębie kladu domowych wielbłądów baktryjskich, mimo że większość wielbłądów z różnych regionów geograficznych była mieszana, dwa mtDNA z IRANu tworzyły najbardziej podstawowe gałęzie populacji domowych (ryc. 3c). Bayesowska binarna analiza MCMC ponownie poparła środkowoazjatyckie pochodzenie domowych wielbłądów baktryjskich (prawdopodobieństwo = 76,43%).

Historia demograficzna wielbłądów baktryjskich

Wykonaliśmy kilka parametrycznych analiz modelowania, aby wywnioskować dynamikę demograficzną wielbłądów w historii. Zgodnie z poprzednimi badaniami10, długoterminowe trajektorie wielbłądów baktryjskich oparte na modelu PSMC (pairwise sequentially Markovian coalescent)41 ujawniły ogromny spadek efektywnej wielkości populacji przodków wielbłądów wcześniej niż milion lat temu (Supplementary Fig. 15). Chociaż długoterminowe trajektorie dzikich i domowych wielbłądów baktryjskich były generalnie podobne, oczywiste było, że różnią się one od siebie już 0,4 miliona lat temu, co wyklucza te pierwsze jako bezpośrednich przodków tych drugich, zgodnie z poprzednimi analizami mtDNA9,14.

Aby zbadać czas dywergencji między populacjami wielbłądów bardziej szczegółowo, użyliśmy uogólnionego próbnika koalescencji filogenetycznej (G-PhoCS)42. Biorąc pod uwagę filogenezę populacji wielbłądów, G-PhoCS mógł oszacować skalowany mutacjami rozmiar populacji i czas dywergencji populacji w oparciu o niepowiązane neutralne loci w poszczególnych genomach z każdej populacji. Aby zredukować złożoność modelu, uwzględniliśmy tylko dromadery, dzikie wielbłądy baktryjskie i trzy reprezentatywne populacje (IRAN, KAZA i MG) domowych wielbłądów baktryjskich (Supplementary Fig. 16 i Supplementary Table 12). Zgodnie z Rys. 3b, IRAN i KAZA były pierwszymi dwiema populacjami środkowoazjatyckimi, które się oddzieliły, a podział MG może wskazywać na dyspersję ze środkowej do wschodniej Azji. Wiek skalibrowano, przyjmując dywergencję wielbłąda baktryjskiego i dromadera na 5,73 mln lat, zgodnie z bazą danych Timetree43. Gdy nie uwzględniono pasma migracji, łatwo osiągnięto zbieżność wszystkich oszacowań parametrów (Supplementary Fig. 17 i Supplementary Table 13). Podobnie jak w przypadku wyników PSMC, efektywna wielkość populacji była generalnie zmniejszona od populacji przodków do współczesnych (ryc. 4). Czas dywergencji między dzikimi i domowymi wielbłądami baktryjskimi oszacowano na 0,43 mln lat temu (95% przedział ufności: 0,13-0,73 Mya; ryc. 4), co było nieco później niż na podstawie mtDNA (0,714 lub 1,1 Mya9). Wśród populacji rodzimych, IRAN został oddzielony od innych około 4,45 tys. lat temu (95% CI: 0,07-17,6 Kya), a następnie populacje środkowo- i wschodnioazjatyckie zostały oddzielone około 2,40 tys. lat temu (95% CI: 0,01-7,84 Kya; ryc. 4).

Aby umożliwić przepływ genów, próbowaliśmy również wprowadzić pasma migracji z dromaderów do populacji wielbłądów baktryjskich (Supplementary Fig. 16 i Supplementary Table 12). Oszacowania mogły być zbieżne tylko wtedy, gdy wprowadzono pasmo migracji z irańskich dromaderów do IRAN i pasmo migracji z populacji widmowej do KAZA (Uzupełniające Rys. 18 i Uzupełniająca Tabela 13). Chociaż czas dywergencji między dzikimi i domowymi wielbłądami baktryjskimi nie uległ zmianie w modelu migracji (0,46 Mya, 95% CI: 0,24-0,71 Mya), to czas pierwszej dywergencji populacji domowych (0,19 Mya, 95% CI: 0,08-0,31 Mya) stał się nierealistyczny, ponieważ znacznie wykraczał poza znaną historię udomowienia zwierząt gospodarskich (11,5 Kya44). Poza tym całkowity wskaźnik migracji oszacowano jedynie na 0,27% i 0,16% dla pasma migracji do IRANU i KAZA, znacznie niższy niż ten oszacowany za pomocą Admixture (1-10%). Możliwym powodem słabego oszacowania może być bardziej złożona historia domieszek niż ciągły model migracji ze stałymi wskaźnikami założony przez G-PhoCS.