Muestreo

Se tomaron muestras de ovejas domésticas (O. aries) y cabras (C. hircus) en Irán (grupos IROA e IRCH, respectivamente) y Marruecos (grupos MOOA y MOCH, respectivamente) para un total de 20 animales por grupo (Fig. 6 suplementaria). Estas muestras se recogieron entre enero de 2008 y marzo de 2012 en la parte norte de Marruecos y entre agosto de 2011 y julio de 2012 en el noroeste de Irán, en el marco del proyecto europeo Nextgen (Acuerdo de subvención nº 244356) de acuerdo con las normas éticas de la Directiva 86/609/CEE de la Unión Europea. Se recogieron clips de oreja de la parte distal de la oreja de animales elegidos al azar, y se almacenaron inmediatamente en etanol al 96% durante un día antes de ser transferidos en perlas de gel de sílice hasta la extracción de ADN.

Las especies silvestres muflón asiático (O. orientalis) e íbice de Bezoar (C. aegagrus) fueron muestreadas en el noroeste de Irán dentro de la cuna de domesticación21,22. Se recogieron trece tejidos de muflones asiáticos y 18 de íbices de Bezoar (respectivamente, grupos IROO e IRCA, Fig. Suplementaria 6) de animales cautivos o recientemente cazados, y de muestras congeladas disponibles en el Departamento de Medio Ambiente iraní. Este enfoque de muestreo basado en el individuo está diseñado para minimizar el sesgo potencial evitando la sobrerrepresentación de los efectos locales (por ejemplo,

Datos adicionales

Además, se reunió un panel mundial de razas de ovejas y cabras (wpOA y wpCH, respectivamente). wpOA incluía 20 muestras de resecuenciación del genoma completo (WGS) con una cobertura de 12x que representaban 20 razas diferentes de todo el mundo proporcionadas por el Consorcio Internacional del Genoma Ovino. wpCH consistía en 14 muestras WGS secuenciadas con una cobertura de 12x que representaban a 9 individuos europeos, es decir, 2 Alpinos franceses y 2 Saanen franceses secuenciados por el INRA, 5 Saanen italianos proporcionados por el Parco Tecnologico Padano, y 5 individuos australianos, es decir 2 muestras de Boer, 2 de Rangeland y 1 de Cashmere proporcionadas por el CSIRO (Datos Suplementarios 5).

Producción de datos WGS

El ADN genómico se extrajo con éxito de todas las muestras de tejido utilizando el kit Macherey Nagel NucleoSpin 96 Tissue, adaptando el protocolo del fabricante. El muestreo de tejidos se realizó en bloques de pozos cuadrados de MN para obtener 25 mg de fragmentos por muestra. Se prepararon tres bloques cuadrados-96 de MN y la extracción se realizó con un manipulador de líquidos Tecan Freedom Evo siguiendo el protocolo del fabricante. Se realizó un paso de prelisis para homogeneizar las muestras con 180 µl de tampón T1 y 25 µl de proteinasa K durante la noche a 56 °C. Para ajustar las condiciones de unión, se añadieron 200 µl de tampón BQ1 y la placa de muestras se incubó 1 h a 70 °C; posteriormente se añadieron 200 µl de etanol al 100%. Los lisados se transfirieron a la placa de unión Nucleospin Tissue y se aplicó un vacío (-0,2 bar, 5 min) para eliminar el flow-through. Se realizaron tres pasos de lavado con los tampones BW y B5, respectivamente, y se aplicó de nuevo el vacío para descartar el flujo. Antes de la elución del ADN genómico, se secó una membrana de sílice de la placa de unión de tejidos Nucleospin al vacío con al menos -0,6 bares durante 10 minutos. El paso de elución se realizó con 100 µl de tampón BE precalentado (70 °C) y un paso de centrifugación a 3700 rcf durante 5 min en los pocillos de la 96-PCR. El ADN genómico se almacenó a 4 °C para evitar la congelación-descongelación y se comprobó su concentración (como ng/μl) mediante el método Picogreen y utilizando un Nanodrop.

Los genomas completos se resecuenciaron a partir de 500 ng de ADN genómico que se cortaron en un rango de 150-700 pb utilizando el instrumento Covaris® E210 para cada muestra y se utilizaron para la preparación de bibliotecas Illumina® mediante un protocolo semiautomatizado. La reparación de extremos, la cola A y la ligadura de adaptadores compatibles con Illumina® (BioScientific) se realizaron utilizando el sistema de preparación de bibliotecas SPRIWorks y el instrumento SPRI TE (Beckmann Coulter) siguiendo el protocolo del fabricante. Se aplicó una selección de tamaño de 300-600 pb para recuperar la mayoría de los fragmentos. Los fragmentos de ADN se amplificaron mediante 12 ciclos de PCR utilizando el kit de polimerasa Platinum Pfx Taq (Life® Technologies) y cebadores específicos para adaptadores de Illumina®. Las bibliotecas se purificaron con perlas AMPure XP 0,8x (Beckmann Coulter) y se analizaron con el bioanalizador Agilent 2100 (Agilent® Technologies) y la cuantificación qPCR. Las bibliotecas se secuenciaron utilizando la química de lecturas de 100 bases en celda de flujo de fin de par en el Illumina® HiSeq2000.

Las lecturas de fin de par de Illumina para Ovis se asignaron al genoma de referencia de oveja (OAR v3.1, ensamblaje GenBank GCA_000298735.146), y para Capra al genoma de referencia de cabra (CHIR v1.0, ensamblaje GenBank GCA_000317765.147), utilizando BWA-MEM48. El archivo BAM producido para cada individuo fue ordenado usando Picard SortSam y mejorado usando secuencialmente Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator y GATK IndelRealigner49, y SAMtools calmd50.

El descubrimiento de variantes se llevó a cabo usando tres algoritmos diferentes: Samtools mpileup50, GATK UnifiedGenotyper51, y Freebayes52. Los sitios de las variantes se identificaron de forma independiente para cada uno de los seis grupos, utilizando los modos multimuestra de los algoritmos de llamada: (i) 162 muestras de MOOA; (ii) 20 muestras de IROA; (iii) 14 muestras de IROO; (iv) 162 muestras de MOCH; (v) 20 muestras de IRCH; (vi) 19 muestras de IRCA. Para algunos grupos, se disponía de las WGS de más individuos como parte del proyecto NextGen (véase más arriba). Las muestras utilizadas en el presente estudio se seleccionaron para obtener grupos equilibrados de 20 individuos siempre que fue posible. Para los grupos IRCA e IROO, se dispuso de muestras adicionales en una fase posterior y se añadieron para los análisis posteriores. Se eliminaron los animales con baja calidad de alineación y llamadas para obtener el conjunto de datos final (Datos suplementarios 5).

Dentro de cada grupo, hubo dos rondas sucesivas de filtrado de la calidad del sitio de la variante. La etapa de filtrado 1 fusionó las llamadas de los tres algoritmos, filtrando al mismo tiempo las llamadas de menor confianza. Un sitio de variante pasó si fue llamado por al menos dos algoritmos de llamada diferentes con calidad de variante phred >30. Un alelo alternativo en un sitio pasó si fue llamado por cualquiera de los algoritmos de llamada, y el recuento de genotipos fue >0. La etapa de filtrado 2 utilizó la Recalibración de la Puntuación de Calidad de la Variante por GATK. En primer lugar, generamos un conjunto de entrenamiento de los sitios de variantes de mayor confianza dentro del grupo en el que (i) el sitio es llamado por los tres algoritmos de llamada de variantes con una calidad de variante escalada en fase >100, (ii) el sitio es bialélico, (iii) el recuento de alelos menores es al menos 3 mientras se cuentan sólo las muestras con una calidad escalada en fase de genotipo >30. El conjunto de entrenamiento se utilizó para construir un modelo gaussiano con la herramienta GATK VariantRecalibrator utilizando las siguientes anotaciones de variantes de UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Se aplicó un modelo gaussiano al conjunto de datos completo, generando un VQSLOD (log odds ratio de ser una variante verdadera). Los sitios se filtraron si VQSLOD <valor de corte. El valor de corte se estableció para cada grupo de la siguiente manera VQSLOD mínimo = {el valor medio de VQSLOD para las variantes del conjunto de entrenamiento}-3 × {la desviación absoluta media de VQSLOD de las variantes del conjunto de entrenamiento}. La relación transición/transversión SNP sugirió que el criterio de corte elegido daba el mejor equilibrio entre selectividad y sensibilidad.

Se generaron conjuntos de llamadas SNPs para seis grupos de animales Ovis y Capra (es decir, domésticos iraníes y marroquíes, y salvajes para cada género). Como los análisis realizados en este estudio requerían comparaciones entre grupos, creamos conjuntos de llamadas de genotipos en un conjunto consistente de sitios SNP para todos los animales de cualquier grupo. Para cada género, fusionamos los sitios de llamada de variantes de sus tres grupos, y sólo retuvimos las posiciones bialélicas sin datos perdidos. Los genotipos se volvieron a llamar en cada sitio SNP bialélico para todos los individuos de interés por GATK UnifiedGenotyper utilizando la opción GENOTYPE_GIVEN_ALLELES. En esta etapa, la lista de individuos se amplió para incluir los animales pertenecientes a los paneles de razas mundiales de ovejas y cabras (wpOA y wpCH) y las muestras silvestres adicionales que estuvieron disponibles en esta etapa (4 O. orientalis y 4 C. aegagrus). Los genotipos fueron mejorados y clasificados dentro de los grupos por Beagle 453, y luego filtrados cuando la probabilidad del genotipo era inferior a 0,95. Finalmente, filtramos los sitios que eran monomórficos a través de los diferentes subconjuntos de individuos utilizados en este estudio (ver más abajo).

Para comparar las señales de selección detectadas entre Ovis y Capra, realizamos un alineamiento cruzado entre los dos genomas de referencia. En primer lugar, utilizamos la línea de alineación por pares de la base de código Ensembl versión 6954 para alinear los genomas de referencia de oveja (OARv3.1) y cabra (CHIR1.0). Este proceso utiliza LastZ55 para alinear a nivel de ADN, seguido de un post-procesamiento en el que los bloques alineados se encadenan de acuerdo con su ubicación en ambos genomas. El proceso de alineación por pares LastZ es ejecutado rutinariamente por Ensembl para todas las especies soportadas, pero la cabra aún no está incluida en Ensembl. Para evitar el sesgo hacia cualquiera de las especies, produjimos dos alineaciones inter-específicas diferentes. Uno utilizó la oveja como genoma de referencia y la cabra como no referencia, mientras que el otro utilizó la cabra como genoma de referencia y la oveja como no referencia. La diferencia radica en que las regiones genómicas de la especie de referencia se ven obligadas a mapear de forma única a un único loci de la especie no de referencia, mientras que a las regiones genómicas no de referencia se les permite mapear a múltiples localizaciones de la especie de referencia. Obtuvimos para segmentos de cromosomas de un genoma de referencia las coordenadas en el genoma de no referencia. Finalmente, para los SNPs descubiertos en un género, utilizamos la alineación del genoma completo con el genoma de referencia del otro género para identificar las posiciones correspondientes (Tabla Suplementaria 6).

Estructura genética

Para describir la diversidad genética dentro de los grupos, utilizamos VCFtools56 para calcular las estadísticas de resumen de la variación genética en los 73 individuos para Ovis (es decir, 13 IROO, 20 IROA, 20 MOOA, y 20 wpOA) y 72 individuos para Capra (es decir, 18 IRCA, 20 IRCH, 20 MOCH, y 14 wpCH). Los estadísticos medidos fueron el número total de variantes polimórficas (S) para el conjunto de individuos de cada género y dentro de cada grupo, la diversidad nucleotídica media (π) dentro de cada grupo y el coeficiente de endogamia (F) para cada individuo. Dentro de cada género, las diferencias entre el grupo silvestre y cada grupo doméstico se comprobaron mediante una prueba t unilateral para los valores individuales de endogamia y carga genética, y una prueba Mann-Whitney bilateral para la diversidad de nucleótidos por sitio.

La divergencia global entre los cuatro grupos dentro de cada género (es decir, silvestre, domésticos iraníes y marroquíes, y panel mundial) se estimó utilizando todos los SNPs bialélicos y el promedio ponderado de Fst por pares siguiendo a Weir y Cockerham57 según lo implementado en VCFtools56. La estructura genética entre los grupos se evaluó con el método de agrupación sNMF26, después de podar el conjunto de datos para eliminar los SNP con desequilibrio de enlace (r²) superior a 0,2 utilizando VCFtools. El desequilibrio de ligamiento (r²) se calculó entre pares de SNPs dentro de ventanas deslizantes de 50 SNPs, con un SNP por par eliminado aleatoriamente cuando r² era superior a 0,2. Para cada análisis sNMF, se realizaron cinco ejecuciones del mismo número de clusters (K) con valores de K de 1 a 10. Se utilizó el criterio de entropía cruzada para identificar la solución de agrupación más probable, sin embargo, también se exploraron particiones alternativas para diferentes números de K con el fin de evaluar cómo se dividían los individuos entre los clústeres.

Para desentrañar entre la ascendencia compartida y la mezcla, ejecutamos TreeMix27 para estimar conjuntamente las divisiones de la población y los subsiguientes eventos de mezcla utilizando el conjunto de datos podados utilizados para sNMF. Ejecutamos TreeMix con la opción -global para refinar nuestras inferencias de máxima probabilidad. Enraizamos el árbol TreeMix con la división entre individuos salvajes y domésticos. El tamaño del bloque para el jackknifing fue de -k 500 SNPs, que corresponde aproximadamente a 150 kb, superando los bloques medios de LD encontrados tanto en ovejas como en cabras. Generamos un árbol de Máxima Verosimilitud sin migración y luego añadimos eventos de migración y examinamos el cambio incremental en la varianza explicada por el modelo y los valores residuales entre individuos. El objetivo era detectar cualquier valor residual potencialmente alto o borde de migración entre individuos salvajes y domésticos. Para explorar aún más la relevancia estadística de los posibles vectores de mezcla identificados por TreeMix (Tabla Suplementaria 3), calculamos la prueba de tres poblaciones f328 como una prueba formal de introgresión genética, utilizando el programa qp3Pop de la suite ADMIXTOOLS58 para cada combinación de grupos. Para Capra, el grupo wpCH se dividió entre las razas australianas, francesas e italianas. Los resultados se reportan en los Datos Suplementarios 2.

Inferencia demográfica

Para cada género, realizamos análisis de inferencia demográfica ancestral utilizando el modelo MSMC implementado en el software MSMC225. MSMC se basa en la coalescencia secuencial markoviana por pares59; sin embargo, utiliza como entrada los haplotipos de los datos de la secuencia del genoma por fases. Para cada análisis utilizamos dos individuos de un grupo, es decir, 4 haplotipos. Cada análisis se repitió para otro conjunto aleatorio de dos individuos, es decir, una réplica del análisis por grupo. Los archivos de entrada y salida se generaron y analizaron con los scripts de python proporcionados con el software MSMC y que se encuentran en https://github.com/stschiff/msmc-tools. Los parámetros del análisis se mantuvieron por defecto, excepto la tasa de mutación que se fijó en 2,5×10-8 y la longitud de la generación se fijó en 2 años. Para estimar la incertidumbre en las estimaciones de tiempo, variamos estos parámetros (tasa de mutación de 2,5×10-8 y 1,0×10-8 en combinación con la longitud de la generación de 2 y 4 años) y proporcionamos una estimación aproximada del período de domesticación (véase la Fig. 2 suplementaria).

Carga genética

La carga genética se estimó de dos maneras. En primer lugar, calculando la carga genética para cada individuo como la suma de los efectos deletéreos de fitness sobre todas las posiciones genómicas que codifican proteínas siguiendo el método de Librado et al.60. Brevemente, como un proxy para la restricción evolutiva, utilizamos las puntuaciones PhyloP del alineamiento de 46 mamíferos (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). A partir de este alineamiento, identificamos los sitios de codificación de proteínas que evolucionan bajo restricciones funcionales (puntuación PhyloP ≥1,5). Para cada genoma de Ovis o Capra, investigamos entonces si estos sitios estaban mutados. En caso afirmativo, sumamos las puntuaciones phyloP de todos los sitios mutados, de modo que las mutaciones en sitios altamente restringidos contribuyen proporcionalmente más a la estimación de la carga total. Esto proporcionó una estimación de la carga para cada genoma de oveja/cabra. Finalmente, para obtener una carga media por sitio, la dividimos por el número total de posiciones analizadas. Cabe destacar que condicionamos los sitios homocigotos para evitar modelar el coeficiente de dominancia de las mutaciones en sitios heterocigotos (por ejemplo, recesivo, intermedio, dominante). En segundo lugar, comparamos gen a gen la carga deletérea genética en los grupos de Ovis silvestres y domesticados realizando una prueba de Wilcoxon, siendo la hipótesis alternativa que los animales domésticos tienen más carga que los parientes silvestres. Los valores p fueron corregidos para pruebas múltiples61 y aplicamos un umbral de valores p ajustados < 0,05. Realizamos un análisis de enriquecimiento de la ontología genética en el conjunto de genes que mostraban un aumento significativo de la carga genética utilizando WebGestalt62,63. Como los genomas de referencia están mal anotados para los genes, nos basamos en los ortólogos de una sola copia entre nuestra especie y el ser humano y el ratón. Los genes del cromosoma X se excluyeron del conjunto de antecedentes. No llevamos a cabo este análisis en Capra debido a la mayor endogamia observada en las muestras silvestres.

Detección de firmas de selección

Para detectar firmas de selección relacionadas con la domesticación, utilizamos todos los SNPs bialélicos que mostraban una frecuencia alélica menor superior a 0,10 en al menos uno de los tres grupos analizados (es decir, los grupos domésticos iraní y marroquí, y el grupo silvestre para cada género). Dado que esperábamos que las firmas de selección relacionadas con el proceso de domesticación estuvieran presentes en todos los animales domésticos, adoptamos la siguiente estrategia general: probamos con hapFLK29 (véase la Nota Suplementaria 5 y las Figuras Suplementarias 9, 10 y 11) para cada género el grupo silvestre frente a cada uno de los grupos domésticos manejados tradicionalmente (es decir, el iraní y el marroquí) y nos centramos en aquellas regiones comunes putativamente bajo selección que se detectaron en ambos casos. El tamaño de las muestras de los grupos (n = 13-20) era compatible con los requisitos del método29. Comprobamos visualmente si las señales consistentes de selección encontradas con hapFLK estaban también presentes en el correspondiente conjunto de paneles mundiales de cada género, pero no incluimos estos grupos en la prueba estadística debido a su composición multirracial. Por último, buscamos señales de selección compartidas entre Ovis y Capra utilizando un enfoque FDR estratificado. La estrategia se representa en la Fig. Suplementaria 4.

Realizamos pruebas hapFLK para contrastar el grupo salvaje con cada uno de los grupos iraníes y marroquíes de cada género. La matriz de parentesco se calculó a partir de las distancias genéticas de Reynold64 entre pares de grupos, utilizando un subconjunto aleatorio del uno por ciento de las variantes. El árbol poblacional inferido se construyó utilizando el algoritmo de unión de vecinos. Para cada SNP, realizamos la prueba hapFLK que incorpora información haplotípica para aumentar el poder de detección de barridos selectivos. Para cada SNP probado, el estadístico hapFLK calculó la desviación de las frecuencias haplotípicas con respecto al modelo neutral estimado por la matriz de parentesco65. Para explotar la información de desequilibrio de ligamiento, hapFLK utiliza el modelo multipunto de Scheet y Stephens’66 para genotipos multilocus que puede ajustarse a datos sin fase. Una de las principales aplicaciones de este modelo es realizar la estimación de fase (software fastPHASE66). En nuestro análisis, el modelo fue entrenado con datos sin fase, y por lo tanto nuestro análisis tiene en cuenta la incertidumbre de la fase. El método se utilizó para reagrupar los haplotipos locales a lo largo de los cromosomas en un número especificado de clusters K fijado en 25, utilizando un modelo de Markov oculto.

Para identificar las regiones comunes putativamente bajo selección en los dos grupos domésticos tradicionalmente manejados para cada género, combinamos los dos análisis hapFLK anteriores. Para cada análisis las puntuaciones hapFLK se ajustaron a una distribución χ2 para obtener los valores p (script disponible en https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Los resultados de los dos contrastes entre el grupo salvaje y cada uno de los grupos domésticos se combinaron utilizando el método de Stouffer67 para obtener valores p únicos para la comparación de los animales salvajes frente a los domésticos. Por último, se aplicó el marco FDR68 a todo el conjunto de SNPs para convertir los valores p combinados en valores q. Los SNPs que mostraban valores q < 10-2 fueron retenidos y agrupados en regiones genómicas cuando estaban a menos de 50 kb de distancia unos de otros.

Para investigar si la señal de selección era compartida entre Ovis y Capra, primero utilizamos la alineación cruzada de los dos genomas de referencia para identificar segmentos homólogos. A continuación, aplicamos un marco FDR estratificado69. Este enfoque se basa en el hecho de que existe una estratificación inherente en las pruebas dada la información previa en los datos genéticos69, porque la distribución subyacente de las verdaderas hipótesis alternativas podría ser diferente según las diferentes dinámicas de varias regiones genómicas, lo que lleva a diferentes distribuciones de valores p. Esto requiere obtener valores p ajustados por FDR (es decir, valores q) por separado para los diferentes estratos. Buscamos convergencias en cada género separando las regiones homólogas a las detectadas en el otro género (denominadas estrato compartido) y el resto del genoma (denominado estrato general). Extrajimos los valores p por separado para cada uno de los dos estratos definidos y luego calculamos los valores q mediante el marco FDR. Estos valores q estratificados fueron las cantidades finales consideradas para la significación estadística (<10-2) para detectar los SNP bajo selección y fusionarlos en las regiones genómicas correspondientes.

Para probar las firmas convergentes de selección que diferencian a los animales salvajes de los domésticos en ambos géneros, examinamos la relación entre el umbral de significación aplicado a los valores q (que hicimos variar de 0,2 a 0,002) en un género y la probabilidad estimada de que un SNP sea seleccionado en el estrato compartido del otro género utilizando el enfoque de Storey et al.70 …y la probabilidad estimada de que un SNP sea seleccionado en el estrato compartido del otro género. Un aumento de la probabilidad inferida con una disminución del umbral aplicado al valor q (aumento de la rigurosidad) indica que cuanto más significativa es la región en un género, más probable es que encontremos SNPs significativos en el otro género.

Filtramos las señales de selección que no eran consistentes entre los tres grupos domésticos. Para cada región detectada, utilizamos los haplotipos en fase de los individuos que se agruparon utilizando árboles Neighbor-Joining basados en el porcentaje de identidad entre las secuencias. Sólo se mantuvieron las regiones que mostraban señales consistentes (Fig. Suplementaria 5).

Para inferir si las señales de selección detectadas con hapFLK indicaban relajación de la selección o selección positiva en los domésticos, estimamos la diferencia en la diversidad de nucleótidos (π) en cada región putativa bajo selección entre los grupos salvajes y domésticos. Expresamos esta diferencia como el índice Δπ, que se calculó para cada región genómica como la diferencia entre π calculada para el grupo silvestre y la media de π para los grupos domésticos iraní y marroquí, menos la diferencia de π entre estos dos grupos calculada sobre todo el genoma: