mai 20, 2020 Par sexe2 0

L’origine des gènes de domestication chez les chèvres

Abstrait

La domestication des chèvres a été fondamentale pour l’agriculture et la civilisation, mais ses changements génétiques et ses régimes d’élevage restent flous. Ici, nous analysons les génomes des chèvres domestiques à travers le monde, les espèces de chèvres sauvages et les vestiges historiques, fournissant la preuve d’un ancien événement d’introgression d’une espèce turque dans le Caucase occidental à l’ancêtre des chèvres domestiques. Un lieu introgresso avec une signature forte de la sélection accueille le MUC6 gène, qui code pour une mucine sécrétée par le tractus gastro-intestinal. Des expériences ont révélé que l’haplotype introgressé presque fixe confère une plus grande résistance immunitaire aux pathogènes gastro-intestinaux. Un autre locus avec un signal de sélection fort peut être lié au comportement. Les allèles sélectionnés dans ces deux loci ont émergé chez les chèvres domestiques il y a au moins 7200 et 8100 ans, respectivement, et ont augmenté à des fréquences élevées en conjonction avec l’expansion de l’haplogroupe mitochondrial moderne omniprésent A. Le suivi de ces transformations évolutionnaires archéologiquement cryptiques fournit de nouvelles perspectives sur les mécanismes de domestication des animaux.

INTRODUCTION

La domestication présente un changement extrême du stress physiologique et comportemental chez les animaux en liberté (1) dans un environnement anthropique à haute densité et propice aux maladies, en particulier pour les herbivores. La chèvre (Hircus de chèvre) a été l’une des premières espèces animales domestiques domestiquées, démontrant une adaptabilité et une polyvalence remarquables (2, 3) et a été intimement associée à la dispersion humaine (4). Des études récentes ont identifié des candidats cibles à sélectionner lors de la domestication des chèvres, y compris des loci liés à la pigmentation, au métabolisme xénobiotique et à la production de lait (5, 6); cependant, la dynamique évolutive des gènes clés impliqués dans l’adaptation pendant la phase initiale de domestication reste incertaine.

La domestication des chèvres aurait commencé il y a environ 11 000 ans à partir d’une mosaïque de populations de bézoards sauvages (Chèvre aegagrus) (2, 6). Cependant, d’autres sauvages Chèvre les espèces sont largement réparties et nombre d’entre elles peuvent s’hybrider avec des chèvres domestiques (79). Leur contribution au processus de domestication des chèvres reste inexplorée. L’introgression des variantes adaptatives a été largement reconnue comme un phénomène évolutif significatif chez l’homme (dix), mouton (11) et le bétail (12) et peut être particulièrement efficace pour augmenter la forme physique sans effets pléiotropes négatifs, comme l’ont démontré d’autres espèces (13). Ici, nous avons mené une étude génomique complète sur la population des populations de chèvres modernes réparties dans le monde, six Chèvre , et des génomes de chèvres antiques publiés précédemment pour étudier les changements temporels dans les principales variantes génétiques sélectionnées pendant la domestication des chèvres.

RÉSULTATS

Structure de la population et origine des chèvres domestiques

Nous avons séquencé 101 Chèvre génomes (couverture 3 à 47 fois; moyenne, 12 fois), dont 88 chèvres domestiques récoltées sur trois continents différents, un bézoard, un bouquetin des Alpes (Chèvre bouquetin), trois bouquetins de Sibérie (Chèvre Sibirica), trois Markhors (Falconeri de chèvre), un Turc du Caucase occidental (Chèvre du Caucase) et quatre bouquetins de Nubie × hybrides de chèvres domestiques (Chèvre nubienne × C. hircus) (Fig. 1A et texte S1). Nous avons également séquencé cinq échantillons de chèvres anciennes avec une couverture d’environ 0,04 à 13,44 fois (figure 1A et texte S1), y compris les premiers échantillons archéologiques chinois connus de la période néolithique (14). Avec des séquences génomiques accessibles au public pour la chèvre et le bézoard modernes (tableau S1) (5, 15), ainsi que les anciens génomes de (tableau S4) (6), nous avons collecté une collection mondiale de 164 chèvres domestiques modernes, 52 chèvres anciennes, 24 bézoards modernes et 4 bézoards anciens.

Fig. 1 Répartition géographique et affinités génétiques des chèvres sauvages et domestiques utilisées dans cette étude.

(UNE) Positions des chèvres sauvages (carrés), domestiques (cercles) et antiques (pentagone) pour chaque groupe géographique principal. Les chèvres domestiques sont colorées pour refléter leur origine géographique. BP, avant le présent. (B) Un arbre qui unit les voisins des séquences du génome utilisées dans cette étude. Les branches sont colorées suivant le même code couleur utilisé en (A). (C est ) APC de bézoards et de chèvres domestiques (C) et uniquement de chèvres domestiques (D). (EST) ADMIXTURE résultats pour K = 3 e K = 6, avec l’erreur de validation croisée faible.

« data-icon-position = » « data-hide-link-title = » 0 « >

Fig. 1 Répartition géographique et affinités génétiques des chèvres sauvages et domestiques utilisées dans cette étude.

(UNE) Positions des chèvres sauvages (carrés), domestiques (cercles) et antiques (pentagone) pour chaque groupe géographique principal. Les chèvres domestiques sont colorées pour refléter leur origine géographique. BP, avant le présent. (B) Un arbre qui unit les voisins des séquences du génome utilisées dans cette étude. Les branches sont colorées suivant le même code couleur utilisé en (A). (C est ) APC de bézoards et de chèvres domestiques (C) et uniquement de chèvres domestiques (D). (EST) ADMIXTURE résultats pour K = 3 e K = 6, avec l’erreur de validation croisée faible.

Un arbre phylogénétique qui rejoint le génome voisin entier révèle que toutes les chèvres domestiques forment une lignée monophylétique sœur vers les bézoards (Fig.1B et fig.S4), confirmant que les chèvres domestiques modernes sont les descendants des ancêtres bézoards (16). Les quatre autres sauvages Chèvre espèces (C. bouquetin, C. sibirica, C. falconeri, est C. caucasica), appelées espèces de type bouquetins (17), tombent exclusivement dans une autre branche divergente de la chèvre bezoaria (Fig. 1B et Fig. S4). L’analyse en composantes principales (ACP) montre que les trois populations de bézoards, qui ont été collectées au Moyen-Orient (figure 1A) près du centre de domestication (2, 18), sont structurées et regroupées correspondant à leurs origines géographiques (Fig. 1C). Au sein des chèvres domestiques, l’ACP et le clustering basé sur des modèles (K = 3) montre que les chèvres asiatiques sont génétiquement distinctes des échantillons européens (EUR) et africains (AFR) (Fig. 1, D et E). À K = 6, les chèvres d’Asie sont divisées en deux sous-groupes géographiques: Asie du Sud-Ouest-Asie du Sud (SWA-SAS) et Asie de l’Est (EAS) (figure 1E). Cette structure géographique est également soutenue par TreeMix (fig. S8) et des statistiques basées sur l’haplotype (fig. S9), conformément au scénario selon lequel les ancêtres des populations de chèvres domestiques d’aujourd’hui ont suivi des voies distinctes de dispersion le long de la axe est-ouest de l’Afro-Eurasie (figure 2A et tableau S6) (4). Ce scénario de dispersion a été confirmé par l’observation que les populations éloignées du centre de domestication (non SWA) présentaient un déséquilibre de connexion élevé et une diversité génétique réduite (fig. S10).

Fig.2 Flux de gènes dans la phase initiale de domestication et de propagation de la chèvre.

(UNE) Répartition géographique de la nature Chèvre espèces et voies de dispersion des chèvres domestiques en dehors de leurs zones de domestication (3, 7). Les positions d’échantillonnage des espèces semblables à celles des bouquetins sont indiquées par des carrés. (B) Partage des allèles entre les bézoards modernes (H3) et les chèvres domestiques. (C) Partage des allèles entre les chèvres domestiques (H3) et les espèces semblables à des bouquetins (H1). Voir également le tableau S8 pour les autres populations (H3) étudiées. Résultats statistiquement significatifs, définis par |Z scores | ≥ 3, sont marqués d’un astérisque rouge pour (B) et (C). () Un diagramme de dispersion de la fréquence des haplotypes introgressifs dans les populations de chèvres EUR-AFR et asiatiques. Les points rouges représentent les loci liés au système immunitaire. (EST) Une carte thermique de des tests statistiques d’affinité différentielle entre SWA (bleu) et H2 (rouge), dans lesquels H2 représente l’individu ou la population de bézoards et de chèvres domestiques situés sur la carte. Les chèvres d’Asie de l’Est sont représentées dans l’encadré dans le coin supérieur droit. (F) Réseau d’haplotypes du nombre de différences par paires à MUC6 région non répétée. Toutes les chèvres domestiques sont de couleur verte et sauvage Chèvre les espèces sont divisées en cinq sous-groupes, y compris les haplotypes d’hybrides bouquetins nubiens × chèvres domestiques. Le rayon du graphique à secteurs et la largeur des bords étaient log2-transformé.

« data-icon-position = » « data-hide-link-title = » 0 « >

Fig.2 Flux de gènes dans la phase initiale de domestication et de propagation de la chèvre.

(UNE) Répartition géographique de la nature Chèvre espèces et voies de dispersion des chèvres domestiques en dehors de leurs zones de domestication (3, 7). Les positions d’échantillonnage des espèces semblables à celles des bouquetins sont indiquées par des carrés. (B) Partage des allèles entre les bézoards modernes (H3) et les chèvres domestiques. (C) Partage des allèles entre les chèvres domestiques (H3) et les espèces semblables à des bouquetins (H1). Voir également le tableau S8 pour les autres populations (H3) étudiées. Résultats statistiquement significatifs, définis par |Z scores | ≥ 3, sont marqués d’un astérisque rouge pour (B) et (C). () Un diagramme de dispersion de la fréquence des haplotypes introgressifs dans les populations de chèvres EUR-AFR et asiatiques. Les points rouges représentent les loci liés au système immunitaire. (EST) Une carte thermique de des tests statistiques d’affinité différentielle entre SWA (bleu) et H2 (rouge), dans lesquels H2 représente l’individu ou la population de bézoards et de chèvres domestiques situés sur la carte. Les chèvres d’Asie de l’Est sont représentées dans l’encadré dans le coin supérieur droit. (F) Réseau d’haplotypes du nombre de différences par paires à MUC6 région non répétée. Toutes les chèvres domestiques sont de couleur verte et sauvage Chèvre les espèces sont divisées en cinq sous-groupes, y compris les haplotypes d’hybrides bouquetins nubiens × chèvres domestiques. Le rayon du graphique à secteurs et la largeur des bords étaient log2-transformé.

Nos analyses démographiques utilisent plusieurs coalescents Markoviens séquentiellement (MSMC), SMC ++ et ∂uneje indiquent que les périodes de divergence entre les populations modernes de chèvres asiatiques et européennes ont précédé le temps estimé sur le plan archéologique de domestication (il y a environ 11 000 ans) (fig. S12 et texte S2). Aussi, les statistiques ont révélé que les bézoards de Zagros montrent une plus grande affinité génétique avec les populations domestiques orientales, tandis que les bézoards d’Azerbaïdjan montrent une plus grande affinité avec les populations domestiques occidentales (Fig. 2B et fig. S15A). Par conséquent, les temps de coalescence profonde enregistrés entre les paires de populations est-ouest peuvent s’expliquer par leur origine dans les populations de bézoards ancestraux structurés (6) ou par recrutement post-domestication dans différentes populations de bézoards locaux.

Flux de gènes d’espèces semblables à des bouquetins aux bézoards prédomestifs et aux chèvres modernes

les statistiques révèlent que les quatre espèces semblables à des bouquetins ont des signes significatifs de partage de l’allèle avec des chèvres anciennes et modernes, indiquant un mélange (figure 2C et tableau S8). Nous avons donc examiné ce modèle génétique de mélange entre les espèces de bouquetins et les chèvres domestiques en utilisant statistiques et identité par statut dans des fenêtres coulissantes de 20 ko. Nous avons également vérifié les régions introgressives candidates en utilisant des arbres phylogénétiques Sprime et le maximum de vraisemblance (ML). En utilisant un critère conservateur (c’est-à-dire en ne conservant que des haplotypes introgressifs putatifs avec une fréquence supérieure à 0,1 chez les chèvres), nous avons identifié 112 segments génomiques superposés à 81 gènes codant pour des protéines avec des signatures d’introgression d’espèces semblables à des bouquetins (Fig. 2D, fig. S16 et fichiers de données S1). Une analyse d’enrichissement de la route des gènes et génomes de Kyoto (KEGG) pour ces gènes montre que la catégorie la plus enrichie est l’amibiase (test hypergéométrique adapté P <5,28 × 10-3: Tableau S10), qui est lié à l’invasion parasitaire et à l’immunosuppression, y compris quatre gènes (SERPINB3, SERPINB4, CD1b, est COL4A4). Trois gènes supplémentaires (BPI, MAN2A1, est CD2AP) participent également à la fonction immunitaire (1921). Dans ces segments, nous avons observé une prononciation prononcée des allèles introduits putativement par le turbo du Caucase occidental (fig. S17), cohérente avec cette espèce qui montre le plus grand partage d’allèles du génome avec les chèvres domestiques (Fig. 2C).

En étudiant l’histoire du mélange turc du Caucase de l’Ouest, nous avons constaté que le tur du Caucase de l’Ouest partageait plusieurs allèles avec les quatre bézoards pré-domestiqués (un bézoard arménien datant de> 47 000 ans et trois bézoards anatoliens datant d’environ 13 000 ans). et en particulier le bézoard arménien, par rapport aux Bezoars d’aujourd’hui et aux animaux domestiques (Fig. 2E). Plus loin F3 statistiques sous forme de F3 (bezoar moderne, turbo de l’ouest du Caucase; cible) soutient le flux de gènes du turbo de l’ouest du Caucase vers l’ancien bézoar arménien (tableau S9). Avec le fait que trois Bezoars d’Anatolie pré-domestication avaient un haplotype mitochondrial semblable à un tur6), nos résultats fournissent la preuve d’un événement de mélange pré-domestication de bézoard avec une espèce turbulente. Cet ancien événement d’addition est également confirmé par les analyses TreeMix (fig. S8C) et par la proximité géographique étroite entre le tur du Caucase occidental et l’ancien bézoard arménien, les deux étant à leur tour plus proches du centre de domestication des chèvres. que toute autre espèce de bouquetin vivant similaire (figure 2A). Bien que tous les génomes de chèvre modernes suggèrent une combinaison avec tur, les chèvres néolithiques des Balkans et les chèvres modernes EUR et AFR présentent un plus grand partage d’allèles avec tur que les chèvres asiatiques (fig. S15B), probablement en raison de leur plus grande affinité génétique avec les anciens Arméniens et Anatoliens (Fig.2E) (6). Par conséquent, l’introgression de la pré-domestication d’une espèce similaire à la turque peut avoir contribué des allèles aux anciens bézoards et donc aux premières chèvres domestiques et à leurs populations dérivées modernes.

En particulier, une région introgressée sur le chromosome 29 a un haplotype introgressif presque fixe (fréquence = 95,7%) chez les chèvres domestiques modernes partout dans le monde (Fig. 2F et Fig. S22A). Cette région contient un gène codant pour les protéines intact, MUC6, qui code pour une mucine sécrétée par le tractus gastro-intestinal qui forme un revêtement protecteur de glycoprotéine impliquée dans les réponses immunitaires innées de l’hôte à l’invasion de multiples pathogènes gastro-intestinaux (2224). Nous avons recherché des populations de donneurs potentiels dans nos génomes de caprins sauvages séquencés. Le réseau haplotypique construit avec la région non répétée de MUC6 montre que l’haplotype domestique le plus courant (MUC6) est similaire au tur du Caucase occidental, différant par une seule mutation, bien qu’il soit différent des haplotypes de fréquence plus faible des chèvres domestiques et de l’haplotype commun chez le bezoar (MUC6B) (Fig.2G). Pour valider davantage ce signal d’introgression, nous avons estimé l’ancêtre commun le plus récent des haplotypes divergents et étudié si la sélection qui agit sur une mutation de novo ou sur une variation permanente pourrait éventuellement conduire au modèle dans cette région chez les chèvres domestiques. Les estimations du temps de coalescence (fig. S18A) et les simulations neutres (fig. S18B) suggèrent que le modèle observé de différenciation des haplotypes est hautement improbable en l’absence de flux de gènes interspécifiques. Par conséquent, le presque résolu MUC6 chez les chèvres, il a probablement été introgressé par une descente près du Turc du Caucase occidental, conformément au signal de l’ensemble du génome.

Sélection médiée par la bonne des gènes immunitaires et neuronaux

Pour identifier les balayages sélectifs clés dans la domestication des chèvres, nous avons comparé les populations de chèvres domestiques du monde entier avec les 24 bézoards en estimant la différenciation génétique par paires (FST), les différences de diversité nucléotidique (rapport π bézoard / domestique) et l’homozygotie de l’haplotype étendu à la population (XP-EHH) dans des fenêtres glissantes de 50 kb le long du génome (fig. S19 et S20). Nous avons défini des fenêtres avec des valeurs significatives (Z tester, P <0,005) dans les trois statistiques (FST > 0,195, rapport π ln> 0,395 et XP-EHH> 2,10) en tant que balayages sélectifs putatifs. Après avoir rejoint des fenêtres anormales consécutives, 105 régions uniques contenant 403 gènes codant pour les protéines ont été identifiées (fichiers de données des figures 3A et S2). Dix-huit de ces régions ont été identifiées lors d’analyses à domicile précédentes, y compris celles associées aux effets phénotypiques liés à l’immunité, aux voies ou processus neuronaux, à la pigmentation et aux caractéristiques de productivité associées à la composition du lait et aux caractéristiques des cheveux (fichier de données S2). Le chevauchement modeste avec les études précédentes est principalement dû aux différences dans les ensembles d’échantillons, dans la méthodologie de balayage sélectif et dans les différentes versions du génome de référence (testo S3).

0,195, rapport π ln> 0,395 et XP-EHH> 2,1) utilisé pour l’extraction des valeurs aberrantes. (B) Voies KEGG identifiées comme significativement surreprésentées (test hypergéométrique, P ajusté <0,01). (C et D) Distribution sélective de balayage et de recombinaison (RHO) sur le chromosome 29 (46,22 à 46,31 Mo) (C) et balayage sélectif sur le chromosome 15 (32,24 à 32,37 Mo) (D). La région de balayage putatif est également validée par les tests FST, Tajima D et CLR. Les lignes pointillées horizontales représentent la moyenne de l'ensemble du génome pour les paramètres correspondants. Les annotations des gènes dans la région de balayage et SNP presque fixe pour les allèles dérivés chez les chèvres domestiques (C et D) sont indiquées ci-dessous. "Class =" fragment-images colorbox-load "rel =" gallery-fragment-images-1140014933 "data -the figure-caption ="
Fig. 3 Régions génomiques avec signaux de sélection chez les chèvres domestiques.

(UNE) Répartition de l’indice de fixation par paires (FST) (X axe), rapport π ln (y axe) et la valeur XP-EHH (couleur) entre les bézoards et les chèvres domestiques. Les lignes verticales et horizontales en pointillés indiquent le seuil de signification (correspondant à Z tester, P <0,005, où FST > 0,195, rapport π ln> 0,395 et XP-EHH> 2,1) utilisé pour l’extraction des valeurs aberrantes. (B) Voies KEGG identifiées comme significativement surreprésentées (test hypergéométrique, ajusté P <0,01). (C est ) Distribution sélective de balayage et de recombinaison (RHO) sur le chromosome 29 (46,22 à 46,31 Mo) (C) et balayage sélectif sur le chromosome 15 (32,24 à 32,37 Mo) (D). La région de balayage putatif est également validée par FST, De Tajima et tests CLR. Les lignes pointillées horizontales représentent la moyenne de l’ensemble du génome pour les paramètres correspondants. Les annotations des gènes dans la région de balayage et SNP presque fixe pour les allèles dérivés chez les chèvres domestiques (C et D) sont indiquées ci-dessous.

« data-icon-position = » « data-hide-link-title = » 0 « >

Fig. 3 Régions génomiques avec signaux de sélection chez les chèvres domestiques.

(UNE) Répartition de l’indice de fixation par paires (FST) (X axe), rapport π ln (y axe) et la valeur XP-EHH (couleur) entre les bézoards et les chèvres domestiques. Les lignes verticales et horizontales en pointillés indiquent le seuil de signification (correspondant à Z tester, P <0,005, où FST > 0,195, rapport π ln> 0,395 et XP-EHH> 2,1) utilisé pour l’extraction des valeurs aberrantes. (B) Voies KEGG identifiées comme significativement surreprésentées (test hypergéométrique, ajusté P <0,01). (C est ) Distribution sélective de balayage et de recombinaison (RHO) sur le chromosome 29 (46,22 à 46,31 Mo) (C) et balayage sélectif sur le chromosome 15 (32,24 à 32,37 Mo) (D). La région de balayage putatif est également validée par FST, De Tajima et tests CLR. Les lignes pointillées horizontales représentent la moyenne de l’ensemble du génome pour les paramètres correspondants. Les annotations des gènes dans la région de balayage et SNP presque fixe pour les allèles dérivés chez les chèvres domestiques (C et D) sont indiquées ci-dessous.

