Signatures génomiques convergentes de la domestication chez les moutons et les chèvres | Nature Communications

Échantillonnage

Des moutons (O. aries) et des chèvres (C. hircus) domestiques ont été échantillonnés en Iran (groupes IROA et IRCH, respectivement) et au Maroc (groupes MOOA et MOCH, respectivement) pour un total de 20 animaux par groupe (figure supplémentaire 6). Ces échantillons ont été collectés entre janvier 2008 et mars 2012 dans la partie nord du Maroc et entre août 2011 et juillet 2012 dans le nord-ouest de l’Iran, dans le cadre du projet européen Nextgen (Grant Agreement no. 244356), conformément aux règles éthiques de la Directive 86/609/CEE de l’Union européenne. Des clips d’oreille ont été prélevés dans la partie distale de l’oreille d’animaux choisis au hasard, et immédiatement stockés dans de l’éthanol à 96% pendant une journée avant d’être transférés dans des billes de gel de silice jusqu’à l’extraction de l’ADN.

Les espèces sauvages mouflon d’Asie (O. orientalis) et bouquetin de Bézoard (C. aegagrus) ont été échantillonnées dans le nord-ouest de l’Iran au sein du berceau de domestication21,22. Treize mouflons d’Asie et 18 bouquetins du Bézoard (respectivement, groupes IROO et IRCA, Fig. 6 supplémentaire) ont été prélevés sur des animaux captifs ou récemment chassés, et sur des échantillons congelés disponibles au Département iranien de l’environnement. Cette approche d’échantillonnage basée sur les individus est conçue pour minimiser les biais potentiels en évitant la surreprésentation des effets locaux (par ex, consanguinité locale).

Données supplémentaires

En outre, un panel mondial de races a été assemblé pour les moutons et les chèvres (wpOA et wpCH, respectivement). wpOA comprenait 20 échantillons de reséquençage du génome entier (WGS) à couverture 12x représentant 20 races mondiales différentes fournies par le Consortium international du génome ovin. wpCH comprenait 14 échantillons WGS séquencés à couverture 12x représentant 9 individus européens, à savoir , 2 échantillons de moutons alpins français et 2 échantillons de moutons de race saanen français séquencés par l’INRA, 5 échantillons de moutons de race saanen italiens fournis par le Parco Tecnologico Padano, et 5 individus australiens, à savoir, 2 échantillons de Boer, 2 de Rangeland, et 1 de Cashmere fournis par le CSIRO (Données supplémentaires 5).

Production de données WGS

L’ADN génomique a été extrait avec succès de tous les échantillons de tissus à l’aide du kit Macherey Nagel NucleoSpin 96 Tissue, en adaptant le protocole du fabricant. Les prélèvements de tissus ont été effectués dans des blocs MN à puits carrés afin d’obtenir 25 mg de fragments par échantillon. Trois blocs carrés MN-96 et demi ont été préparés, et l’extraction a été réalisée à l’aide d’un manipulateur de liquide Tecan Freedom Evo en suivant le protocole du fabricant. Une étape de pré-lyse a été effectuée pour homogénéiser les échantillons avec 180 µl de tampon T1 et 25 µl de protéinase K pendant une nuit à 56 °C. Pour ajuster les conditions de liaison, 200 µl de tampon BQ1 ont été ajoutés et la plaque d’échantillons a été incubée 1 h à 70 °C ; 200 µl d’éthanol à 100 % ont ensuite été ajoutés. Les lysats ont été transférés sur la plaque de liaison Nucleospin Tissue et un vide (-0,2 bar, 5 min) a été appliqué pour éliminer le flux. Trois étapes de lavage ont été effectuées avec les tampons BW et B5, respectivement, et un vide a de nouveau été appliqué pour éliminer le flux. Avant l’élution de l’ADN génomique, une membrane de silice de la plaque de fixation des tissus Nucleospin a été séchée sous vide avec au moins -0,6 bar pendant 10 minutes. L’étape d’élution a été réalisée avec 100 µl de tampon BE préchauffé (70 °C) et une étape de centrifugation à 3700 rcf pendant 5 min dans des puits 96-PCR. L’ADN génomique a été stocké à 4 °C pour éviter la congélation-décongélation et testé pour la concentration (en ng/μl) en utilisant la méthode Picogreen et en utilisant un Nanodrop.

Des génomes entiers ont été reséquencés à partir de 500 ng d’ADN génomique qui ont été cisaillés à une gamme de 150-700 pb en utilisant l’instrument Covaris® E210 pour chaque échantillon et utilisés pour la préparation de bibliothèque Illumina® par un protocole semi-automatisé. La réparation des extrémités, la queue en A et la ligature des adaptateurs compatibles avec Illumina® (BioScientific) ont été effectuées à l’aide du système de préparation de bibliothèque SPRIWorks et de l’instrument SPRI TE (Beckmann Coulter) selon le protocole du fabricant. Une sélection de taille de 300-600 pb a été appliquée pour récupérer la plupart des fragments. Les fragments d’ADN ont été amplifiés par 12 cycles de PCR à l’aide de la trousse Platinum Pfx Taq Polymerase (Life® Technologies) et des amorces spécifiques à l’adaptateur Illumina®. Les librairies ont été purifiées avec des billes AMPure XP 0,8x (Beckmann Coulter), et analysées avec le bioanalyseur Agilent 2100 (Agilent® Technologies) et la quantification qPCR. Les librairies ont été séquencées en utilisant la chimie des lectures de 100 bases de long en cellule de flux paired-end sur le Illumina® HiSeq2000.

Les lectures paired-end d’Illumina pour Ovis ont été mappées au génome de référence du mouton (OAR v3.1, assemblage GenBank GCA_000298735.146), et pour Capra au génome de référence de la chèvre (CHIR v1.0, assemblage GenBank GCA_000317765.147), en utilisant BWA-MEM48. Le fichier BAM produit pour chaque individu a été trié en utilisant Picard SortSam et amélioré en utilisant séquentiellement Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator et GATK IndelRealigner49, et SAMtools calmd50.

La découverte de variants a été effectuée en utilisant trois algorithmes différents : Samtools mpileup50, GATK UnifiedGenotyper51, et Freebayes52. Les sites variants ont été identifiés indépendamment pour chacun des six groupes, en utilisant les modes multi-échantillons des algorithmes d’appel : (i) 162 échantillons de MOOA ; (ii) 20 échantillons d’IROA ; (iii) 14 échantillons d’IROO ; (iv) 162 échantillons de MOCH ; (v) 20 échantillons d’IRCH ; (vi) 19 échantillons d’IRCA. Pour certains groupes, les WGS de plus d’individus étaient disponibles dans le cadre du projet NextGen (voir ci-dessus). Les échantillons utilisés dans la présente étude ont été sélectionnés pour obtenir des groupes équilibrés de 20 individus dans la mesure du possible. Pour les groupes IRCA et IROO, des échantillons supplémentaires sont devenus disponibles à un stade ultérieur et ont été ajoutés pour les analyses en aval. Les animaux présentant une faible qualité d’alignement et d’appel ont été supprimés pour obtenir l’ensemble de données final (Données supplémentaires 5).

Dans chaque groupe, il y a eu deux tours successifs de filtrage de la qualité des sites variants. L’étape de filtrage 1 a fusionné les appels des trois algorithmes, tout en filtrant les appels de plus faible confiance. Un site variant est accepté s’il a été appelé par au moins deux algorithmes d’appel différents avec une qualité de variant phred >30. Un allèle alternatif sur un site est accepté s’il a été appelé par l’un des algorithmes d’appel et que le nombre de génotypes était >0. L’étape de filtrage 2 a utilisé le recalibrage du score de qualité des variants par GATK. Tout d’abord, nous avons généré un ensemble d’entraînement des sites de variante les plus fiables au sein du groupe où (i) le site est appelé par les trois algorithmes d’appel de variante avec une qualité de variante à échelle réduite >100, (ii) le site est bialélique, (iii) le nombre d’allèles mineurs est d’au moins 3 tout en ne comptant que les échantillons avec une qualité de génotype à échelle réduite >30. L’ensemble d’entraînement a été utilisé pour construire un modèle gaussien à l’aide de l’outil GATK VariantRecalibrator en utilisant les annotations de variants suivantes de UnifiedGenotyper : QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Un modèle gaussien a été appliqué à l’ensemble des données, générant un VQSLOD (log odds ratio d’être un vrai variant). Les sites ont été filtrés si VQSLOD < valeur seuil. La valeur seuil a été fixée pour chaque groupe de la manière suivante : VQSLOD minimum = {la valeur médiane de VQSLOD pour les variants de l’ensemble de formation}-3 × {la déviation absolue médiane de VQSLOD des variants de l’ensemble de formation}. Le rapport SNP de transition/transversion a suggéré que le critère de coupure choisi donnait le meilleur équilibre entre sélectivité et sensibilité.

Des ensembles d’appels SNP pour six groupes d’animaux Ovis et Capra ont été générés (c’est-à-dire les domestiques iraniens et marocains, et les sauvages pour chaque genre). Comme les analyses effectuées dans cette étude nécessitaient des comparaisons entre groupes, nous avons créé des ensembles d’appels de génotypes à un ensemble cohérent de sites SNP pour tous les animaux de n’importe quel groupe. Pour chaque genre, nous avons fusionné les sites d’appel de variants de ses trois groupes, et n’avons retenu que les positions bialéliques sans données manquantes. Les génotypes ont été réappelés à chaque site SNP bialélique pour tous les individus d’intérêt par GATK UnifiedGenotyper en utilisant l’option GENOTYPE_GIVEN_ALLELES. À ce stade, la liste des individus a été élargie pour inclure les animaux appartenant aux panels mondiaux des races ovine et caprine (wpOA et wpCH) et des échantillons sauvages supplémentaires devenus disponibles à ce stade (4 O. orientalis et 4 C. aegagrus). Les génotypes ont été améliorés et mis en phase au sein des groupes par Beagle 453, puis filtrés lorsque la probabilité du génotype était inférieure à 0,95. Enfin, nous avons filtré les sites qui étaient monomorphes dans les différents sous-ensembles d’individus utilisés dans cette étude (voir ci-dessous).

Afin de comparer les signaux de sélection détectés entre Ovis et Capra, nous avons réalisé un alignement croisé entre les deux génomes de référence. Tout d’abord, nous avons utilisé le pipeline d’alignement par paire de la base de code Ensembl release 6954 pour aligner les génomes de référence du mouton (OARv3.1) et de la chèvre (CHIR1.0). Ce pipeline utilise LastZ55 pour aligner au niveau de l’ADN, suivi d’un post-traitement dans lequel les blocs alignés sont enchaînés en fonction de leur emplacement dans les deux génomes. Le pipeline d’alignement par paires LastZ est exécuté de manière routinière par Ensembl pour toutes les espèces prises en charge, mais la chèvre n’est pas encore incluse dans Ensembl. Pour éviter tout biais en faveur de l’une ou l’autre espèce, nous avons produit deux alignements interspécifiques différents. L’un a utilisé le mouton comme génome de référence et la chèvre comme non-référence, tandis que l’autre a utilisé la chèvre comme génome de référence et le mouton comme non-référence. La différence réside dans le fait que les régions génomiques de l’espèce de référence sont forcées de s’associer de façon unique à des loci uniques de l’espèce non de référence, alors que les régions génomiques non de référence sont autorisées à s’associer à des emplacements multiples de l’espèce de référence. Nous avons obtenu pour les segments de chromosomes d’un génome de référence les coordonnées sur le génome de non-référence. Enfin, pour les SNP découverts dans un genre, nous avons utilisé l’alignement du génome entier avec le génome de référence de l’autre genre pour identifier les positions correspondantes (tableau supplémentaire 6).

Structure génétique

Afin de décrire la diversité génétique au sein des groupes, nous avons utilisé VCFtools56 pour calculer les statistiques sommaires de variation génétique sur les 73 individus pour Ovis (soit , 13 IROO, 20 IROA, 20 MOOA, et 20 wpOA) et 72 individus pour Capra (c’est-à-dire 18 IRCA, 20 IRCH, 20 MOCH, et 14 wpCH). Les statistiques mesurées étaient le nombre total de variants polymorphes (S) pour l’ensemble des individus de chaque genre et au sein de chaque groupe, la diversité nucléotidique moyenne (π) au sein de chaque groupe et le coefficient de consanguinité (F) pour chaque individu. Au sein de chaque genre, les différences entre le groupe sauvage et chaque groupe domestique ont été testées à l’aide d’un test t unilatéral pour les valeurs individuelles de consanguinité et de charge génétique, et d’un test de Mann-Whitney bilatéral pour la diversité nucléotidique par site.

