Recogida de muestras y secuenciación del genoma completo
Un total de 105 camellos bactrianos domésticos de Asia, 19 camellos bactrianos salvajes de la región de Gobi-Altai en MG, así como 4 dromedarios de IRÁN se reunieron para este estudio (Fig. Suplementaria 1 y Tabla suplementaria 1). Los camellos bactrianos domésticos se eligieron para cubrir el mayor número posible de regiones geográficas importantes, incluyendo 55 de MG interior (IMG), Xinjiang (XJ) y Qinghai de China, 28 de MG, 6 de KAZA, 10 de Rusia (RUS) y 6 de IRÁN. Debido a la amplia utilización de los camellos bactrianos domésticos en China y MG se formó una variedad de razas locales, por lo que se eligieron ocho razas representativas diferentes de las regiones. Los otros camellos bactrianos domésticos de Asia Central vivían alrededor del Mar Caspio.
Después de la extracción de ADN, los genomas individuales fueron secuenciados con un promedio de cobertura de 13× (Fig. suplementaria 2 y Tabla suplementaria 2). Las lecturas de la secuencia se alinearon con nuestro anterior ensamblaje del genoma del camello bactriano29 para la llamada de variantes. Tras un riguroso filtrado (Fig. 3), identificamos en total 13,83 millones de polimorfismos de un solo nucleótido (SNP) y 1,41 millones de pequeños indels. En particular, la relación transición-transversión de los SNP sin procesar (2,29) era inferior a la registrada en los dromedarios (2,31-2,34)31 , pero aumentó a 2,44 tras los procedimientos de filtrado, lo que sugiere una mejora de la calidad de los SNP identificados. La anotación funcional de las variantes indicó que aproximadamente el 63,10% de ellas eran intergénicas, el 33,62% intrónicas y el 0,94% exónicas (Tabla suplementaria 3). Se identificaron 13,73 millones, 6,39 millones y 10,55 millones de variantes en los camellos bactrianos domésticos, los camellos bactrianos salvajes y los dromedarios, respectivamente. Aunque los dromedarios eran más divergentes de las dos especies de camellos bactrianos en la filogenia, los camellos bactrianos domésticos compartían más variantes con los dromedarios (66,73%) que con los camellos bactrianos salvajes (39,31%) (Fig. Suplementaria 4) debido a la tremenda reducción de variantes genéticas observada en el camello bactriano salvaje existente y al flujo de genes entre los dromedarios y los camellos bactrianos domésticos. Entre los camellos bactrianos domésticos, se identificaron 12,68 millones y 11,61 millones de variantes en las poblaciones de Asia oriental y central, respectivamente (Fig. 4 suplementaria). Aunque los camellos domésticos muestreados de Asia Oriental eran más que los de Asia Central, el recuento de variantes privado de cada población no mostró un sesgo significativo entre las dos áreas (valor P = 0,77, prueba t de dos colas; Tabla Suplementaria 4).
Diversidad genética y diferenciación
Para una comparación más detallada de la diversidad genética entre las diferentes poblaciones, primero eliminamos 14 individuos que mostraban una estrecha relación genética con los restantes (Tabla Suplementaria 5). La diversidad de nucleótidos por pares π (Fig. 1a) de los dromedarios (1,54 × 10-3) fue significativamente mayor que la de los camellos bactrianos de todas las regiones geográficas (0,88 × 10-3-1,11 × 10-3; Tabla Suplementaria 6), lo que contrasta con las estimaciones previas de heterocigosidad basadas en genomas individuales10. Una razón importante podría ser la práctica de hibridación entre dromedarios y camellos bactrianos en Asia Central30. Entre los camellos bactrianos, la población salvaje mostró la π más baja (0,88 × 10-3) en comparación con todas las poblaciones domésticas (Fig. 1a y Tabla Suplementaria 6). Aunque este resultado infringe muchos casos en los que los animales salvajes suelen tener una mayor diversidad genética que sus homólogos domésticos, como los perros24, los cerdos27 y los conejos32 , esto ocurriría en el caso de animales salvajes en peligro de extinción con un tamaño de población extremadamente pequeño, como los caballos33. Además, las poblaciones domésticas que viven en Asia Central mostraron en general una mayor diversidad (1,03 × 10-3-1,11 × 10-3) que las que viven en Asia Oriental (0,95 × 10-3-1,02 × 10-3; Fig. 1a). Esta tendencia también fue confirmada por la θ de Watterson (Fig. 5 suplementaria). Una vez más, la hibridación con dromedarios en Asia Central podría explicar la mayor diversidad en la región.
a Diversidad de nucleótidos π. El boxplot muestra π para 2,0 × 105 ventanas de 10 kb de deslizamiento a través del genoma. El origen geográfico y el tamaño de la muestra de cada población se muestran a la izquierda y la media de π se muestra a la derecha. Se tomaron muestras de múltiples razas locales para MG, IMG y XJ. Se eliminaron los individuos con relaciones genéticas estrechas. Los elementos del boxplot se definen de la siguiente manera: línea central, mediana; límites de la caja, el tercer y el primer cuartil; bigotes, rango intercuartílico de 1,5. b Diferenciación poblacional por pares Fst. El mapa de calor representa el promedio de Fst para ventanas deslizantes de 2,0 × 105 10 kb. El dendrograma representa la agrupación jerárquica de las poblaciones basada en Fst.
A continuación, medimos la distancia genética por pares entre las poblaciones de camellos mediante Fst de Weir (Fig. 1b). El resultado coincidía con la filogenia conocida, que indicaba que los dromedarios tenían la mayor Fst con los camellos bactrianos (0,54-0,64) y los camellos bactrianos salvajes tenían la segunda mayor Fst con los domésticos (0,27-0,31). La diferenciación entre los camellos bactrianos domésticos fue mucho menor, en consonancia con su reciente origen único. Curiosamente, entre los camellos bactrianos domésticos, los de IRÁN mostraron la mayor divergencia con todos los demás (0,05-0,06). Para validar la diferenciación de la población, construimos un árbol de unión de vecinos (NJ) para todos los individuos basado en su matriz de identidad por estado (IBS) (Fig. 6 suplementaria). El árbol NJ también apoyó un clado monofilético de todos los camellos bactrianos domésticos, dentro del cual IRAN formó la rama más profunda.
Un problema potencial con las estimaciones genómicas de la población fue el sesgo de referencia, donde el uso de un único genoma de referencia conduciría a una baja eficiencia en la llamada de variantes para los individuos que difieren mucho de él34. Para investigar el sesgo, comparamos el recuento de variantes faltantes entre las tres especies, tomando la profundidad de secuenciación como covariable (Tabla Suplementaria 7). El análisis de la varianza (ANOVA) mostró que, aunque los camellos bactrianos domésticos no tenían diferencias significativas con los salvajes (valor P = 0,50), sí tenían un menor recuento de faltas que los dromedarios (valor P = 4,38 × 10-3). Para evaluar el impacto del sesgo en nuestras estimaciones, volvimos a calcular la diversidad genética y la Fst sólo con los SNPs sinónimos (Fig. Suplementaria 7), ya que es más probable que las secuencias codificantes sean invariables entre especies. Como resultado, las estimaciones basadas en los SNPs sinónimos estaban en buen acuerdo con todo el genoma para todas las especies, lo que sugiere que el sesgo de referencia sólo tenía efectos menores en nuestras estimaciones genómicas de la población.
Estructura de la población con mezcla
Para revelar la estructura general de la población con mezcla potencial, podamos los SNPs mediante la eliminación de aquellos con alto desequilibrio de enlace y potenciales efectos funcionales. El análisis de escala multidimensional (MDS) basado en el subconjunto podado reprodujo un resultado similar al del conjunto completo (Fig. 2a y Fig. 8 suplementaria). Como se esperaba, los dromedarios y los camellos bactrianos salvajes podían separarse por la primera y la segunda coordenadas, respectivamente. Cuando se incorporó la tercera coordenada en el MDS, IRÁN se separó de todos los demás camellos bactrianos domésticos (Fig. 2a).
a Gráfico de escalamiento multidimensional (MDS) con las coordenadas 1-4 (C1-C4). b Análisis de admixtura asumiendo diferente número de ancestros K. La proporción del genoma de un individuo asignada a cada ascendencia está representada por diferentes colores. c Análisis TreeMix con diferentes supuestos de eventos de migración m. El peso de la migración es la proporción de ascendencia recibida de la población donante.
Para estimar diferentes proporciones ancestrales, realizamos un análisis de estructura poblacional con Admixture35 asumiendo K poblaciones ancestrales (Fig. 2b). El procedimiento de validación cruzada apoyó que K = 3 era el óptimo (Supplementary Fig. 9), mostrando una clara división entre los dromedarios, los camellos bactrianos salvajes y los camellos bactrianos domésticos. Se observó una introgresión evidente de los camellos bactrianos domésticos en los dromedarios iraníes, al menos en un dromedario. Por consiguiente, la ascendencia de los dromedarios era prevalente en las poblaciones de camellos bactrianos de Asia Central, incluidas las de IRÁN, KAZA y RUS, con una proporción estimada del 1 al 10%. Además, observamos ancestros de camellos bactrianos domésticos en varios individuos salvajes, con una proporción del 7-15%. Esto podría deberse a polimorfismos ancestrales, pero también a la hibridación introgresiva, que se ha observado con el ADNmt36 y los cromosomas Y15, y que se ha propuesto como una amenaza para la distinción genética de las especies salvajes. Sorprendentemente, los camellos salvajes no aportaron casi nada a la ascendencia de las poblaciones domésticas, ni siquiera a las poblaciones de Mongolia, que comparten hábitats cercanos con los camellos salvajes. Aunque la mayoría de los camellos bactrianos domésticos carecían de diferenciación cuando K crecía, IRÁN fue la primera población en separarse con una ascendencia única (K = 5; Fig. 2b).
Como otro método para examinar la estructura de la población con mezcla, inferimos el árbol de la población para los camellos utilizando TreeMix37 (Fig. 2c). Cuando no se añadió ninguna pista de migración, la topología del árbol volvió a indicar que IRÁN era la primera población separada entre todos los camellos bactrianos domésticos (Fig. 10 suplementaria). El aumento de las pistas de migración (m = 1-3) pudo mejorar en gran medida el ajuste del modelo (Fig. Suplementaria 11), que identificó los flujos de genes de los dromedarios a los camellos bactrianos domésticos en Asia Central, incluyendo KAZA, RUS e IRÁN con pesos de migración que oscilaban entre el 4% y el 9% (Tabla Suplementaria 8). Cabe mencionar que aunque la pista de migración apuntaba al final de la rama dromedaria hacia IRÁN, apuntaba al centro de la rama dromedaria hacia KAZA y RUS (Fig. 2c). Esto podría implicar una población fantasma relacionada con el dromedario iraní que contribuyó a la ascendencia de KAZA y RUS. Una pista de migración adicional (m = 4) podría seguir mejorando el ajuste del modelo, que indicaba la migración del dromedario a una raza XJ (Fig. 2c). Aunque TreeMix no detectó ninguna señal fuerte de migración entre los camellos bactrianos salvajes y los domésticos, los residuos mostraron una admixtura moderada entre las razas salvajes y las de Asia oriental (Fig. 11 suplementaria). A continuación, utilizamos la prueba de tres y cuatro poblaciones (F3/F4) menos paramétrica38 para evaluar la importancia estadística de estos eventos de mezcla. De nuevo, la prueba F3 apoyó fuertemente la mezcla de dromedarios y camellos bactrianos en KAZA, RUS e IRÁN (Tabla Suplementaria 9). La prueba F4, más sensible, confirmó un grado de mezcla significativamente mayor entre los dromedarios y los camellos bactrianos en Asia Central, en comparación con los de Asia Oriental (Tabla Suplementaria 10). Entre estos últimos, se detectó una mayor extensión de la admixtura con los dromedarios en XJ que en MG/IMG.
Evidencia para el origen de Asia Central eliminando la introgresión
Asia Oriental y Central fueron las dos regiones alternativas de domesticación para los camellos bactrianos basadas en evidencias arqueológicas1,12,17, pero la más probable quedó sin resolver. Aunque observamos la mayor diferenciación genética entre la población iraní y todas las demás domésticas, la existencia de mezcla entre dromedarios y camellos bactrianos en Asia Central debilitaría el apoyo a la inferencia de origen. Para reducir este efecto, intentamos eliminar los segmentos introgresados de los dromedarios de los genomas de los camellos bactrianos utilizando la prueba «BABA/ABBA «39. Agrupamos las poblaciones de Asia Oriental y Central, y comparamos el intercambio de alelos entre los dos grupos con los dromedarios (Fig. 3a). Como la ascendencia de los camellos bactrianos en un dromedario, así como la ascendencia de los camellos domésticos en tres salvajes (Fig. 2b), serían factores de confusión, eliminamos los cuatro individuos en el análisis. Utilizamos el estadístico fd, una versión robusta de la D de Patterson, para localizar los segmentos introgresados40 y aplicamos un nivel de significación estricto de Z-score = 2 mediante el procedimiento Jackknife (Fig. suplementaria 12). En un total de 21.153 segmentos de 100 kb no solapados a lo largo del genoma, había muchos más segmentos con señales putativas de introgresión en las poblaciones de Asia Central (11.711, Z-score > 2) que en las poblaciones de Asia Oriental (3891, Z-score < -2) como se esperaba. A continuación, realizamos el análisis de admixtura basado en los segmentos restantes y confirmamos que la introgresión de los dromedarios se redujo efectivamente (Fig. 13 suplementaria). El recálculo de la Fst por pares tras eliminar la introgresión seguía indicando que IRAN era la más diferenciada (0,04-0,06) entre todas las poblaciones domésticas (Fig. 3b). Para obtener más información sobre la filogenia de la población, reconstruimos el árbol NJ basado en el Fst por pares y realizamos la prueba bootstrap (Fig. 3b). Se confirmó que IRÁN fue la primera en separarse entre todas las poblaciones bactrianas domésticas, seguida de KAZA y RUS. El análisis bayesiano binario de Markov Chain Monte Carlo (MCMC) basado en la filogenia apoyó fuertemente a Asia Central como el área ancestral de los camellos bactrianos domésticos (probabilidad = 99,78%) y una ruta de dispersión posterior desde Asia Central a Asia Oriental (Fig. Suplementaria 14).
Un análisis BABA/ABBA para la introgresión de los dromedarios en los camellos bactrianos domésticos. Para centrarse en esta introgresión, se eliminó un dromedario con la ascendencia de los camellos bactrianos y tres camellos salvajes con la ascendencia de los domésticos. Se muestra el número de segmentos de 100 kb con fd significativo (|Z-score| > 2) para cada configuración de árbol. b Árbol de unión de vecinos (NJ) de las poblaciones después de eliminar los segmentos introgresados. El mapa de calor representa el promedio de Fst por pares para ventanas de deslizamiento de 5,1 × 104 10 kb. Los valores Bootstrap del árbol NJ se calcularon mediante el muestreo aleatorio de cinco mil ventanas de 10 kb durante 100 veces. c Árbol de máxima verosimilitud de ADNmt de longitud completa. Las poblaciones se representan con diferentes colores y las secuencias de Genbank se indican con puntos. Los valores de Bootstrap para las ramas principales están etiquetados.
Como evidencia independiente, también reconstruimos el árbol de máxima verosimilitud de ADNmt de longitud completa basado en las 128 muestras que secuenciamos en este estudio, así como 39 muestras adicionales disponibles en Genbank (Fig. 3c y Tabla Suplementaria 11). La introgresión de los ADNmt pudo ser fácilmente identificada y excluida del árbol. Por ejemplo, dos ADNmt recientemente secuenciados de KAZA y RUS se agruparon con dromedarios. Dentro del clado de los camellos bactrianos domésticos, aunque la mayoría de los camellos de diferentes regiones geográficas estaban mezclados, dos ADNmt de IRÁN formaban las ramas más basales de las poblaciones domésticas (Fig. 3c). El análisis bayesiano binario MCMC apoyó de nuevo el origen centroasiático de los camellos bactrianos domésticos (probabilidad = 76,43%).
Historia demográfica de los camellos bactrianos
Realizamos varios análisis de modelado paramétrico para inferir la dinámica demográfica de los camellos en la historia. En consonancia con un estudio anterior10, las trayectorias a largo plazo de los camellos bactrianos basadas en el modelo de coalescencia secuencial markoviana por pares (PSMC)41 revelaron una tremenda disminución del tamaño efectivo de la población de los camellos ancestrales hace más de un millón de años (Fig. 15 suplementaria). Aunque las trayectorias a largo plazo de los camellos bactrianos salvajes y domésticos fueron en general similares, resultaba obvio que divergían entre sí hace tan solo 0,4 millones de años, excluyendo a los primeros como progenitores directos de los segundos como los análisis previos de ADNmt9,14.
Para explorar el tiempo de divergencia entre las poblaciones de camellos con más detalle, utilizamos el muestreador coalescente filogenético generalizado (G-PhoCS)42. Dada la filogenia de las poblaciones de camellos, G-PhoCS pudo estimar el tamaño de la población a escala de mutaciones y el tiempo de divergencia de la población basándose en loci neutros no vinculados en genomas individuales de cada población. Para reducir la complejidad del modelo, sólo incluimos dromedarios, camellos bactrianos salvajes y tres poblaciones representativas (IRÁN, KAZA y MG) de camellos bactrianos domésticos (Fig. 16 y Tabla 12 suplementarias). Según la Fig. 3b, IRAN y KAZA fueron las dos primeras poblaciones de Asia Central que se separaron y la división de MG podría indicar la dispersión desde Asia Central a Asia Oriental. La edad se calibró asumiendo la divergencia entre el camello bactriano y el dromedario de 5,73 millones de años según la base de datos Timetree43. Cuando no se incorporó ninguna banda de migración, se logró fácilmente la convergencia de todas las estimaciones de los parámetros (Fig. 17 suplementaria y Tabla 13 suplementaria). De forma similar a los resultados del PSMC, el tamaño efectivo de la población disminuyó en general desde las poblaciones ancestrales hasta las modernas (Fig. 4). El tiempo de divergencia entre los camellos bactrianos salvajes y los domésticos se estimó en 0,43 millones de años (intervalo de confianza del 95%: 0,13-0,73 Mya; Fig. 4), ligeramente posterior al basado en el ADNmt (0,714 o 1,1 Mya9). Entre las poblaciones domésticas, IRÁN se separó de las demás hace unos 4,45 mil años (IC 95%: 0,07-17,6 Kya) y luego las poblaciones de Asia Central y Oriental se separaron hace unos 2,40 mil años (IC 95%: 0,01-7,84 Kya; Fig. 4).
El cambio en el tamaño efectivo de la población a escala de mutaciones θ está representado por colores de calor. El tiempo en años fue calibrado por el tiempo de divergencia entre dromedarios y camellos bactrianos. Los intervalos de confianza del 95% se muestran mediante barras en el eje temporal. La barra roja y azul indican la divergencia IRÁN-MG y KAZA-MG, respectivamente. Estas estimaciones se basan en el modelo sin migración.
Para permitir el flujo de genes, también intentamos introducir bandas de migración desde los dromedarios a las poblaciones de camellos bactrianos (Fig. Suplementaria 16 y Tabla Suplementaria 12). Las estimaciones sólo pudieron converger cuando se introdujo una banda de migración de los dromedarios iraníes a IRÁN y una banda de migración de una población fantasma a KAZA (Fig. 18 suplementaria y Tabla 13 suplementaria). Aunque el tiempo de divergencia entre los camellos bactrianos salvajes y los domésticos no se modificó con el modelo de migración (0,46 Mya, IC 95%: 0,24-0,71 Mya), el primer tiempo de divergencia de las poblaciones domésticas (0,19 Mya, IC 95%: 0,08-0,31 Mya) pasó a ser poco realista, porque estaba muy por encima de la historia conocida de la domesticación del ganado (11,5 Kya44). Además, la tasa de migración total sólo se estimó en un 0,27% y un 0,16% para la franja de migración a IRÁN y KAZA, respectivamente, mucho más baja que la estimada con Admixture (1-10%). Una posible razón para la mala estimación sería una historia de admixtura más compleja que el modelo de migración continua con tasas constantes asumido por G-PhoCS.