Celogenomové sekvenování 128 velbloudů napříč Asií odhalilo původ a migraci domácích velbloudů baktrijských

Sběr vzorků a celogenomové sekvenování

Celkem 105 domácích velbloudů baktrijských napříč Asií, 19 divokých velbloudů baktrijských z oblasti Gobi-Altaj v MG a také 4 dromedáři z IRÁNU byli shromážděni pro tuto studii (Supplementary Fig. 1 a doplňková tabulka 1). Domácí velbloudi baktrijští byli vybráni tak, aby pokryli co nejvíce hlavních zeměpisných oblastí, včetně 55 velbloudů z Vnitřní MG (IMG), Sin-ťiangu (XJ) a Čching-chaj v Číně, 28 z MG, 6 z KAZA, 10 z Ruska (RUS) a 6 z Íránu. Vzhledem k tomu, že v důsledku širokého využívání domácích velbloudů baktrijských v Číně a MG vznikla různorodá místní plemena, bylo z těchto regionů vybráno osm různých reprezentativních plemen. Ostatní domácí baktrijští velbloudi ze Střední Asie žili v okolí Kaspického moře.

Po extrakci DNA byly jednotlivé genomy sekvenovány s průměrným 13× pokrytím (doplňkový obr. 2 a doplňková tabulka 2). Čtení sekvencí byla zarovnána k naší předchozí sestavě genomu velblouda baktrijského29 pro volání variant. Po přísném filtrování (doplňkový obr. 3) jsme celkem identifikovali 13,83 milionu jednonukleotidových polymorfismů (SNP) a 1,41 milionu malých indelů. Pozoruhodné je, že poměr přechodů k transverzím u surových SNP (2,29) byl nižší než poměr uváděný u dromedárů (2,31-2,34)31 , ale po filtračních procedurách se zvýšil na 2,44, což naznačuje zlepšení kvality identifikovaných SNP. Funkční anotace variant ukázala, že přibližně 63,10 % z nich bylo intergenických, 33,62 % intronických a 0,94 % exonických (doplňková tabulka 3). U domácích velbloudů baktrijských, divokých velbloudů baktrijských a dromedárů bylo identifikováno 13,73 milionu, 6,39 milionu a 10,55 milionu variant. Ačkoli se dromedáři ve fylogenezi více lišili od obou druhů velbloudů baktrijských, domácí velbloudi baktrijští sdíleli více variant s dromedáry (66,73 %) než s divokými velbloudy baktrijskými (39,31 %) (doplňkový obr. 4), což bylo způsobeno obrovskou redukcí genetických variant pozorovanou u existujících divokých velbloudů baktrijských a tokem genů mezi dromedáry a domácími velbloudy baktrijskými. Mezi domácími velbloudy baktrijskými bylo ve východoasijské a středoasijské populaci identifikováno 12,68 milionu a 11,61 milionu variant (doplňkový obr. 4). Přestože velbloudů domácích odebraných z východní Asie bylo více než velbloudů ze střední Asie, počet variant privátních pro každou populaci nevykazoval mezi oběma oblastmi žádnou významnou odchylku (P-hodnota = 0,77, dvouvýběrový t-test; Doplňková tabulka 4).

Genetická diverzita a diferenciace

Pro podrobnější srovnání genetické diverzity mezi různými populacemi jsme nejprve odstranili 14 jedinců vykazujících úzkou genetickou příbuznost se zbývajícími (Doplňková tabulka 5). Párová nukleotidová diverzita π (obr. 1a) dromedárů (1,54 × 10-3) byla výrazně vyšší než diverzita velbloudů baktrijských ze všech geografických oblastí (0,88 × 10-3-1,11 × 10-3; doplňková tabulka 6), což bylo v rozporu s předchozími odhady heterozygotnosti na základě jednotlivých genomů10. Jedním z důležitých důvodů by mohla být praxe hybridizace mezi dromedáry a velbloudy baktrijskými ve střední Asii30. Mezi velbloudy baktrijskými vykazovala divoká populace nejnižší π (0,88 × 10-3) ve srovnání se všemi domácími populacemi (obr. 1a a doplňková tab. 6). Ačkoli tento výsledek porušil mnoho případů, že divoká zvířata mají obvykle vyšší genetickou diverzitu než jejich domácí protějšky, jako jsou psi24, prasata27 a králíci32 , stalo by se tak v případě ohrožených divokých zvířat s extrémně malou velikostí populace, jako jsou koně33. Domácí populace žijící ve střední Asii navíc obecně vykazovaly vyšší diverzitu (1,03 × 10-3-1,11 × 10-3) než populace žijící ve východní Asii (0,95 × 10-3-1,02 × 10-3; obr. 1a). Tuto tendenci potvrdil i Wattersonův θ (doplňkový obr. 5). Za vyšší diverzitu v tomto regionu může opět hybridizace s dromedáry ve Střední Asii.