La divergence globale entre les quatre groupes de chaque genre (c’est-à-dire, sauvage, domestique iranien et marocain, et panel mondial) a été estimée à l’aide de tous les SNP bialéliques et du Fst moyen pondéré par paire suivant Weir et Cockerham57 tel qu’implémenté dans VCFtools56. La structure génétique entre les groupes a été évaluée à l’aide de la méthode de regroupement sNMF26, après avoir élagué l’ensemble des données pour éliminer les SNP dont le déséquilibre de liaison (r²) était supérieur à 0,2 en utilisant VCFtools. Le déséquilibre de liaison (r²) a été calculé entre des paires de SNP dans des fenêtres glissantes de 50 SNP, un SNP par paire étant retiré de manière aléatoire lorsque r² était supérieur à 0,2. Pour chaque analyse sNMF, cinq passages du même nombre de clusters (K) ont été effectués avec des valeurs de K de 1 à 10. Nous avons utilisé le critère d’entropie croisée pour identifier la solution de clustering la plus probable, cependant, des partitions alternatives pour différents nombres de K ont également été explorées pour évaluer comment les individus étaient répartis entre les clusters.

Pour démêler l’ascendance partagée et le mélange, nous avons exécuté TreeMix27 pour estimer conjointement les divisions de population et les événements de mélange ultérieurs en utilisant l’ensemble de données élagué utilisé pour sNMF. Nous avons exécuté TreeMix avec l’option -global pour affiner nos inférences de vraisemblance maximale. Nous avons enraciné l’arbre TreeMix avec la division entre les individus sauvages et domestiques. La taille du bloc pour le jackknifing était de -k 500 SNP, ce qui correspond approximativement à 150 kb, dépassant les blocs moyens de DL trouvés chez les moutons et les chèvres. Nous avons généré un arbre de Maximum de Vraisemblance sans migration, puis nous avons ajouté des événements de migration et nous avons examiné le changement incrémental de la variance expliquée par le modèle et les valeurs résiduelles entre les individus. L’objectif était de détecter toute valeur résiduelle élevée potentielle ou tout bord de migration entre les individus sauvages et domestiques. Afin d’explorer davantage la pertinence statistique des vecteurs d’admixture possibles identifiés par TreeMix (Tableau supplémentaire 3), nous avons calculé le test des trois populations f328 comme un test formel d’introgression génétique, en utilisant le programme qp3Pop de la suite ADMIXTOOLS58 pour chaque combinaison de groupes. Pour Capra, le groupe wpCH a été divisé entre les races australiennes, les races françaises et les races italiennes. Les résultats sont rapportés dans les données supplémentaires 2.

Inférence démographique

Pour chaque genre, nous avons effectué des analyses d’inférence démographique ancestrale en utilisant le modèle MSMC implémenté dans le logiciel MSMC225. MSMC est basé sur la coalescence séquentielle Markovienne par paire59 ; cependant, il utilise en entrée les haplotypes des données de séquences génomiques phasées. Pour chaque analyse, nous avons utilisé deux individus d’un groupe, donc 4 haplotypes. Chaque analyse a été répétée pour un autre ensemble aléatoire de deux individus, c’est-à-dire une réplique de l’analyse par groupe. Les fichiers d’entrée et de sortie ont été générés et analysés avec les scripts python fournis avec le logiciel MSMC et trouvés à https://github.com/stschiff/msmc-tools. Les paramètres d’analyse ont été conservés par défaut, à l’exception du taux de mutation qui a été fixé à 2,5×10-8 et de la durée de la génération qui a été fixée à 2 ans. Afin d’estimer l’incertitude sur les estimations de temps, nous avons fait varier ces paramètres (taux de mutation de 2,5×10-8 et 1,0×10-8 en combinaison avec une longueur de génération de 2 et 4 ans) et fourni une estimation approximative de la période de domestication (voir la figure supplémentaire 2).

Charge génétique