Une analyse KEGG a montré que 9 des 14 voies significativement enrichies (test hypergéométrique, ajusté P <0,01) sont d'origine immunologique (test hypergéométrique, ajusté P = 5 × 10-3 à 2,35 × 10-7) (Figure 3B et tableau S12), y compris la grippe A, le paludisme, la trypanosomiase africaine, la cytotoxicité à médiation cellulaire naturelle, la rougeole, l’infection par l’herpès simplex, l’interaction des récepteurs des cytokines et des cytokines, l’hépatite C et la voie de signalisation de l’adipocytocine. Nous avons examiné la littérature et identifié 40 gènes de fonction immunitaire supplémentaires (tableau S13), obtenant un total de 41 régions contenant 77 gènes liés à l’immunité. Entre ceux-ci, FST il est particulièrement élevé MUC6 région, et nous avons détecté 16 mutations de polymorphisme faux-sens simple nucléotide (SNP) montrant FST > 0,88 dans ce gène (fenêtre FST = 0,89, rapport π ln = 1,01 et XP-EHH = 5,25; Fig.3C et tableau S14). On s’attendait à ce que deux des 16 mutations faux-sens soient délétères et ces allèles délétères se produisent à une fréquence très basse (moyenne = 0,038) dans la population caprine domestique (fig. S21). Globalement, cette région contient 228 SNP presque fixes pour l’allèle dérivé de chèvres modernes (fréquence> 95%, absent chez les bézoards) qui représentent 93,8% de ces variantes dans l’ensemble du génome (pour un total de 243, rangées de données S3), qui illustre la nature singulière de ce signal de sélection.

L’analyse d’enrichissement de sélection identifie également 12 gènes impliqués dans l’interaction entre le ligand et le récepteur neuroactif (test hypergéométrique, adapté P = 3,70 × 10-5). Nous avons observé 37 autres gènes liés à d’autres fonctions du système nerveux qui étaient en cours de sélection lors de la domestication de la chèvre (tableau S15). Une région génomique située sur le chromosome 15 présente le signal de sélection combiné le plus puissant (FST = 0,90, rapport π ln = 3,20 et XP-EHH = 8,09) (Fig.3A). Essayez de grands Tajimas négatifs les valeurs, les scores élevés du rapport de vraisemblance (CLR) et le large partage des haplotypes chez les chèvres domestiques suggèrent une forte sélection positive dans ce locus (Fig. 3D et Fig. S22B). Ce locus contient 15 SNP presque fixes pour l’allèle dérivé de chèvres domestiques et abrite deux gènes codant pour des protéines, STIM1 est RRM1 (Fichier de données S3). STIM1 est un capteur de calcium du réticulum endoplasmique impliqué dans la régulation du Ca2+ et la signalisation du récepteur métabotrope du glutamate dans le système neuronal (25, 26). Il est connu pour son rôle dans la modulation de l’excitabilité des neurones de Purkinje dans le cervelet, qui jouent un rôle clé dans l’apprentissage moteur et dans l’intégration de l’information sensori-motrice (27). RRM1 code la grande sous-unité ribonucléoside-diphosphate réductase, une enzyme essentielle à la production de désoxyribonucléotides et influence la sensibilité aux défauts du tube neural induits par l’acide valproïque chez la souris (28). Par conséquent, nous émettons l’hypothèse que ce locus fortement sélectionné peut être lié à la fonction neuronale et / ou au comportement.

introgressé MUC6 joue un rôle dans la résistance aux agents pathogènes

le MUC6 La région est la seule intersection entre les analyses de sélection et d’introgression. Chez les ovins et les bovins, MUC6 est associée à la résistance du parasite gastro-intestinal (29, 30). Les résultats du séquençage du transcriptome, de la réaction de polymérisation en chaîne quantitative en temps réel (qPCR) et de l’immunohistochimie ont révélé que MUC6 elle s’exprime expressément dans la caillette caprine et le duodénum (Fig. 4A et fig. S23). Structurellement, le MUC6 possède un nombre important et variable de répétitions en tandem (VNTR) riches en résidus Pro, Thr et Ser, qui peuvent influencer l’attachement covalent des O-glycanes (31). Nous avons donc utilisé le séquençage du transcriptome PacBio SMRT pour étudier la différence entre MUC6 est MUC6B transcriptionnellement. En plus d’un certain nombre de petits indels, une suppression de 82 acides aminés contenant trois copies de VNTR après le 2789e acide aminé a été identifiée de manière unique dans MUC6 par rapport à MUC6B (fig. S24 et S25). Le nombre différent de VNTR dans ces deux haplotypes peut influencer la fonction de MUC6, qui est la clé pour générer du gel muqueux gastro-intestinal (23, 32). Examiner la différence de potentiel de résistance pathogène MUC6 est MUC6B variantes, nous avons réalisé une enquête épidémiologique dans une population caprine polymorphe. Évaluation du nombre d’œufs fécaux (FEC) pour les nématodes gastro-intestinaux chez 143 MUC6/MUC6, 111 MUC6/MUC6Bet 14 MUC6B/MUC6B moutons à 13 mois, nous avons constaté que les chèvres portent MUC6/MUC6 présentaient une FEC plus faible que ceux portant les deux autres génotypes (fichiers de données des figures 4B et S4), ce qui implique que les chèvres portant MUC6 il peut être plus résistant aux nématodes gastro-intestinaux. Pour vérifier le fond génétique, nous avons également séquencé un sous-ensemble de 117 animaux à 5 profondeurs moyennes (fichier de données S4). Une analyse d’association à l’échelle du génome pour la FEC par incorporation des covariables des composants principaux et corrélation par paires a trouvé un signal d’association fort qui contient MUC6 locus sur le chromosome 29 (de 42,77 à 51,33 Mo) (Fig.4C). Ces résultats soutiennent introgresso MUC6 chez les chèvres domestiques, il est très probablement sous sélection en raison de son avantage dans la réponse immunitaire innée de l’hôte aux pathogènes gastro-intestinaux potentiels (23).

Fig.4 L’association de MUC6 avec résistance aux nématodes gastro-intestinaux.

(UNE) Immunohistochimie pour MUC6 dans un bulbe d’aboma pylorique et duodénal de chèvre. Des micrographies 4 × et 20 × sont présentées à gauche et au centre. Les contrôles négatifs (20 ×) sont indiqués à droite. (B) L’association statistique entre MUC6 génotype et FEC pour les nématodes gastro-intestinaux. Le test de somme de rang de Wilcoxon a été utilisé pour calculer la P valeurs. (C) Diagramme de Manhattan des résultats de l’association au niveau du génome pour la FEC. Les SNP importants à l’intérieur MUC6 les locus sont surlignés en vert foncé. La ligne horizontale en pointillés indique le seuil (−Logdix(P) = 7,65).

« data-icon-position = » « data-hide-link-title = » 0 « >

Fig.4 L’association de MUC6 avec résistance aux nématodes gastro-intestinaux.

(UNE) Immunohistochimie pour MUC6 dans un bulbe d’aboma pylorique et duodénal de chèvre. Des micrographies 4 × et 20 × sont présentées à gauche et au centre. Les contrôles négatifs (20 ×) sont indiqués à droite. (B) L’association statistique entre MUC6 génotype et FEC pour les nématodes gastro-intestinaux. Le test de somme de rang de Wilcoxon a été utilisé pour calculer la P valeurs. (C) Diagramme de Manhattan des résultats de l’association au niveau du génome pour la FEC. Les SNP importants à l’intérieur MUC6 les locus sont surlignés en vert foncé. La ligne horizontale en pointillés indique le seuil (−Logdix(P) = 7,65).

L’origine et la diffusion du domestique STIM1-RRM1 est MUC6 allèles

Une étude génétique approfondie des génomes antiques à travers le Proche-Orient a révélé que deux chèvres des Balkans datant d’il y a 8100 ans portaient STIM1RRM1 (fig. fichiers de données S28 et S5). MUC6 il est apparu plus tard il y a environ 7200 ans dans le sud-ouest de l’Iran (fig. S28 et fichier de données S6), une période caractérisée par une plus grande densité de peuplements dans le croissant fertile (33). L’élevage de chèvres de densité plus élevée dans un environnement anthropique surpeuplé et sujet aux maladies aurait probablement exercé une pression plus sélective pour la résistance des agents pathogènes du bétail (34). Les premiers animaux anciens détectés ont STIM1RRM1 ou MUC6 les deux portent l’haplogroupe mitochondrial A, bien que cet haplogroupe mitochondrial ait une faible fréquence et une distribution étroite avant il y a 7500 ans (6). Il y a environ 6500 ans, la fréquence des STIM1RRM1 est MUC6 augmenté à ~ 60% au Proche-Orient (figure 5A) et s’est étendu à la Chine il y a environ 3900 ans (figure 5B), conjointement avec la consolidation et la propagation des économies d’élevage dans toute l’Eurasie (35, 36). L’expansion de ces deux loci sélectionnés était également contemporaine de la propagation écrasante de l’haplogroupe mitochondrial A. Au contraire, les fréquences des deux haplogroupes principaux du chromosome Y sont restées relativement inchangées au fil du temps (figure 5B).

Fig. 5 L’émergence et la diffusion du domestique STIM1RRM1 est MUC6 les haplotypes sont simultanés à la propagation de l’haplogroupe A. ADNmt

(UNE) Modifications temporelles de la fréquence des STIM1RRM1, MUC6et l’haplogroupe A de l’ADNmt des bézoards à la pré-domestication des chèvres domestiques modernes. Les dates sont exprimées en cal. BP. (B) Génotypes de STIM1RRM1 est MUC6et les haplogroupes de l’ADNmt et du chromosome Y. Les présences dans les états homozygotes ou hétérozygotes sont montrées respectivement en vert et vert clair. L’absence de l’allèle domestique est montrée en rose clair. La couleur gris clair symbolise les informations manquantes.

« data-icon-position = » « data-hide-link-title = » 0 « >

Fig. 5 L’émergence et la diffusion du domestique STIM1RRM1 est MUC6 les haplotypes sont simultanés à la propagation de l’haplogroupe A. ADNmt

(UNE) Modifications temporelles de la fréquence des STIM1RRM1, MUC6et l’haplogroupe A de l’ADNmt des bézoards à la pré-domestication des chèvres domestiques modernes. Les dates sont exprimées en cal. BP. (B) Génotypes de STIM1RRM1 est MUC6et les haplogroupes de l’ADNmt et du chromosome Y. Les présences dans les états homozygotes ou hétérozygotes sont montrées respectivement en vert et vert clair. L’absence de l’allèle domestique est montrée en rose clair. La couleur gris clair symbolise les informations manquantes.

DISCUSSION

La présente étude a généré des données génomiques à partir d’un nombre considérable d’animaux domestiques et sauvages Chèvre en particulier pour caractériser l’introgression adaptative et les changements génétiques lors de la domestication des chèvres. Collectivement, les analyses d’enrichissement sélectif par balayage et les deux signaux de sélection exceptionnels trouvés dans MUC6 est STIM1RRM1 chez les chèvres domestiques sont compatibles avec l’hypothèse selon laquelle les gènes liés à la résistance aux agents pathogènes (37) et le comportement a été particulièrement visé lors de la domestication des chèvres.

Numerosi studi hanno dimostrato che l’introgressione adattiva può fornire alleli benefici che consentono alle popolazioni destinatarie di adattarsi a nuovi ambienti (12, 38, 39). Nell’uomo, i loci immuno-correlati hanno acquisito alleli sostanzialmente vantaggiosi dagli esseri umani moderni immigrati, mescolandosi con umani arcaici che probabilmente erano ben adattati agli ambienti locali e ai patogeni attraverso la loro lunga esposizione ad essi (4042). Qui riportiamo numerosi segmenti genomici putativamente introgressi da una specie simile allo stambecco. Questi segmenti introgressi sono arricchiti di geni con una funzione immunitaria, suggerendo che il flusso genico storico da specie selvatiche ha svolto un ruolo importante nel modellare la diversità dei fenotipi del sistema immunitario nella capra domestica e forse ha aumentato il suo potenziale adattativo. Quando si confronta l’insieme di alleli putativamente introgressi in ciascun segmento con i nostri genomi simili a stambecchi campionati, troviamo che la maggior parte dei tassi di corrispondenza vanno da 0,5 a 0,8 (fig. S17), indicando il grado di somiglianza tra l’introgressione e il nostro simile stambecco in sequenza individui. Ulteriori selvaggio Capra le sequenze del genoma in futuro potrebbero chiarire l’origine di questi alleli, sebbene sia anche possibile che l’introgressione provenga da un ramo ormai estinto sul Capra albero.

In particolare, una serie di linee di prova confermano che il flusso genico tra il turbo del Caucaso occidentale e le capre si è verificato prima dell’inizio dell’addomesticamento. Il tur del Caucaso occidentale è distribuito intorno alla costa del Mar Nero (Fig. 2A), geograficamente vicino al centro di domesticazione delle capre. Abita in un ambiente subtropicale umido in cui l’esposizione ai patogeni e il carico di parassiti previsti possono essere considerevolmente più elevati rispetto alle regioni più aride del sud-ovest asiatico. La rigorosa identificazione dei segmenti di introgressione nelle capre domestiche mostra che il quasi corretto MUC6D fu introgredito da una specie turbolenta del Caucaso occidentale. Sebbene il primo segno di introgressione dal turbo del Caucaso occidentale sia vecchio (> 47.000 anni fa), non abbiamo trovato MUC6D in antichi e moderni bezoar campionati. Due possibilità possono spiegare questa osservazione. Innanzitutto, dato il flusso genico pervasivo tra i diversi Capra specie (8, 9, 43), non possiamo escludere la possibilità di introgressione multipla da specie simili allo stambecco prima e dopo l’addomesticamento delle capre. In secondo luogo, a causa della scarsità di bezoari antichi e moderni, può darsi che le varianti introgresse si separassero a basse frequenze in queste popolazioni. Tuttavia, i nostri risultati indicano che l’introgressione di vantaggioso MUC6D potrebbe migliorare la reattività immunitaria innata contro potenziali agenti patogeni gastrointestinali in una nicchia agro-ecologica al pascolo (44) con pressioni di malattie infettive precedentemente sconosciute.

Il nostro set di dati che include diversi campioni di capre antiche ci ha permesso di tracciare provvisoriamente l’emergenza e la diffusione delle varianti vantaggiose in questi due importanti loci domestici (Fig. 5). Le prime occorrenze di STIM1RRM1D est MUC6D il locus coincideva con l’aplogruppo A del DNA mitocondriale (mtDNA) altrimenti raro, e i tre aumentavano di frequenza approssimativamente contemporaneamente (Fig. 5A), diventando quasi fissi nelle popolazioni di capre moderne. Questi risultati suggeriscono che la diffusione globale delle variazioni centrali del processo di domesticazione è stata probabilmente mediata da un sottogruppo di capre femmine che trasportano l’aplogruppo A. mtDNA. Nonostante ciò, le popolazioni di capre moderne rimangono differenziate l’una dall’altra. The most likely explanation we can think of is that even Neolithic goat husbandry entailed some kind of breeding strategy under which immigrant matrilines containing globally advantageous alleles were deliberately crossbred with local populations, which probably carried local genetic adaptations, rather than simply replacing local populations.

Overall, our results provide evidence for adaptive introgression and the genetic basis of selected traits during the domestication of goats. We highlight that livestock domestication is a dynamic evolutionary process, with adaptive leaps driven by introgression and selection. Our results indicate that domestication may have profound impacts on neural traits and pathogen resistance, which helped manage herds to adapt to an anthropogenic environment.

MATERIALS AND METHODS

Sample collection

We collected 106 samples including 88 domestic goats (C. hircus), one bezoar (C. aegagrus), three Markhors (C. falconeri), three Siberian ibex (C. sibirica), one Alpine ibex (C. ibex), one West Caucasian tur (C. caucasica), five ancient goats spanning from the Neolithic to the Middle Ages, and four Nubian ibex hybrids (C. nubiana × C. hircus) for whole-genome sequencing. Among them, the Nubian ibex hybrids offer a potential source for determining whether the introgressed haplotypes could originate from the Nubian ibex or closely related ibex. All procedures involving sample collection were approved by the Northwest A&F University Animal Care Committee (permit number: NWAFAC1019), making all efforts to minimize invasiveness. le C. caucasica (tooth) sample was obtained from a zoo specimen donated to Ménagerie du Jardin des Plantes in 1976 of unknown provenance (wild or captive-born) but with no evidence of recent admixture (text S1). The individual died in 1982 and was subsequently sampled by the Muséum National d’Histoire Naturelle (MNHN, sample ID MNHN ZM AC 1982-1092). In addition, we also obtained whole-genome data and their geographical coordinates of the sampling sites from previous studies (5, 15). Details of the samples used in this study are given in tables S1 to S4.

DNA extraction and sequencing

For modern samples, genomic DNA was extracted from whole blood using the standard phenol-chloroform method and checked for quality and quantity on the Qubit 2.0 fluorometer (Invitrogen). Library construction for resequencing was performed with 1 to 3 μg of genomic DNA using standard library preparation protocols and insert sizes from 300 to 500 base pair (bp). All libraries were sequenced on either the Illumina HiSeq 2500 or X Ten platforms.

For the historical sample (C. caucasica), tooth pulverization, DNA extraction, and Illumina sequencing library preparation were performed in Trinity College Dublin, Ireland as described in (6). The double-stranded DNA sequencing library was constructed without uracil-DNA glycosylase treatment. The sample library was then sequenced using a HiSeq 2500 platform (single end, 1 × 101 bp). Details on the processing of five ancient samples are described in the Supplementary Materials.

Read alignment and variant calling

Filtered reads from all individuals were aligned to the latest goat reference genome by the Burrows-Wheeler Aligner (BWA v0.7.15). To obtain high-quality SNPs, we carried out SNP calling and genotyping at two separate stages. Variants were found on a population scale (without outgroup argali Ovis ammon) using a SAMtools model implemented in the analysis of next-generation sequencing data (ANGSD) (45) with “-only_proper_pairs 1 -uniqueOnly 1 -remove_bads 1 -minQ 20 -minMapQ 30 -C 50 -doMaf 1 -doMajorMinor 1 -GL 1 -setMinDepth [Total_read_depth/3] -setMaxDepth [Total_read_depth*3] -doCounts 1 -dosnpstat 1 -SNP_pval 1.” The sites that could be called in at least 90% of the samples and had a strand bias score below 90% were retained. UN P value threshold of 1 × 10−6 and minor allele frequency value > 1/2n were used as SNP discovery criteria, where n was the sample size. We also kept the sites at which the sampled individuals were homozygous for the alternative alleles. To obtain the hard-called genotypes in all Capra species and outgroups, we used a workflow adapted from GATK v3.7.0 HaplotypeCaller (46) best practice. Briefly, the genome variants, in genomic variant call format (gVCF) for each sample, were identified with the HaplotypeCaller. After all the gVCF files were merged, a raw population genotype file with the SNPs and indels was created. Optional variant hard filter values were applied for SNPs using GATK VariantFiltration as Quality by Depth (QD) < 2.0, root mean square of Mapping Quality (MQ) < 40.0, Fisher Strand (FS) > 60.0, HaplotypeScore >13.0, and MQRankSum ≤ 12.5. Last, only biallelic SNPs identified by both GATK and ANGSD that all show no more than 10% missing data were extracted using VCFtools v0.1.15 (47) with “–min-alleles 2 –max-alleles 2 –max-missing 0.9.” Then, all sample sets of filtered variant calls were used for imputation and phasing using Beagle v4.1 (48, 49) with default parameters, except for “niterations,” which was set to 10.

To retrieve the nucleotides that were accessible to variant discovery (used in the demographic estimation), hard filtering of the invariant sites was carried out with thresholds based on sequencing coverage (one-third– to threefold of corresponding mean depth) and missing data rate (no more than 10%) following the same strategies adopted for variant sites.

Population structure and phylogenetic analysis

Neighbor-joining tree was constructed for the whole-genome SNP set (table S5) using MEGA v6.0 based on pairwise genetic distances matrix calculated by PLINK v1.9. We also made use of the 5,043,096 fourfold degenerate sites to build an ML tree using RAxML v8.2.9. PCA was performed using smartpca, part of the EIGENSOFT v6.1. A Tracy-Widom test was used to determine the significance level of the eigenvectors. ADMIXTURE was used to perform an unsupervised clustering analysis. We increased the number of predefined genetic clusters from k = 2 to k = 7 and ran the analysis 20 times for each k. We further compared the individual-level haplotype similarity using fineSTRUCTURE v2.1.3. TreeMix v1.12 was also used to infer a population-level phylogeny using the ML approach.

Demographic reconstruction

To infer patterns of effective population size and population separations over historical time, we used MSMC2, using the statistically phased autosomal genotypes. We also used the SMC++, which does not rely on haplotype phase information. In addition, we calculated the likelihood of the observed site frequency spectrum conditioned on several demographic models using ∂ài (fig. S13). For the best model, nonparametric bootstrapping (100 replicates) was performed to determine the variance of each parameter (fig. S14 and table S7). To obtain reliable time estimates, we used a phylogenetic comparison (50) of goat and sheep (Ovis aries) to calibrate the mutation rate for goats. Using an average generation time of 2 years for goats (51), we estimated a mutation rate of 4.32 × 10−9 per site per generation, which is similar to that obtained in (52).

Gene flow analysis

We computed f3 statistics and D statistics to formally test hypotheses of gene flow. To identify regions introgressed into modern domestic goats from ibex-like species, we calculated D statistics using nonoverlapping 20-kb sliding windows across the genome. Windows within top 5% of the D statistics were considered as outliers and were further examined. An introgressed segment should have low sequence divergence to the putative donor species, whereas sequence divergence between the donor species and bezoar should be higher. Therefore, we also calculated the identity by state distance matrix of the outlier regions using ANGSD. Differences in sequence divergence were determined by a t test–based approach using a P value of 0.01 as cutoff. The windows which passed the above two filtering steps were merged as the significantly diverged regions.

