Il sequenziamento dell’intero genoma di 128 cammelli in tutta l’Asia rivela l’origine e la migrazione dei cammelli bactriani domestici

Raccolta di campioni e sequenziamento dell’intero genoma

Un totale di 105 cammelli bactriani domestici in Asia, 19 cammelli bactriani selvatici della regione di Gobi-Altai in MG e 4 dromedari dell’IRAN sono stati raccolti per questo studio (Fig. 1 e Tabella 1 supplementare). I cammelli bactriani domestici sono stati scelti per coprire il maggior numero possibile di regioni geografiche, tra cui 55 da Inner MG (IMG), Xinjiang (XJ), e Qinghai della Cina, 28 da MG, 6 da KAZA, 10 dalla Russia (RUS) e 6 dall’IRAN. Poiché si è formata una varietà di razze locali a causa dell’ampio utilizzo di cammelli bactriani domestici in Cina e MG, sono state scelte otto diverse razze rappresentative dalle regioni. Gli altri cammelli bactriani domestici dell’Asia centrale vivevano intorno al Mar Caspio.

Dopo l’estrazione del DNA, i singoli genomi sono stati sequenziati con una copertura media di 13× (Fig. 2 supplementare e Tabella 2 supplementare). Le letture di sequenza sono state allineate al nostro precedente assemblaggio del genoma del cammello bactriano29 per la chiamata delle varianti. Dopo un filtraggio rigoroso (Fig. 3 supplementare), abbiamo identificato totalmente 13,83 milioni di polimorfismi a singolo nucleotide (SNPs) e 1,41 milioni di piccoli indel. In particolare, il rapporto transizione-trasversione degli SNP grezzi (2.29) era inferiore a quello riportato nei dromedari (2.31-2.34)31, ma è stato aumentato a 2.44 dopo le procedure di filtraggio, suggerendo un miglioramento della qualità degli SNP identificati. L’annotazione funzionale delle varianti ha indicato che circa il 63,10% di esse erano intergeniche, il 33,62% erano introniche e lo 0,94% erano esoniche (tabella 3 supplementare). C’erano 13,73 milioni, 6,39 milioni, e 10,55 milioni di varianti identificate nei cammelli bactriani domestici, nei cammelli bactriani selvatici e nei dromedari, rispettivamente. Sebbene i dromedari fossero più divergenti da entrambe le specie di cammello bactriano nella filogenesi, i cammelli bactriani domestici hanno condiviso più varianti con i dromedari (66,73%) che con il cammello bactriano selvatico (39,31%) (Fig. 4 supplementare) a causa dell’enorme riduzione delle varianti genetiche osservate nel cammello bactriano selvatico esistente e al flusso genico tra i dromedari e i cammelli bactriani domestici. Tra i cammelli bactriani domestici, c’erano 12,68 milioni e 11,61 milioni di varianti identificate nelle popolazioni dell’Asia orientale e dell’Asia centrale, rispettivamente (Fig. 4 supplementare). Anche se i cammelli domestici campionati dall’Asia orientale erano più di quelli dell’Asia centrale, il conteggio delle varianti private di ciascuna popolazione non ha mostrato alcuna distorsione significativa tra le due aree (P-value = 0,77, t-test a due code; Tabella supplementare 4).

Diversità genetica e differenziazione

Per un confronto più dettagliato della diversità genetica tra le diverse popolazioni, abbiamo prima rimosso 14 individui che mostravano una stretta relazione genetica con gli altri rimasti (Tabella supplementare 5). La diversità nucleotidica a coppie π (Fig. 1a) dei dromedari (1,54 × 10-3) era significativamente più alta di quella dei cammelli bactriani di tutte le regioni geografiche (0,88 × 10-3-1,11 × 10-3; Tabella supplementare 6), che era in contrasto con le precedenti stime di eterozigosi basate sui genomi individuali10. Una ragione importante potrebbe essere la pratica di ibridazione tra dromedari e cammelli bactriani in Asia centrale30. Tra i cammelli bactriani, la popolazione selvatica ha mostrato il più basso π (0,88 × 10-3) rispetto a tutte le popolazioni domestiche (Fig. 1a e Tabella 6 supplementare). Anche se questo risultato ha violato molti casi che gli animali selvatici di solito hanno una maggiore diversità genetica rispetto alle loro controparti domestiche come cani24, maiali27, e conigli32, sarebbe successo nel caso di animali selvatici in pericolo con una dimensione della popolazione estremamente piccola come i cavalli33. Inoltre, le popolazioni domestiche che vivono in Asia centrale hanno generalmente mostrato una maggiore diversità (1,03 × 10-3-1,11 × 10-3) rispetto a quelle che vivono in Asia orientale (0,95 × 10-3-1,02 × 10-3; Fig. 1a). La tendenza è stata confermata anche dal θ di Watterson (Fig. 5 supplementare). Ancora una volta, l’ibridazione con i dromedari in Asia centrale potrebbe spiegare la maggiore diversità nella regione.

Abbiamo poi misurato la distanza genetica a coppie tra le popolazioni di cammelli mediante Fst di Weir (Fig. 1b). Il risultato era ben in accordo con la filogenesi conosciuta, che indicava che i dromedari avevano il più alto Fst con i cammelli bactriani (0,54-0,64) e i cammelli bactriani selvatici avevano il secondo più alto Fst con quelli domestici (0,27-0,31). La differenziazione tra i cammelli bactriani domestici era molto più bassa, in linea con la loro recente origine unica. È interessante notare che tra i cammelli bactriani domestici, quelli provenienti dall’IRAN hanno mostrato la maggiore divergenza con tutti gli altri (0,05-0,06). Per convalidare la differenziazione della popolazione, abbiamo costruito un albero neighbor-joining (NJ) per tutti gli individui basato sulla loro matrice di identità a coppie (IBS) (Fig. 6 supplementare). L’albero NJ ha anche supportato un clade monofiletico di tutti i cammelli bactriani domestici, all’interno del quale IRAN formava il ramo più profondo.

Un potenziale problema con le stime genomiche della popolazione era il bias di riferimento, dove l’utilizzo di un singolo genoma di riferimento avrebbe portato ad una bassa efficienza nella chiamata delle varianti per gli individui che differiscono molto da esso34. Per indagare il bias, abbiamo confrontato il conteggio mancante di varianti tra le tre specie, prendendo la profondità di sequenziamento come covariata (Tabella supplementare 7). L’analisi della varianza (ANOVA) ha mostrato che anche se i cammelli bactriani domestici non avevano alcuna differenza significativa con quelli selvatici (P-value = 0.50), hanno effettivamente avuto un numero inferiore di mancanti rispetto ai dromedari (P-value = 4.38 × 10-3). Per valutare l’impatto del bias sulle nostre stime, abbiamo ricalcolato la diversità genetica e Fst con solo gli SNPs sinonimi (Fig. 7 supplementare), poiché le sequenze codificanti hanno maggiori probabilità di essere invarianti tra le specie. Come risultato, le stime basate sugli SNP sinonimi erano in buon accordo con l’intero genoma per tutte le specie, suggerendo che il bias di riferimento ha avuto solo effetti minori sulle nostre stime genomiche della popolazione.

Struttura della popolazione con commistione

Per rivelare la struttura complessiva della popolazione con potenziale commistione, abbiamo sfrondato gli SNP rimuovendo quelli con alto linkage disequilibrium e potenziali effetti funzionali. L’analisi di scala multidimensionale (MDS) basata sul sottoinsieme sfrondato ha riprodotto un risultato simile al set completo (Fig. 2a e Fig. 8 supplementare). Come previsto, i dromedari e i cammelli bactriani selvatici potevano essere separati dalla prima e dalla seconda coordinata, rispettivamente. Quando la terza coordinata è stata incorporata nel MDS, IRAN è stato separato da tutti gli altri cammelli bactriani domestici (Fig. 2a).

Per stimare diverse proporzioni ancestrali, abbiamo eseguito l’analisi della struttura della popolazione con Admixture35 assumendo K popolazioni ancestrali (Fig. 2b). La procedura di convalida incrociata ha sostenuto che K = 3 era ottimale (Fig. 9 supplementare), mostrando una chiara divisione tra i dromedari, i cammelli bactriani selvatici e i cammelli bactriani domestici. È stata osservata un’evidente introgressione dei cammelli bactriani domestici nei dromedari iraniani, almeno in un dromedario. Di conseguenza, l’ascendenza dromedaria era prevalente nelle popolazioni di cammelli bactriani dell’Asia centrale tra cui IRAN, KAZA e RUS, con una proporzione stimata come 1-10%. Inoltre, abbiamo osservato l’ascendenza dei cammelli bactriani domestici in diversi individui selvatici con una percentuale del 7-15%. Questo potrebbe derivare da polimorfismi ancestrali, ma potrebbe anche essere causato dall’ibridazione introgressiva, che è stata osservata con mtDNA36 e cromosomi Y15, ed è stata proposta per minacciare la distintività genetica delle specie selvatiche. Sorprendentemente, i cammelli selvatici non hanno contribuito quasi per niente all’ascendenza delle popolazioni domestiche, anche alle popolazioni mongole, che condividono habitat vicini a quelli dei cammelli selvatici. Anche se la maggior parte dei cammelli bactriani domestici mancava di differenziazione quando K cresceva, IRAN è stata la prima popolazione a separarsi con un antenato unico (K = 5; Fig. 2b).