La charge génétique a été estimée de deux façons. Tout d’abord, en calculant la charge génétique pour chaque individu comme la somme des effets de fitness délétères sur toutes les positions génomiques codant pour des protéines en suivant la méthode de Librado et al60. En bref, en tant qu’indicateur de la contrainte évolutive, nous avons utilisé les scores PhyloP de l’alignement 46 voies des mammifères (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). À partir de cet alignement, nous avons identifié les sites codant pour des protéines évoluant sous contraintes fonctionnelles (score phyloP ≥1,5). Pour chaque génome d’Ovis ou de Capra, nous avons ensuite recherché si ces sites étaient mutés. Si c’est le cas, nous avons additionné les scores phyloP sur tous les sites mutés, de sorte que les mutations dans les sites fortement contraints contribuent proportionnellement plus à l’estimation de la charge totale. Nous avons ainsi obtenu une estimation de la charge pour chaque génome de mouton/chèvre. Enfin, pour obtenir une charge moyenne par site, nous avons divisé par le nombre total de positions analysées. Il convient de noter que nous avons conditionné aux sites homozygotes pour éviter de modéliser le coefficient de dominance des mutations sur les sites hétérozygotes (par exemple, récessif, intermédiaire, dominant). Deuxièmement, nous avons comparé gène par gène la charge délétère génétique dans les groupes d’Ovis sauvages et domestiqués en effectuant un test de Wilcoxon, l’hypothèse alternative étant que les animaux domestiques ont une charge plus importante que les parents sauvages. Les valeurs p ont été corrigées pour les tests multiples61 et nous avons appliqué un seuil de valeurs p ajustées < 0,05. Nous avons effectué une analyse d’enrichissement de l’ontologie génétique sur l’ensemble des gènes montrant une augmentation significative de la charge génétique en utilisant WebGestalt62,63. Comme les génomes de référence sont mal annotés pour les gènes, nous nous sommes appuyés sur les orthologues à copie unique entre notre espèce et l’homme et la souris. Les gènes du chromosome X ont été exclus de l’ensemble de référence. Nous n’avons pas effectué cette analyse sur Capra en raison de la consanguinité plus élevée observée dans les échantillons sauvages.

Détection des signatures de sélection

Pour détecter les signatures de sélection liées à la domestication, nous avons utilisé tous les SNP bialéliques présentant une fréquence d’allèle mineur supérieure à 0,10 dans au moins un des trois groupes testés (c’est-à-dire les groupes domestiques iranien et marocain, et le groupe sauvage pour chaque genre). Comme nous nous attendions à ce que les signatures de sélection liées au processus de domestication soient présentes chez tous les animaux domestiques, nous avons adopté la stratégie générale suivante : nous avons testé avec hapFLK29 (voir la note supplémentaire 5 et les figures supplémentaires 9, 10 et 11) pour chaque genre le groupe sauvage par rapport à chacun des groupes domestiques gérés traditionnellement (c’est-à-dire l’Iran et le Maroc) et nous nous sommes concentrés sur les régions communes putativement sous sélection qui ont été détectées dans les deux cas. La taille des échantillons des groupes (n = 13-20) était compatible avec les exigences de la méthode29. Nous avons vérifié visuellement si les signatures cohérentes de sélection trouvées avec hapFLK étaient également présentes dans l’ensemble du panel mondial correspondant de chaque genre, mais nous n’avons pas inclus ces groupes dans le test statistique en raison de leur composition multi-races. Enfin, nous avons recherché des signaux de sélection partagés entre Ovis et Capra en utilisant une approche FDR stratifiée. La stratégie est décrite dans la figure supplémentaire 4.

Nous avons effectué des tests hapFLK pour contraster le groupe sauvage à chacun des groupes iraniens et marocains dans chaque genre. La matrice de parenté a été calculée à partir des distances génétiques de Reynold64 entre les paires de groupes, en utilisant un sous-ensemble aléatoire de 1 % des variants. L’arbre de population déduit a été construit en utilisant l’algorithme de jonction des voisins. Pour chaque SNP, nous avons effectué le test hapFLK qui incorpore l’information haplotypique pour augmenter la puissance de détection des balayages sélectifs. Pour chaque SNP testé, la statistique hapFLK a calculé la déviation des fréquences haplotypiques par rapport au modèle neutre estimé par la matrice de parenté65. Pour exploiter les informations sur le déséquilibre de liaison, hapFLK utilise le modèle multipoint de Scheet et Stephens66 pour les génotypes multilocus qui peut être ajusté aux données non phasées. L’une des principales applications de ce modèle est d’effectuer une estimation de la phase (logiciel fastPHASE66). Dans notre analyse, le modèle a été entraîné sur des données non phasées, et notre analyse tient donc compte de l’incertitude de la phase. La méthode a été utilisée pour regrouper les haplotypes locaux le long des chromosomes dans un nombre spécifié de clusters K fixé à 25, en utilisant un modèle de Markov caché.

