Helgenomsekvensering av 128 kameler i Asien avslöjar ursprung och migration av tama baktriska kameler

Provsamling och helgenomsekvensering

Totalt 105 tama baktriska kameler i Asien, 19 vilda baktriska kameler från Gobi-Altai-regionen i MG samt 4 dromedarer från Iran samlades in för den här studien (kompletterande fig. 1 och kompletterande tabell 1). De tama baktriska kamelerna valdes ut för att täcka så många större geografiska områden som möjligt, inklusive 55 från Inre MG (IMG), Xinjiang (XJ) och Qinghai i Kina, 28 från MG, 6 från KAZA, 10 från Ryssland (RUS) och 6 från IRAN. Eftersom en mängd olika lokala raser bildades på grund av det breda utnyttjandet av tama baktriska kameler i Kina och MG, valdes åtta olika representativa raser ut från regionerna. De övriga tama baktriska kameler från Centralasien levde runt Kaspiska havet.

Efter DNA-extraktion sekvenserades enskilda genomer till en genomsnittlig 13× täckning (kompletterande figur 2 och kompletterande tabell 2). Sekvensavläsningarna anpassades till vår tidigare genomsamling av den baktriska kamelen29 för att kunna identifiera varianter. Efter sträng filtrering (kompletterande figur 3) identifierade vi totalt 13,83 miljoner enskilda nukleotidpolymorfismer (SNP) och 1,41 miljoner små indels. Det är värt att notera att förhållandet mellan övergång och transversion för råa SNP (2,29) var lägre än det som rapporterats för dromedarer (2,31-2,34)31 , men det ökade till 2,44 efter filtreringsprocedurerna, vilket tyder på en kvalitetsförbättring av de identifierade SNP:erna. Den funktionella annoteringen av varianterna visade att cirka 63,10 % av dem var intergeniska, 33,62 % var introniska och 0,94 % var exoniska (kompletterande tabell 3). Det fanns 13,73 miljoner, 6,39 miljoner och 10,55 miljoner identifierade varianter hos tama baktriska kameler, vilda baktriska kameler respektive dromedarer. Även om dromedarer var mer avvikande från båda de baktriska kamelarterna i fylogeni, delade de tama baktriska kamelerna fler varianter med dromedarer (66,73 %) än med vilda baktriska kameler (39,31 %) (kompletterande figur 4) på grund av den enorma minskningen av genetiska varianter som observerats hos de vilda baktriska kamelerna och på grund av genflödet mellan dromedarer och tama baktriska kameler. Bland de tama baktriska kamelerna identifierades 12,68 miljoner och 11,61 miljoner varianter i de östasiatiska respektive centralasiatiska populationerna (kompletterande figur 4). Även om de tama kameler som provtagits från Östasien var fler än de från Centralasien, visade variantantalet privat för varje population ingen signifikant bias mellan de två områdena (P-värde = 0,77, tvåsidigt t-test; kompletterande tabell 4).

Genisk mångfald och differentiering

För en mer detaljerad jämförelse av den genetiska mångfalden mellan olika populationer tog vi först bort 14 individer som uppvisade ett nära genetiskt släktskap med de återstående andra (kompletterande tabell 5). Den parvisa nukleotiddiversiteten π (fig. 1a) hos dromedarer (1,54 × 10-3) var betydligt högre än hos baktriska kameler från alla geografiska områden (0,88 × 10-3-1,11 × 10-3; kompletterande tabell 6), vilket stod i kontrast till tidigare uppskattningar av heterozygositeten baserade på individuella genomer10. En viktig orsak kan vara hybridiseringspraxis mellan dromedarer och baktriska kameler i Centralasien30. Bland de baktriska kamelerna uppvisade den vilda populationen den lägsta π (0,88 × 10-3) jämfört med alla domesticerade populationer (fig. 1a och kompletterande tabell 6). Även om detta resultat strider mot många fall där vilda djur vanligtvis har högre genetisk mångfald än sina tama motsvarigheter, t.ex. hundar24, grisar27 och kaniner32 , skulle det inträffa i fallet med utrotningshotade vilda djur med extremt liten populationsstorlek, t.ex. hästar33. Dessutom uppvisade de husdjurspopulationer som lever i Centralasien i allmänhet en högre mångfald (1,03 × 10-3-1,11 × 10-3) än de som lever i Östasien (0,95 × 10-3-1,02 × 10-3; fig. 1a). Tendensen bekräftades också av Wattersons θ (kompletterande figur 5). Återigen skulle hybridisering med dromedarer i Centralasien kunna förklara den högre mångfalden i regionen.