We then applied Sprime (53) for investigating whether introgression can explain the significantly diverged regions. The first step of Sprime is to read the phased genotypes of the target and outgroup individuals. The outgroup is a population that is closely related to the ancestor(s) of the target population but that is not expected to have experienced admixture from ibex-like species that contributed to the target population. The 24 modern bezoars are used as the outgroup population, although they may themselves contain introgressed sequence from ibex-like species. This may lead to the occurrence of false negatives but will minimize the risk of false positives. We consider this an acceptable trade-off given our objective. Sprime first groups variants into three classes: those common in the outgroup, those not seen in the outgroup, and those uncommon but present in the outgroup. The threshold frequency used to distinguish common versus uncommon variants is 0.01 (53) for the analyses presented here. The common alleles (those with frequency > 0.01 in 24 modern bezoars) are excluded from further consideration because these variants in domestic goats are more likely inherited from their wild progenitor (bezoar). We assumed a score threshold of 150,000 (53) and a mutation rate of 4.32 × 10−9 per site per generation to identify introgressed alleles in domestic goat genomes. The next step is to find putatively introgressed segments (i.e., sets of alleles) for each significantly diverged region. We process one region at a time and find the set of alleles with the highest score. These comprise the putative introgressed segment. The segment boundaries were defined by the position of first and last introgressed alleles. We consider only segments with at least 100 putatively introgressed alleles that are of nonbezoar origin. To determine whether a target individual carries the introgressed haplotype, at least 50% of introgressed alleles must be found in that haplotype. We primarily focus on the most common introgressed regions (i.e., introgressed haplotype frequency > 0.1) from ibex-like species to domestic goats.

We further checked the introgression status by constructing ML trees with MEGA v6.0 and searching for domestic goats located within the ibex-like clade (fig. S16). For each region, we also reported a match if the genotype of ibex like species included the putative introgressed allele and a mismatch otherwise. The match rate was calculated as the number of matches divided by the number of compared sites (matches and mismatches) (fig. S17). Note that the match rate should not imply that the sampled ibex-like genomes represent the direct source population(s), but the degree of similarity between the introgressing and sampled ibex-like species.

Selective sweep analysis

We screened for sweeps selected during goat domestication by parsing specific 50-kb windows that showed low diversity in domestic goats and had high divergence (a high fixation index, FST, and haplotype differences, XP-EHH) between domestic goats and wild goats. After calculating all tests, the windows with P values less than 0.005 (Z test) were considered significant signals. Candidate genes under selection were defined as those overlapped by sweep regions or within 20 kb of the signals. Two additional statistics, the Tajima’s D and CLR test were applied to confirm the top signals.

Functional enrichment analyses

We characterized the most relevant functions of the protein-coding genes overlapping with the introgressed regions and selective sweeps by searching for overrepresented KEGG pathways. Goat protein sequences were used to conduct functional enrichment tests on the target genes using the KOBAS 3.0 server (54). le P value was calculated using a hypergeometric distribution. False discovery rate (FDR) correction was performed to adjust for multiple testing. Pathways with an FDR-corrected P value of <0.01 were considered statistically significantly enriched (tables S10 and S12).

Genotyping loci underlying domestication in ancient genomes

To investigate the genotypes of the STIM1RRM1 est MUC6 loci in ancient samples, we used a total of 243 SNPs (15 SNPs within STIM1RRM1 locus and 228 SNPs within MUC6 locus, respectively), which showed derived allele frequency > 0.95 in 164 modern domestic goats (defined as domestic haplotypes: STIM1RRM1D est MUC6D) and was absent in 24 modern bezoars and 4 ancient bezoars (Hovk1, Direkli1-2, Direkli6, and Direkli5) (defined as bezoar haplotypes: STIM1RRM1B est MUC6B). Because of the low coverage for some ancient genomes, we used allele reads count at SNP positions to determine the genotypes of ancient goats at STIM1RRM1 locus and MUC6 locus. For each ancient goat and SNP, we determined the numbers of reads corresponding to the ancestral and derived allele using GATK v3.7.0 UnifiedGenotyper (–min_base_quality_score 15–output_mode EMIT_ALL_SITES –standard_min_confidence_threshold_for_calling 20) (46). For each locus in each ancient goat, we calculated the proportion of the summed derived allele reads count versus total allele reads count. If the proportion was lower than 10%, the genotype of ancient goat was considered to be homozygous for ancestral allele (bezoar like). The genotype was set to heterozygous if the proportion was between 10 and 90%. If the proportion was more than 90%, it was classified as homozygous for derived allele (domestic like). We set the genotype to be missing when the ancient goats only had one single informative read (mapping to the predefined SNPs) mapped to the locus.

Expression and epidemiologic survey related to MUC6 locus