Pour identifier les régions communes putativement sous sélection dans les deux groupes domestiques gérés traditionnellement pour chaque genre, nous avons combiné les deux analyses hapFLK précédentes. Pour chaque analyse, les scores hapFLK ont été ajustés à une distribution χ2 pour obtenir les valeurs p (script disponible à https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Les résultats des deux contrastes entre le groupe sauvage et chacun des groupes domestiques ont été combinés à l’aide de la méthode de Stouffer67 pour obtenir des valeurs p uniques pour la comparaison entre les animaux sauvages et les animaux domestiques. Enfin, le cadre FDR68 a été appliqué à l’ensemble des SNP pour convertir les p-valeurs combinées en q-valeurs. Les SNP présentant des valeurs q < 10-2 ont été retenus et regroupés en régions génomiques lorsqu’ils étaient distants de moins de 50 kb les uns des autres.

Pour étudier si le signal de sélection était partagé entre Ovis et Capra, nous avons d’abord utilisé l’alignement croisé des deux génomes de référence pour identifier les segments homologues. Nous avons ensuite appliqué un cadre FDR stratifié69. Cette approche est basée sur le fait qu’il y a une stratification inhérente dans les tests étant donné l’information préalable dans les données génétiques69, parce que la distribution sous-jacente des vraies hypothèses alternatives pourrait être différente selon la dynamique différente de diverses régions génomiques, conduisant à des distributions différentes des p-values. Cela nécessite d’obtenir des valeurs p ajustées par le FDR (c’est-à-dire des valeurs q) séparément pour les différentes strates. Nous avons recherché les convergences dans chaque genre en séparant les régions homologues à celles détectées dans l’autre genre (appelées strate partagée) et le reste du génome (appelé strate générale). Nous avons extrait les valeurs p séparément pour chacune des deux strates définies et avons ensuite calculé les valeurs q par le biais du cadre FDR. Ces valeurs q stratifiées ont été les quantités finales considérées pour la signification statistique (<10-2) afin de détecter les SNP sous sélection et de les fusionner dans les régions génomiques correspondantes.

Pour tester les signatures convergentes de sélection différenciant les animaux sauvages des animaux domestiques dans les deux genres, nous avons examiné la relation entre le seuil de signification appliqué aux valeurs q (que nous avons fait varier de 0,2 à 0,002) dans un genre et la probabilité estimée qu’un SNP soit sélectionné dans la strate partagée de l’autre genre en utilisant l’approche de Storey et al70. de Storey et al.70 . Une augmentation de la probabilité inférée avec une diminution du seuil appliqué à la valeur q (augmentation de la rigueur) indique que plus la région est significative dans un genre, plus nous sommes susceptibles de trouver des SNP significatifs dans l’autre genre.

Nous avons filtré les signaux de sélection qui n’étaient pas cohérents entre les trois groupes domestiques. Pour chaque région détectée, nous avons utilisé les haplotypes en phase des individus qui ont été regroupés à l’aide d’arbres Neighbor-Joining basés sur le pourcentage d’identité entre les séquences. Seules les régions présentant des signaux cohérents ont été conservées (figure supplémentaire 5).

Afin de déduire si les signaux de sélection détectés avec hapFLK indiquaient un relâchement de la sélection ou une sélection positive chez les domestiques, nous avons estimé la différence de diversité nucléotidique (π) sur chaque région putative sous sélection entre les groupes sauvages et domestiques. Nous avons exprimé cette différence sous la forme de l’indice Δπ, qui a été calculé pour chaque région génomique comme la différence entre π calculé pour le groupe sauvage et la moyenne de π pour les groupes domestiques iranien et marocain, moins la différence de π entre ces deux groupes calculée sur l’ensemble du génome :

$$\Delta \pi = \left( {\pi _{{\rm wilds}} – \pi _{{\mathrm{\i} nom de l’opérateur{{\i} iranien}}}}}} \right)_{{\mathrm{\operatorname{{genomic-region}}}} – \left( {\pi _{{\rmathrm wilds}}) – \pi _{{{mathrm{\operatorname{iran-morocco}}}}} \right)_{{\mathrm{\operatorname{whole-genome}}}}$$

Une valeur négative indiquerait que la diversité nucléotidique est plus faible dans le groupe sauvage par rapport à la moyenne des deux groupes domestiques, et serait considérée comme montrant un relâchement de la sélection dans ces derniers groupes, une sélection diversifiante chez les domestiques ou une sélection positive dans les sauvages. A l’inverse, une valeur positive indiquerait une sélection directionnelle positive ou stabilisatrice qui s’est produite dans les groupes domestiques. Nous avons également utilisé le regroupement d’haplotypes pour vérifier manuellement dans chaque région si le balayage sélectif détecté confirmait les indications données par l’indice Δπ.

Nous avons effectué des interprétations fonctionnelles comme suit. Pour chaque région sélectionnée, nous avons considéré la région plus 50 kb de chaque côté pour identifier les rôles fonctionnels et 5 kb en amont et en aval des gènes et nous avons évalué le chevauchement entre ces coordonnées pour retenir les gènes d’intérêt. Enfin, nous avons considéré qu’un gène était lié à une région détectée donnée lorsque les positions de la région et du gène se chevauchaient. Nous avons ensuite évalué quel gène était le plus probablement ciblé par la sélection en considérant le gène le plus proche du signal supérieur, c’est-à-dire la position de la valeur q la plus basse dans la région. Les gènes ont été annotés fonctionnellement en utilisant Uniprot (http://www.uniprot.org/), en considérant leur implication dans 30 termes enfants (c’est-à-dire les descendants directs des termes) de la catégorie « Processus biologique » (c’est-à-dire GO:0008150). Nous avons récupéré tous les termes GO correspondant à chaque gène (Données supplémentaires 4) pour 30 des 33 catégories, car nous n’avons pas pris en compte trois termes qui n’étaient pas impliqués dans les fonctions des mammifères (c’est-à-dire, GO:0006791 utilisation du soufre, GO:0006794 utilisation du phosphore, GO:0015976 utilisation du carbone). Nous avons effectué deux tests χ2 pour comparer les distributions des gènes dans les catégories GO, c’est-à-dire (i) les gènes sous sélection des régions spécifiques du genre par rapport à ceux des régions homologues, et (ii) tous les gènes sous sélection par rapport aux 18 689 gènes humains associés aux termes GO dans Swiss-Prot. Afin d’interpréter les fonctions des gènes dans un contexte d’élevage, nous avons également récupéré les informations disponibles dans la littérature sur leurs effets phénotypiques.

Enfin, pour trouver les SNP au sein des régions précédemment détectées qui étaient les plus différenciés entre les groupes sauvages et domestiques, nous avons utilisé la statistique FLK. Comme pour hapFLK, elle représente la déviation des fréquences alléliques mono-marqueur par rapport au modèle neutre estimé par la matrice de parenté65. La même procédure a été utilisée pour ajuster les scores des deux analyses à une distribution χ2 et combiner les valeurs p obtenues que celle utilisée pour le test hapFLK. Cependant, la distribution non uniforme des valeurs p a empêché l’application du cadre FDR et nous avons sélectionné les SNP dans les régions détectées avec hapFLK montrant des valeurs p <10-4. Pour ces SNP, nous avons utilisé les annotations du Variant Effect Predictor (VEP)71 qui ont été générées à partir de l’annotation du génome ovin OARv3.1 d’Ensembl v74 pour Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) et de l’annotation du génome caprin CHIR1.0 produite par le pipeline d’annotation du génome eucaryote du NCBI pour Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). Les SNP ont été classés en positions intergéniques, amont et aval (y compris les UTR), et introniques et exoniques. Les différences entre les distributions des SNP avec des p-values FLK <10-4 et tous les SNP utilisés pour détecter les signatures de sélection ont été examinées avec un test χ2.

Disponibilité des données

Les séquences et les données de métadonnées générées pour les 73 échantillons d’Ovis et 72 échantillons de Capra utilisés dans ces analyses sont disponibles publiquement. Des informations générales et tous les fichiers vcf peuvent être trouvés sur le site Web d’Ensembl (http://projects.ensembl.org/nextgen/). Tous les fichiers Fastq, les fichiers Bam et les assemblages de novo de O. orientalis et C. aegagrus peuvent être trouvés sur l’European Nucleotide Archive (https://www.ebi.ac.uk/ena) sous le code d’accession du projet Nextgen (PRJEB7436).