Vi mätte sedan det parvisa genetiska avståndet mellan kamelpopulationerna genom Weirs Fst (Fig. 1b). Resultatet stämde väl överens med den kända fylogenin, som visade att dromedarerna hade den högsta Fst med baktriska kameler (0,54-0,64) och de vilda baktriska kamelerna hade den näst högsta Fst med de tama kamelerna (0,27-0,31). Differentieringen bland de tama baktriska kamelerna var mycket lägre, vilket är i linje med att de nyligen hade ett enda ursprung. Intressant nog var det bland de tama baktriska kamelerna de från IRAN som uppvisade den största divergensen med alla andra (0,05-0,06). För att bekräfta populationsdifferentieringen konstruerade vi ett NJ-träd (neighbor-joining) för alla individer baserat på deras parvisa IBS-matris (identity-by-state) (kompletterande figur 6). NJ-trädet stödde också en monofyletisk klad av alla tama baktriska kameler, inom vilken IRAN bildade den djupaste grenen.

Ett potentiellt problem med de populationsgenomiska uppskattningarna var referensbias, där användningen av ett enda referensgenom skulle leda till låg effektivitet i variantkallning för individer som skiljer sig mycket från det34. För att undersöka bias jämförde vi det saknade antalet varianter mellan de tre arterna, med sekvenseringsdjupet som en kovariat (kompletterande tabell 7). Variansanalysen (ANOVA) visade att även om de tama baktriska kamelerna inte hade någon signifikant skillnad mot de vilda kamelerna (P-värde = 0,50), hade de faktiskt ett lägre antal saknade varianter än dromedarer (P-värde = 4,38 × 10-3). För att utvärdera biasens inverkan på våra uppskattningar räknade vi om den genetiska mångfalden och Fst med enbart synonyma SNP:er (kompletterande figur 7), eftersom kodande sekvenser med större sannolikhet var invarianta mellan arter. Som ett resultat var skattningarna baserade på de synonyma SNP:erna i god överensstämmelse med hela genomet för alla arter, vilket tyder på att referensbiaset endast hade mindre effekter på våra populationsgenomiska skattningar.

Populationsstruktur med blandning

För att avslöja den övergripande populationsstrukturen med potentiell blandning beskärde vi SNP:erna genom att ta bort de SNP:er som hade hög länkningssjukvärdigdom och potentiella funktionella effekter. MDS-analysen (multidimensional scaling) baserad på den beskurna delmängden gav liknande resultat som den fullständiga mängden (fig. 2a och kompletterande fig. 8). Som förväntat kunde dromedarer och vilda baktriska kameler skiljas åt med hjälp av den första respektive andra koordinaten. När den tredje koordinaten införlivades i MDS separerades IRAN från alla andra tama baktriska kameler (fig. 2a).

För att skatta olika proportioner av förfäder utförde vi en analys av befolkningsstrukturer med Admixture35 genom att anta K förfäderpopulationer (fig. 2b). Korsvalideringsproceduren stödde att K = 3 var optimalt (kompletterande figur 9), vilket visar en tydlig uppdelning mellan dromedarer, vilda baktriska kameler och tama baktriska kameler. Tydlig introgression av tama baktriska kameler till iranska dromedarer observerades, åtminstone hos en dromedar. Följaktligen var dromedarfamiljen vanlig i de centralasiatiska populationerna av baktriska kameler, inklusive IRAN, KAZA och RUS, med en andel som uppskattas till 1-10 %. Dessutom observerade vi anor från tama baktriska kameler hos flera vilda individer med en andel på 7-15 %. Detta kan bero på polymorfismer hos förfäderna, men det kan också orsakas av introgressiv hybridisering, som observerades med mtDNA36 och Y-kromosomer15 , och som föreslogs hota de vilda arternas genetiska särprägel. Överraskande nog bidrog de vilda kamelerna nästan inte alls till de inhemska populationernas anor, inte ens till de mongoliska populationerna, som delar nära livsmiljöer med de vilda kamelerna. Även om de flesta tama baktriska kameler saknade differentiering när K växte, var IRAN den första populationen som separerade med en unik härstamning (K = 5; fig. 2b).