Poté jsme měřili párovou genetickou vzdálenost mezi populacemi velbloudů pomocí Weirova Fst (obr. 1b). Výsledek byl v dobré shodě se známou fylogenezí, která ukázala, že dromedáři měli nejvyšší Fst s velbloudy baktrijskými (0,54-0,64) a divocí velbloudi baktrijští měli druhou nejvyšší Fst s domácími (0,27-0,31). Diferenciace mezi domácími baktrijskými velbloudy byla mnohem nižší, což odpovídá jejich nedávnému jednotnému původu. Zajímavé je, že mezi domácími baktrijskými velbloudy vykazovali velbloudi z IRÁNU největší divergenci se všemi ostatními (0,05-0,06). Pro ověření populační diferenciace jsme pro všechny jedince sestrojili strom sousedského spojení (NJ) na základě jejich párové matice identity podle stavu (IBS) (doplňkový obr. 6). Strom NJ rovněž podpořil monofyletický klad všech domácích velbloudů baktrijských, v jehož rámci IRAN tvořil nejhlubší větev.

Potenciálním problémem populačních genomických odhadů bylo referenční zkreslení, kdy použití jediného referenčního genomu by vedlo k nízké účinnosti při volání variant u jedinců, kteří se od něj velmi liší34. Abychom toto zkreslení prozkoumali, porovnali jsme počet chybějících variant mezi třemi druhy, přičemž jsme jako kovariátu vzali hloubku sekvenování (doplňková tabulka 7). Analýza rozptylu (ANOVA) ukázala, že ačkoli se domácí velbloudi baktrijští významně nelišili od divokých (P-hodnota = 0,50), měli skutečně nižší počet chybějících variant než dromedáři (P-hodnota = 4,38 × 10-3). Abychom zhodnotili dopad zkreslení na naše odhady, přepočítali jsme genetickou diverzitu a Fst pouze se synonymními SNP (doplňkový obr. 7), protože u kódujících sekvencí je větší pravděpodobnost, že budou napříč druhy neměnné. Výsledkem bylo, že odhady založené na synonymních SNP byly u všech druhů v dobré shodě s celým genomem, což naznačuje, že referenční zkreslení mělo na naše odhady populační genomiky jen malý vliv.

Struktura populace s příměsí

Pro odhalení celkové struktury populace s potenciální příměsí jsme ořezali SNP odstraněním těch s vysokou vazebnou nerovnováhou a potenciálním funkčním vlivem. Analýza vícerozměrného škálování (MDS) založená na ořezané podmnožině reprodukovala podobný výsledek jako celý soubor (obr. 2a a doplňkový obr. 8). Podle očekávání bylo možné dromedáry a divoké velbloudy baktrijské oddělit pomocí první, resp. druhé souřadnice. Po zahrnutí třetí souřadnice do MDS byl IRAN oddělen od všech ostatních domácích velbloudů baktrijských (obr. 2a).

Pro odhad různých podílů předků jsme provedli analýzu struktury populace pomocí Admixture35 za předpokladu K předků (obr. 2b). Postup křížové validace potvrdil, že K = 3 je optimální (doplňkový obr. 9), což ukazuje jasné rozdělení na dromedáry, divoké velbloudy baktrijské a domácí velbloudy baktrijské. Byla pozorována zjevná introgrese domácích velbloudů baktrijských do íránských dromedárů, přinejmenším u jednoho dromedára. V souladu s tím byl dromedárský původ rozšířen v populacích baktrijských velbloudů ve střední Asii, včetně ÍRÁNU, KAZA a RUSKA, s podílem odhadovaným na 1-10 %. Kromě toho jsme u několika volně žijících jedinců pozorovali předky domácích velbloudů baktrijských s podílem 7-15 %. To by mohlo vyplývat z polymorfismu předků, ale mohlo by to být způsobeno také introgresivní hybridizací, která byla pozorována u mtDNA36 a Y chromozomů15 a podle předpokladu ohrožuje genetickou odlišnost divokých druhů. Překvapivě se divocí velbloudi téměř nijak nepodíleli na původu domácích populací, a to ani mongolských populací, které s divokými velbloudy sdílejí blízké prostředí. Ačkoli většina domácích velbloudů baktrijských postrádala diferenciaci při růstu K, IRAN byla první populací, která se oddělila s jedinečným původem (K = 5; obr. 2b).