Come un altro metodo per esaminare la struttura della popolazione con commistione, abbiamo dedotto l’albero della popolazione per i cammelli usando TreeMix37 (Fig. 2c). Quando nessuna traccia di migrazione è stata aggiunta, la topologia dell’albero ha indicato ancora una volta che IRAN era la prima popolazione separata tra tutti i cammelli bactriani domestici (Fig. 10 supplementare). L’aumento delle tracce di migrazione (m = 1-3) potrebbe migliorare notevolmente l’adattamento del modello (Fig. 11 supplementare), che ha identificato i flussi genici dai dromedari ai cammelli bactriani domestici in Asia centrale, tra cui KAZA, RUS, e IRAN con pesi di migrazione che vanno dal 4% al 9% (Tabella 8 supplementare). Vale la pena menzionare che anche se la traccia di migrazione puntava alla fine del ramo dei dromedari verso l’IRAN, essa puntava al centro del ramo dei dromedari verso KAZA e RUS (Fig. 2c). Questo potrebbe implicare una popolazione fantasma legata al dromedario iraniano che ha contribuito all’ascendenza di KAZA e RUS. Un’ulteriore traccia di migrazione (m = 4) potrebbe continuare a migliorare l’adattamento del modello, che indicava la migrazione del dromedario ad una razza XJ (Fig. 2c). Anche se TreeMix non ha rilevato alcun segnale forte di migrazione tra i cammelli bactriani selvatici e domestici, i residui hanno mostrato una moderata commistione tra le razze selvatiche e dell’Asia orientale (Fig. 11 supplementare). Abbiamo quindi utilizzato il test a tre e quattro popolazioni (F3/F4) meno parametrizzato38 per valutare la significatività statistica di questi eventi di commistione. Ancora una volta, il test F3 ha supportato fortemente la commistione tra dromedari e cammelli bactriani in KAZA, RUS e IRAN (tabella 9 supplementare). Il test F4, più sensibile, ha confermato un grado di commistione significativamente più elevato tra dromedari e cammelli bactriani in Asia centrale rispetto a quelli dell’Asia orientale (Tabella supplementare 10). Tra questi ultimi, una maggiore estensione di commistione con i dromedari è stata rilevata in XJ che in MG/IMG.

Prova dell’origine dell’Asia centrale rimuovendo l’introgressione

L’Asia orientale e centrale erano le due regioni alternative di addomesticamento dei cammelli bactriani basate su prove archeologiche1,12,17, ma la più probabile è rimasta irrisolta. Anche se abbiamo osservato la più grande differenziazione genetica tra la popolazione iraniana e tutte le altre domestiche, l’esistenza di commistione tra dromedari e cammelli bactriani in Asia centrale indebolirebbe il supporto per l’inferenza di origine. Per ridurre questo effetto, abbiamo cercato di rimuovere i segmenti introgressivi dei dromedari dai genomi dei cammelli bactriani utilizzando il test “BABA/ABBA “39. Abbiamo raggruppato le popolazioni dell’Asia orientale e centrale, e abbiamo confrontato la condivisione degli alleli tra i due gruppi con i dromedari (Fig. 3a). Poiché l’ascendenza dei cammelli bactriani in un dromedario, così come l’ascendenza dei cammelli domestici in tre selvatici (Fig. 2b), sarebbero fattori confondenti, abbiamo rimosso i quattro individui nell’analisi. Abbiamo usato la statistica fd, una versione robusta della D di Patterson, per individuare i segmenti introgressi40 e abbiamo applicato un livello di significatività rigoroso di Z-score = 2 utilizzando la procedura Jackknife (Fig. 12 supplementare). In un totale di 21.153 segmenti non sovrapposti di 100 kb attraverso il genoma, c’erano molti più segmenti con segnali putativi di introgressione nelle popolazioni dell’Asia centrale (11.711, Z-score > 2) che nelle popolazioni dell’Asia orientale (3891, Z-score < -2) come previsto. Abbiamo quindi eseguito l’analisi di commistione basata sui segmenti rimanenti e abbiamo confermato che l’introgressione dei dromedari era effettivamente ridotta (Fig. 13 supplementare). Il ricalcolo delle Fst a coppie dopo aver rimosso l’introgressione indicava ancora che IRAN era la più differenziata (0,04-0,06) tra tutte le popolazioni domestiche (Fig. 3b). Per ottenere maggiori informazioni sulla filogenesi della popolazione, abbiamo ricostruito l’albero NJ basato sulle Fst a coppie ed eseguito il test bootstrap (Fig. 3b). Ha confermato che IRAN è stata la prima a separarsi tra tutte le popolazioni domestiche Bactrian, seguita da KAZA e RUS. L’analisi binaria bayesiana Markov Chain Monte Carlo (MCMC) basata sulla filogenesi ha fortemente sostenuto l’Asia centrale come area ancestrale dei cammelli bactriani domestici (probabilità = 99,78%) e un successivo percorso di dispersione dall’Asia centrale a quella orientale (Fig. 14 supplementare).