Som en annan metod för att undersöka populationsstrukturen med inblandning, härledde vi populationsträdet för kamelerna med hjälp av TreeMix37 (fig. 2c). När inget migrationsspår lades till visade trädets topologi återigen att IRAN var den första populationen som separerades bland alla inhemska baktriska kameler (kompletterande fig. 10). Om man ökade antalet migrationsspår (m = 1-3) kunde man avsevärt förbättra modellens anpassning (kompletterande figur 11), vilket identifierade genflöden från dromedarer till tama baktriska kameler i Centralasien, inklusive KAZA, RUS och IRAN med migrationsvikter som varierade från 4 % till 9 % (kompletterande tabell 8). Det var värt att nämna att även om migrationsspåret pekade i slutet av dromedargrenen till IRAN, pekade det i mitten av dromedargrenen till KAZA och RUS (fig. 2c). Detta skulle kunna innebära att det finns en spökpopulation som är relaterad till den iranska dromedären och som bidrog till KAZA:s och RUS:s härstamning. Ytterligare migrationsspår (m = 4) kunde fortsätta att förbättra modellens passform, vilket tydde på att dromedaren migrerade till en XJ-ras (fig. 2c). Även om TreeMix inte upptäckte någon stark signal om migration mellan vilda och tama baktriska kameler, visade resterna en måttlig blandning mellan vilda och östasiatiska raser (kompletterande figur 11). Vi använde sedan det mindre parametriserade tre- och fyrapopulationstestet (F3/F4)38 för att utvärdera den statistiska betydelsen av dessa blandningshändelser. Återigen gav F3-testet starkt stöd för blandningen av dromedarer och baktriska kameler i KAZA, RUS och IRAN (kompletterande tabell 9). Det känsligare F4-testet bekräftade en betydligt högre grad av blandning mellan dromedarer och baktriska kameler i Centralasien jämfört med dem i Östasien (kompletterande tabell 10). Bland de sistnämnda upptäcktes en högre grad av blandning med dromedarer i XJ än i MG/IMG.

Bevis för centralasiatiskt ursprung genom att ta bort introgression

Öst- och Centralasien var de två alternativa regionerna för domesticering av baktriska kameler baserat på arkeologiska bevis1,12,17, men den mest sannolika förblev olöst. Även om vi observerade den största genetiska differentieringen mellan den iranska populationen och alla andra domesticerade populationer, skulle förekomsten av blandning mellan dromedarer och baktriska kameler i Centralasien försvaga stödet för ursprungsslutsatsen. För att minska denna effekt försökte vi avlägsna de introgresserade segmenten från dromedarer från de baktriska kamelernas genomer med hjälp av ”BABA/ABBA”-testet39. Vi grupperade de öst- och centralasiatiska populationerna och jämförde alleldelning mellan de två grupperna och dromedarer (fig. 3a). Eftersom anor av baktriska kameler i en dromedar och anor av tama kameler i tre vilda kameler (fig. 2b) skulle vara förväxlingsfaktorer, tog vi bort de fyra individerna i analysen. Vi använde statistiken fd, en robust version av Pattersons D, för att lokalisera introgresserade segment40 och tillämpade en strikt signifikansnivå på Z-score = 2 med hjälp av Jackknife-proceduren (kompletterande figur 12). I totalt 21 153 icke-överlappande 100 kb-segment över hela genomet fanns det som förväntat betydligt fler segment med förmodade signaler på introgression i de centralasiatiska populationerna (11 711, Z-score > 2) än i de östasiatiska populationerna (3891, Z-score < -2). Vi utförde sedan Admixture-analysen baserad på de återstående segmenten och bekräftade att introgressionen av dromedarer var effektivt reducerad (kompletterande figur 13). En ny beräkning av de parvisa Fst efter att introgressionen tagits bort visade fortfarande att IRAN var den mest differentierade (0,04-0,06) bland alla inhemska populationer (fig. 3b). För att få mer insikt i populationsfylogenin rekonstruerade vi NJ-trädet baserat på det parvisa Fst och utförde bootstrap-testet (fig. 3b). Det bekräftade att IRAN var den första som separerades bland alla inhemska baktriska populationer, följt av KAZA och RUS. Den Bayesianska binära Markov Chain Monte Carlo-analysen (MCMC) baserad på fylogenin stödde starkt Centralasien som det ursprungliga området för tama baktriska kameler (sannolikhet = 99,78 %) och en efterföljande spridningsväg från Central- till Östasien (kompletterande fig. 14).

Som ett oberoende bevis rekonstruerade vi också maximum likelihood-trädet för mtDNA i full längd baserat på de 128 prover som vi sekvenserade i den här studien, samt 39 ytterligare prover som finns tillgängliga från Genbank (fig. 3c och kompletterande tabell 11). Introgression av mtDNAs kunde lätt identifieras och uteslutas från trädet. Till exempel klustrades två nyligen sekvenserade mtDNA från KAZA och RUS med dromedarer. Inom kladen av tama baktriska kameler, även om de flesta kameler från olika geografiska områden var blandade, bildade två mtDNA från IRAN de mest basala grenarna av de inhemska populationerna (fig. 3c). Den Bayesianska binära MCMC-analysen stödde återigen det centralasiatiska ursprunget för tama baktriska kameler (sannolikhet = 76,43 %).

Baktriska kamelers demografiska historia

Vi utförde flera parametriska modelleringsanalyser för att härleda kamelernas demografiska dynamik i historien. I överensstämmelse med tidigare studier10 visade de långsiktiga banorna för baktriska kameler baserade på PSMC-modellen (pairwise sequentially Markovian coalescent)41 en enorm minskning av den effektiva populationsstorleken hos de ursprungliga kamelerna tidigare än för en miljon år sedan (kompletterande fig. 15). Även om de långsiktiga banorna för de vilda och tama baktriska kamelerna i allmänhet var likartade, var det uppenbart att de divergerade från varandra så tidigt som för 0,4 miljoner år sedan, vilket utesluter de förstnämnda som direkta stamfäder till de sistnämnda som tidigare mtDNA-analyser9,14.

För att utforska divergenstiden mellan kamelpopulationerna mer i detalj använde vi den generaliserade fylogenetiska koalescenssamlaren (G-PhoCS)42. Med tanke på kamelpopulationernas fylogeni kunde G-PhoCS uppskatta den mutationsskalerade populationsstorleken och populationsdivergenstiden baserat på obundna neutrala loci i enskilda genomer från varje population. För att minska modellens komplexitet inkluderade vi endast dromedarer, vilda baktriska kameler och tre representativa populationer (IRAN, KAZA och MG) av tama baktriska kameler (kompletterande figur 16 och kompletterande tabell 12). Enligt fig. 3b var IRAN och KAZA de två första centralasiatiska populationerna som skiljdes åt och MG:s uppdelning kan tyda på spridning från Central- till Östasien. Åldern kalibrerades genom att anta en divergens mellan den baktriska kamelen och dromedarerna på 5,73 miljoner år enligt Timetree-databasen43. När inget migrationsband införlivades kunde konvergens för alla parameteruppskattningar lätt uppnås (kompletterande figur 17 och kompletterande tabell 13). I likhet med PSMC-resultaten minskade den effektiva populationsstorleken i allmänhet från förfäder till moderna populationer (fig. 4). Divergenstiden mellan vilda och tama baktriska kameler uppskattades till 0,43 miljoner år sedan (95 % konfidensintervall: 0,13-0,73 Mya; fig. 4), vilket var något senare än den som baseras på mtDNA (0,714 eller 1,1 Mya9). Bland de inhemska populationerna separerades IRAN från andra för cirka 4,45 tusen år sedan (95 % CI: 0,07-17,6 Kya) och därefter separerades de central- och östasiatiska populationerna för cirka 2,40 tusen år sedan (95 % CI: 0,01-7,84 Kya; fig. 4).

För att möjliggöra genflöde försökte vi också införa migrationsband från dromedarer till baktriska kamelpopulationer (kompletterande figur 16 och kompletterande tabell 12). Uppskattningarna kunde endast konvergera när ett migrationsband från iranska dromedarer till IRAN och ett migrationsband från en spökpopulation till KAZA infördes (kompletterande figur 18 och kompletterande tabell 13). Även om divergenstiden mellan vilda och tama baktriska kameler inte ändrades med migrationsmodellen (0,46 Mya, 95 % CI: 0,24-0,71 Mya), blev den första divergenstiden för de tama populationerna (0,19 Mya, 95 % CI: 0,08-0,31 Mya) orealistisk, eftersom den låg långt bortom den kända historien om domesticering av boskap (11,5 Kya44). Dessutom uppskattades den totala migrationstakten endast till 0,27 % och 0,16 % för migrationsbandet till IRAN respektive KAZA, vilket är mycket lägre än vad som uppskattades med hjälp av blandning (1-10 %). En möjlig orsak till den dåliga uppskattningen skulle kunna vara en mer komplex blandningshistoria än den kontinuerliga migrationsmodell med konstanta hastigheter som antas av G-PhoCS.