Jako další metodu zkoumání populační struktury s příměsí jsme odvodili populační strom pro velbloudy pomocí TreeMix37 (obr. 2c). Když nebyla přidána žádná migrační stopa, topologie stromu opět ukázala, že IRAN byl první populací oddělenou mezi všemi domácími velbloudy baktrijskými (doplňkový obr. 10). Zvýšení počtu migračních stop (m = 1-3) mohlo výrazně zlepšit fit modelu (Doplňkový obr. 11), který identifikoval genové toky od dromedárů k domácím velbloudům baktrijským ve střední Asii, včetně KAZA, RUS a IRÁNU s migračními váhami od 4 % do 9 % (Doplňková tabulka 8). Za zmínku stojí, že ačkoli migrační stopa směřovala na konec dromedárské větve do IRÁNU, směřovala doprostřed dromedárské větve do KAZA a RUS (obr. 2c). To by mohlo naznačovat existenci přízračné populace příbuzné íránskému dromedárovi, která přispěla k původu KAZA a RUS. Další migrační stopa (m = 4) by mohla nadále zlepšovat shodu modelu, který naznačoval migraci dromedára do plemene XJ (obr. 2c). Ačkoli TreeMix nezjistil žádný silný signál migrace mezi divokými a domácími velbloudy baktrijskými, zbytky ukázaly mírnou příměs mezi divokými a východoasijskými plemeny (doplňkový obr. 11). K vyhodnocení statistické významnosti těchto příměsí jsme pak použili méně parametrizovaný test tří a čtyř populací (F3/F4)38 . Test F3 opět silně podpořil příměs dromedárů a velbloudů baktrijských v KAZA, RUS a IRÁNU (Doplňková tabulka 9). Citlivější test F4 potvrdil významně vyšší rozsah příměsí dromedárů a velbloudů baktrijských ve střední Asii ve srovnání s velbloudy ve východní Asii (doplňková tabulka 10). Mezi posledně jmenovanými byl zjištěn vyšší rozsah příměsi s dromedáry u XJ než u MG/IMG.

Důkazy pro středoasijský původ odstraněním introgrese

Východní a střední Asie byly na základě archeologických důkazů1,12,17 dvě alternativní oblasti domestikace baktrijských velbloudů, ale ta nejpravděpodobnější zůstala nevyřešena. Přestože jsme pozorovali největší genetickou diferenciaci mezi íránskou populací a všemi ostatními domestikovanými, existence příměsí mezi dromedáry a baktrijskými velbloudy ve Střední Asii by oslabila podporu pro závěr o původu. Abychom tento efekt snížili, pokusili jsme se odstranit introgresní segmenty dromedárů z genomů velbloudů baktrijských pomocí testu „BABA/ABBA „39 . Rozdělili jsme východoasijské a středoasijské populace do skupin a porovnali sdílení alel mezi oběma skupinami s dromedáry (obr. 3a). Vzhledem k tomu, že předci baktrijských velbloudů u jednoho dromedára, stejně jako předci domácích velbloudů u tří divokých (obr. 2b), by byli matoucími faktory, odstranili jsme z analýzy tyto čtyři jedince. K vyhledání introgresovaných segmentů40 jsme použili statistiku fd, robustní verzi Pattersonova D, a pomocí procedury Jackknife jsme použili přísnou hladinu významnosti Z-skóre = 2 (doplňkový obr. 12). V celkem 21 153 nepřekrývajících se 100 kb segmentech napříč genomem bylo podle očekávání mnohem více segmentů s předpokládanými signály introgrese u středoasijských populací (11 711, Z-skóre > 2) než u východoasijských populací (3891, Z-skóre < -2). Poté jsme provedli analýzu příměsí na základě zbývajících segmentů a potvrdili, že introgrese dromedárů byla účinně redukována (doplňkový obr. 13). Přepočet párové Fst po odstranění introgrese stále ukázal, že IRAN je ze všech domácích populací nejvíce diferencovaná (0,04-0,06) (obr. 3b). Abychom získali více informací o fylogenezi populací, rekonstruovali jsme NJ strom na základě párového Fst a provedli bootstrap test (obr. 3b). Ten potvrdil, že mezi všemi domácími populacemi baktrijců se jako první oddělil IRAN, následovaný KAZA a RUS. Bayesovská binární analýza Markov Chain Monte Carlo (MCMC) založená na fylogenezi silně podpořila Střední Asii jako oblast předků domácích velbloudů baktrijských (pravděpodobnost = 99,78 %) a následnou disperzní cestu ze Střední do Východní Asie (doplňkový obr. 14).