qPCR analysis. Total RNA was extracted using goat tissues from gastrointestinal tract (including abomasum, duodenum, jejunum, ileum, and caecum) using Eastep Super Total RNA Extraction Kit (Promega, Shanghai, China). Complementary DNA was generated via PrimerScript RT Reagent Kit with genomic DNA Eraser (Perfect Real Time) (TaKaRa, Beijing, China). qPCR analysis was performed with FastStart Universal SYBR Green Master (Roche, Shanghai, China) on a Bio-Rad instrument. MUC6 primers (forward primer: 5′-CAGCCAGGACAAAATCATGA-3′ and reverse primer: 5′-CTCTGGTCTGGCCTCTGAAC-3′) were designed using National Center for Biotechnology Information Primer-BLAST (http://ncbi.nlm.nih.gov/tools/primer-blast/). GAPDH (forward primer: 5′-ACACCCTCAAGATTGTCAGC-3′ and reverse primer: 5′-ATAAGTCCCTCCACGATGC-3′) was used as internal reference. Gene expression results were calculated using the delta-delta cycle threshold (2−ΔΔCT) method.

Immunohistochemistry analysis. For histological analysis, dissected goat tissues (abomasum pyloric and duodenal bulb) were fixed in 4% paraformaldehyde and embedded in paraffin for sectioning (5 μm thick). The tissues sections were deparaffinized, and rehydration was followed by antigen retrieval using a citrate-buffered solution in a microwave at 100°C for 15 min. After cooling down to room temperature, quenching of endogenous peroxidase, and protein block, the sections were treated with a primary antibody (cat. no. D161001, Sangon Biotech, Shanghai, China) overnight at 4°C. Negative control sections were obtained by incubating with the primary antibody diluted buffer. Subsequently, the secondary antibody (cat. no. D110058, Sangon Biotech, Shanghai, China) was used to detect primary antibody. Specific protein immunoreactivity was visualized with the substrate chromogen 3,3′- diaminobenzidine. Last, the slides were rinsed in water, counterstained with hematoxylin, and mounted with coverslips.

Full-length transcript sequencing. To determine the sequence differences between the domestic (MUC6D) and bezoar (MUC6B) MUC6 haplotypes, total RNA was isolated from the abomasum pyloric with high expression of MUC6 from two heterozygous goats and sequenced by PacBio Sequel. Pacbio SMRTbell libraries were sequenced on two separate SMRT cells (Annoroad Gene Technology Co. Ltd., Beijing, China). A total of 23,338,976 subreads were generated with a mean accuracy of 80% and an average length of 1755 nt. High-quality circular consensus sequences (CCSs) were obtained using the Iso-Seq 3 application in the PacBio SMRT Analysis v6.0.0 (https://github.com/PacificBiosciences/IsoSeq), with parameters “–noPolish –minPasses 1.” Last, we got a total of 960,271 CCSs. These CCSs were aligned to the latest goat reference genome using Minimap2 (55), and we picked out 251 CCSs that mapped to the MUC6 mRNA sequence. We used the CCS that belonged to the bezoar haplotype to manually assemble the mRNA sequence of the MUC6B. Then, we realigned the MUC6B mRNA sequence to MUC6D mRNA sequence (XM_018042766.1) by MEGA v6.0 and carefully checked the indels. The abundance of tandem repeats in MUC6 (XP_017898255.1) was examined by using the rapid automatic detection and alignment of repeat finding program (https://www.ebi.ac.uk/Tools/pfa/radar/) (56). Three major types of repeat units with distinct sequence features were observed. We detected a 246-bp insertion in the MUC6B mRNA sequence located at the 32nd exon, containing three copies of type III units between the 2789th amino acid and the 2871th amino acid (fig. S24). We further validated this insertion by aligning short reads from West Caucasian tur to the MUC6B mRNA sequences (fig. S25).

Test population. To further test the genotype-phenotype associations, we first examined the genotypes at the MUC6 locus using blood samples from breeding rams from a company located in Inner Mongolia, China. This company runs more than 9000 cashmere goats, mainly for superfine cashmere (fiber diameter < 14 μm). The animals are dewormed routinely with oxfendazole, ivermectin, avermectin, and levamisole three times annually. Of the 20 tested breeding rams, we identified 4 animals with MUC6D/MUC6B genotype. We sampled their offspring partly on the basis of the breeding records. Blood samples of 495 animals were collected in January 2019.

Genotypic and phenotypic analysis. DNA was extracted from the blood samples. le MUC6 locus was genotyped by PCR (forward primer: 5′-CAGCACTATCTCCCATACATC-3′ and reverse primer: 5′-GTGGAGCTGAGCTGACACTT-3′) and Sanger sequencing. The frequency of MUC6B was 17.3% (n = 338 MUC6D/MUC6D, n = 143 MUC6D/MUC6B, est n = 14 MUC6B/MUC6B).

After genotyping, we collected fresh fecal samples from a subset of 268 animals from the rectum in April 2019, a time when gastrointestinal nematodes numbers were anticipated to be elevated. The gastrointestinal parasitic infestations were examined using the McMaster’s technique as described in (57). Briefly, 2 g of fecal samples were transferred into a container, and 60 ml of saturated sodium chloride was added. The suspension was thoroughly stirred with a glass stick and poured through a strainer (80-mesh screen) into the new container. The container with the suspension was closed tightly and carefully inverted several times. Then, the suspension was taken up from the container to fill the glass plate with two McMaster counting chambers. The size of the counting chamber was 10 mm by 10 mm by 1.5 mm. The nematode eggs were then counted under microscopic observation at ×100 magnification. Observed nematodes in the feces included the common abomasal and intestinal goat nematodes: Hemonchus contortus est Nematodirus sp. (fig. S26).

The number of nematodes eggs per gram (EPG) was calculated according to the following formula: EPG = (n1 + n2)/2 ÷ 0.15 × 60 ÷ 2, where (n1 + n2)/2 is the average number of eggs per chamber, 0.15 is the effective volume (milliliters) for counting chamber, 60 is the total volume (milliliters) of suspension, 2 is the weight (grams) of feces examined. EPG was measured on three replicates of each fecal sample, and the average of the three replicates was used for analysis.

Data analysis. The average fecal nematode egg counts were not normally distributed; therefore, we used Wilcoxon rank sum test to test the null hypothesis that there was no association between genotype and phenotype.

Genome-wide association study. To have more even representation of sampling of the three genotypes, we used 10 MUC6B/MUC6B samples and randomly chose 46 MUC6B/MUC6D and 61 MUC6D/MUC6D for whole-genome sequencing (data file S4). The 117 individuals were sequenced at approximately fivefold genome coverage using BGISEQ-500 platform. SNP calling and estimation of genotype likelihoods were performed in ANGSD using the following settings: “-baq 1 -C 50 -only_proper_pairs 1 -uniqueOnly 1 -remove_bads 1 -minQ 20 -minMapQ 25 -doMaf 1 -minInd 106 -skipTriallelic 1 -doMajorMinor 1 -GL 1 -setMinDepth 250 -setMaxDepth 1,000 -doCounts 1 -doGlf 2 -SNP_pval 1e-6.” To correct for population stratification and cryptic relatedness, we adopted genotype likelihood–based approaches implemented in PCAngsd v0.973 (58) to perform PCA and estimate pairwise relatedness by supplying the following options: “-kinship -inbreed 2 -minMaf 0.05.” Beagle v3.3.2 (48) was used to convert the genotype likelihoods to actual genotypes. Only sites for which the allelic R2 (larger values of allelic R2 indicate more accurate genotype imputation) was >0.95 were retained. The imputed genotype dataset was converted to tped/tfam format using PLINK v1.9. A total of 12.47 million high-quality SNPs (minor allele frequency > 0.05, using whole-genome sequencing data alone) were used to perform genome-wide association study in EMMAX (59) software by incorporation of the first three principal component covariates (fig. S27, A and B) and kinship matrix. The phenotype was rank-transformed to normality (Shapiro-Wilk test, P = 0.02228) to counteract departures from normality (fig. S27C).

The genome-wide critical value was determined by permutation: The phenotype data were permuted 200 times; for each permutated phenotype, a genome-wide association analysis was conducted and the genome-wide lowest P value was recorded. We then took the 5% lowest tail from the 200 recoded minimal P values as the threshold for genome-wide significance (P = 2.257−8). The Manhattan and quantile-quantile plots (fig. S27D) were generated using the R package CMplot (https://github.com/YinLiLin/R-CMplot).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is ne pas for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

  1. F. Pereira, A. Amorim, in Origin and Spread of Goat Pastoralism (John Wiley & Sons, 2010).

  2. D. M. Shackleton, Wild Sheep and Goats and their Relatives: Status Survey and Conservation Action Plan for Caprinae (International Union for Conservation of Nature and Natural Resources, 1997).

  3. J. A. Lenstra, Evolutionary and demographic history of sheep and goats suggested by nuclear, mtDNA and Y chromosome markers, in The Role of Biotechnology for the Characterization and Conservation of Crop, Forestry, Animal and Fishery Geneic Resources (International Workshop, Turin, Italy, 5 to 7 March 2005), pp. 1–4.

Acknowledgments: We thank members of the NextGen project for sharing data. We thank L. Mohammed for providing Nubian ibex hybrid blood samples. We thank High-Performance Computing (HPC) of NWAFU for providing computing resources. Funding: This project was supported by grants from the National Natural Science Foundation of China (31822052 and 31572381 to Y.J. and 31572369 to Y.Ch.), the Talents Team Construction Fund of Northwestern Polytechnical University (NWPU) and Strategic Priority Research Program of CAS (XDB13000000) to W.W., a DFF Sapere Aude grant to R.H., and the Tang Scholar at Northwest A&F University (NWAFU) to Xiaolong Wang. The Chinese Government contribution to CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources in Beijing (2018-GJHZ-01) is appreciated. Author contributions: Y.J. and Y.Ch. conceived the project and designed the research. Z.Z., M.L., and Xihong Wang performed most of the analysis with contributions from Z.Y., Y.L., Y.Ca., Q.C., J.Li., K.W., X.P., Y.W., S.H., M.G., T.Z., Y.Z., Y.G., Y.X., N.X., and Y.Y. Xiaolong Wang, W.Z., J.H., L.Ch., A.E., and M.O. prepared the modern DNA samples. J.Le. provided the historical sample of West Caucasian tur. D.C. prepared the ancient samples. Xihong Wang, Z.Z., Y.J., and M.L. drafted the manuscripts with input from all authors. Y.J., W.W., R.H., G.Z., K.G.D., D.G.B., L.Co., and T.S.S. revised the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Individual genome sequence data are available at the Sequence Read Archive (https://ncbi.nlm.nih.gov/sra) under accession codes PRJNA387635 and PRJNA361447.