Come prova indipendente, abbiamo anche ricostruito l’albero di massima verosimiglianza dei mtDNA full-length basato sui 128 campioni che abbiamo sequenziato in questo studio, così come 39 campioni aggiuntivi disponibili da Genbank (Fig. 3c e Tabella 11 supplementare). L’introgressione di mtDNA potrebbe essere facilmente identificata ed esclusa dall’albero. Per esempio, due mtDNA recentemente sequenziati da KAZA e RUS sono stati raggruppati con i dromedari. All’interno della clade dei cammelli bactriani domestici, anche se la maggior parte dei cammelli provenienti da diverse regioni geografiche erano misti, due mtDNA dall’IRAN formavano i rami più basali delle popolazioni domestiche (Fig. 3c). L’analisi binaria MCMC bayesiana ha nuovamente supportato l’origine centroasiatica dei cammelli bactriani domestici (probabilità = 76,43%).

Storia demografica dei cammelli bactriani

Abbiamo eseguito diverse analisi di modellazione parametrica per dedurre la dinamica demografica dei cammelli nella storia. Coerentemente con lo studio precedente10, le traiettorie a lungo termine dei cammelli bactriani basate sul modello pairwise sequentially Markovian coalescent (PSMC)41 hanno rivelato una tremenda diminuzione della dimensione effettiva della popolazione dei cammelli ancestrali prima di un milione di anni fa (Fig. 15 supplementare). Anche se le traiettorie a lungo termine dei cammelli selvatici e domestici Bactrian erano generalmente simili, erano evidenti a divergere l’uno dall’altro già da 0,4 milioni di anni fa, escludendo i primi come progenitori diretti dei secondi come precedenti analisi mtDNA9,14.

Per esplorare il tempo di divergenza tra le popolazioni di cammelli in modo più dettagliato, abbiamo usato il campionatore generalizzato coalescente filogenetico (G-PhoCS)42. Data la filogenesi delle popolazioni di cammello, G-PhoCS potrebbe stimare la dimensione della popolazione a scala di mutazione e il tempo di divergenza della popolazione basato su loci neutri non collegati nei singoli genomi di ogni popolazione. Per ridurre la complessità del modello, abbiamo incluso solo i dromedari, i cammelli bactriani selvatici e tre popolazioni rappresentative (IRAN, KAZA e MG) di cammelli bactriani domestici (Fig. 16 supplementare e Tabella 12 supplementare). Secondo la Fig. 3b, IRAN e KAZA sono state le prime due popolazioni dell’Asia centrale a separarsi e la scissione di MG potrebbe indicare la dispersione dall’Asia centrale a quella orientale. L’età è stata calibrata assumendo la divergenza del cammello bactriano e del dromedario di 5,73 milioni di anni secondo il database Timetree43. Quando non è stata incorporata alcuna banda di migrazione, la convergenza di tutte le stime dei parametri è stata facilmente raggiunta (Fig. 17 supplementare e Tabella 13 supplementare). Simile ai risultati PSMC, la dimensione effettiva della popolazione era generalmente diminuita dalle popolazioni ancestrali a quelle moderne (Fig. 4). Il tempo di divergenza tra i cammelli bactriani selvatici e domestici è stato stimato come 0,43 milioni di anni fa (intervallo di confidenza al 95%: 0,13-0,73 Mya; Fig. 4), che era leggermente più tardi di quello basato su mtDNA (0,714 o 1,1 Mya9). Tra le popolazioni domestiche, l’IRAN è stato separato dagli altri circa 4,45 mila anni fa (95% CI: 0,07-17,6 Kya) e poi le popolazioni dell’Asia centrale e orientale sono state separate circa 2,40 mila anni fa (95% CI: 0,01-7,84 Kya; Fig. 4).