Jako nezávislý důkaz jsme také rekonstruovali strom maximální pravděpodobnosti mtDNA plné délky na základě 128 vzorků, které jsme sekvenovali v této studii, a také 39 dalších vzorků dostupných z Genbank (obr. 3c a doplňková tabulka 11). Introgresi mtDNA bylo možné snadno identifikovat a vyloučit ze stromu. Například dvě nově sekvenované mtDNA z KAZA a RUS byly shlukovány s dromedáry. V rámci kladu domácích velbloudů baktrijských, přestože většina velbloudů z různých geografických oblastí byla smíšená, tvořily dvě mtDNA z IRÁNU nejzákladnější větve domácích populací (obr. 3c). Bayesovská binární analýza MCMC opět podpořila středoasijský původ domácích velbloudů baktrijských (pravděpodobnost = 76,43 %).

Demografická historie velbloudů baktrijských

Provedli jsme několik parametrických modelových analýz, abychom odvodili demografickou dynamiku velbloudů v historii. V souladu s předchozí studií10 odhalily dlouhodobé trajektorie velbloudů baktrijských na základě párového sekvenčně markovského koalescenčního modelu (PSMC)41 obrovský pokles efektivní velikosti populace předků velbloudů dříve než před jedním milionem let (doplňkový obr. 15). Přestože dlouhodobé trajektorie divokých a domácích velbloudů baktrijských byly obecně podobné, bylo zřejmé, že se od sebe lišily již před 0,4 milionu let, což vylučuje první z nich jako přímé předchůdce druhých jako předchozí analýzy mtDNA9,14.

Pro podrobnější prozkoumání doby divergence mezi populacemi velbloudů jsme použili zobecněný fylogenetický koalescentní vzorkovač (G-PhoCS)42 . Vzhledem k fylogenezi velbloudích populací dokázal G-PhoCS odhadnout velikost mutačně škálované populace a dobu divergence populace na základě nespojitých neutrálních lokusů v jednotlivých genomech z každé populace. Abychom snížili složitost modelu, zahrnuli jsme pouze dromedáry, divoké velbloudy baktrijské a tři reprezentativní populace (IRAN, KAZA a MG) domácích velbloudů baktrijských (doplňkový obr. 16 a doplňková tabulka 12). Podle obr. 3b byly IRAN a KAZA prvními dvěma středoasijskými populacemi, které se oddělily, a rozdělení MG by mohlo naznačovat rozptyl ze střední do východní Asie. Stáří bylo kalibrováno za předpokladu divergence velblouda baktrijského a dromedára 5,73 milionu let podle databáze Timetree43. Pokud nebylo zahrnuto žádné migrační pásmo, bylo možné snadno dosáhnout konvergence všech odhadů parametrů (doplňkový obr. 17 a doplňková tab. 13). Podobně jako u výsledků PSMC se efektivní velikost populace obecně snižovala od ancestrálních populací k moderním (obr. 4). Doba divergence mezi divokými a domácími velbloudy baktrijskými byla odhadnuta na 0,43 milionu let (95% interval spolehlivosti : 0,13-0,73 Mya; obr. 4), což bylo o něco později než na základě mtDNA (0,714 nebo 1,1 Mya9). Z domácích populací se IRAN oddělil od ostatních asi před 4,45 tisíci lety (95% CI: 0,07-17,6 Kya) a poté se středoasijské a východoasijské populace oddělily asi před 2,40 tisíci lety (95% CI: 0,01-7,84 Kya; obr. 4).

Abychom umožnili tok genů, pokusili jsme se také zavést migrační pásy z dromedárů do populací velbloudů baktrijských (doplňkový obr. 16 a doplňková tabulka 12). Odhady se mohly sblížit pouze tehdy, když byl zaveden migrační pás z íránských dromedárů do IRÁNU a migrační pás z populace ducha do KAZA (doplňkový obr. 18 a doplňková tabulka 13). Ačkoli doba divergence mezi divokými a domácími velbloudy baktrijskými se s migračním modelem nezměnila (0,46 Mya, 95% CI: 0,24-0,71 Mya), první doba divergence domácích populací (0,19 Mya, 95% CI: 0,08-0,31 Mya) se stala nereálnou, protože byla daleko za známou historií domestikace hospodářských zvířat (11,5 Kya44). Kromě toho byla celková míra migrace odhadnuta pouze na 0,27 % a 0,16 % pro migrační pásmo do IRÁNU, resp. do KAZY, což je mnohem méně než míra migrace odhadnutá pomocí příměsi (1-10 %). Možným důvodem tohoto špatného odhadu by mohla být složitější historie příměsí, než jakou předpokládá model kontinuální migrace s konstantními mírami v systému G-PhoCS.