Probenahme

Hausschafe (O. aries) und Hausziegen (C. hircus) wurden im Iran (IROA- bzw. IRCH-Gruppe) und in Marokko (MOOA- bzw. MOCH-Gruppe) mit insgesamt 20 Tieren pro Gruppe beprobt (ergänzende Abb. 6). Diese Proben wurden zwischen Januar 2008 und März 2012 im nördlichen Teil Marokkos und zwischen August 2011 und Juli 2012 im nordwestlichen Iran im Rahmen des europäischen Nextgen-Projekts (Grant Agreement Nr. 244356) in Übereinstimmung mit den ethischen Bestimmungen der EU-Richtlinie 86/609/EWG gesammelt. Ohrclips wurden vom distalen Teil des Ohrs zufällig ausgewählter Tiere entnommen und sofort einen Tag lang in 96%igem Ethanol gelagert, bevor sie bis zur DNA-Extraktion in Kieselgelkügelchen überführt wurden.

Die wildlebenden Arten Asiatisches Mufflon (O. orientalis) und Bezoar-Steinbock (C. aegagrus) wurden im Nordwesten Irans im Rahmen der Domestizierungswiege21,22 beprobt. Dreizehn Gewebe von Asiatischen Mufflons und 18 Gewebe von Bezoar-Steinböcken (IROO- bzw. IRCA-Gruppen, siehe ergänzende Abb. 6) wurden entweder von in Gefangenschaft gehaltenen oder kürzlich erlegten Tieren sowie von gefrorenen Proben entnommen, die im iranischen Umweltministerium verfügbar waren. Dieser auf Einzelpersonen basierende Probenahmeansatz wurde entwickelt, um mögliche Verzerrungen zu minimieren, indem eine Überrepräsentation lokaler Effekte (z. B.,

Zusätzliche Daten

Zusätzlich wurde ein weltweites Rassenpanel für Schafe und Ziegen zusammengestellt (wpOA bzw. wpCH). wpOA umfasste 20 Whole-Genome Re-Sequencing (WGS)-Proben mit 12-facher Abdeckung, die 20 verschiedene weltweite Rassen repräsentieren und vom International Sheep Genome Consortium zur Verfügung gestellt wurden. wpCH bestand aus 14 WGS-Proben, die mit 12-facher Abdeckung sequenziert wurden und 9 europäische Individuen repräsentieren, d.h., 2 französische Alpenschaf- und 2 französische Saanen-Proben, die vom INRA sequenziert wurden, 5 italienische Saanen-Proben, die vom Parco Tecnologico Padano zur Verfügung gestellt wurden, und 5 australische Individuen, d. h., 2 Boer-, 2 Rangeland- und 1 Cashmere-Probe, die vom CSIRO zur Verfügung gestellt wurden (ergänzende Daten 5).

Erstellung von WGS-Daten

Genomische DNA wurde erfolgreich aus allen Gewebeproben mit dem Macherey Nagel NucleoSpin 96 Tissue Kit extrahiert, wobei das Protokoll des Herstellers angepasst wurde. Die Gewebeproben wurden in MN-Quadrat-Well-Blöcken entnommen, um 25 mg Fragmente pro Probe zu erhalten. Es wurden dreieinhalb MN-Quadrat-96-Blöcke vorbereitet, und die Extraktion wurde mit einem Tecan Freedom Evo Liquid Handler nach dem Protokoll des Herstellers durchgeführt. In einem Vorlyseschritt wurden die Proben mit 180 µl T1-Puffer und 25 µl Proteinase K über Nacht bei 56 °C homogenisiert. Zur Einstellung der Bindungsbedingungen wurden 200 µl BQ1-Puffer zugegeben und die Probenplatte 1 h bei 70 °C inkubiert; anschließend wurden 200 µl 100%iges Ethanol zugegeben. Die Lysate wurden auf die Nucleospin-Gewebebindungsplatte übertragen und ein Vakuum (-0,2 bar, 5 min) wurde angelegt, um den Durchfluss zu entfernen. Die Lysate wurden dreimal mit BW- bzw. B5-Puffer gewaschen, und der Durchfluss wurde erneut mit Vakuum entfernt. Vor der Elution der genomischen DNA wurde eine Nucleospin Tissue Binding Plate-Silikamembran 10 Minuten lang unter Vakuum bei mindestens -0,6 bar getrocknet. Der Elutionsschritt wurde mit 100 µl vorgewärmtem BE-Puffer (70 °C) und einem Zentrifugationsschritt bei 3700 rcf für 5 min in 96-PCR-Vertiefungen durchgeführt. Die genomische DNA wurde bei 4 °C gelagert, um ein Einfrieren zu vermeiden, und mit der Picogreen-Methode und einem Nanodrop auf ihre Konzentration (in ng/μl) untersucht.

Ganze Genome wurden aus 500 ng genomischer DNA, die mit dem Covaris® E210-Instrument für jede Probe auf einen Bereich von 150-700 bp geschert wurden, erneut sequenziert und für die Illumina®-Bibliotheksvorbereitung nach einem halbautomatisierten Protokoll verwendet. Die Endreparatur, das A-tailing und die Ligation mit Illumina®-kompatiblen Adaptern (BioScientific) wurden mit dem SPRIWorks Library Preparation System und dem SPRI TE Instrument (Beckmann Coulter) nach dem Protokoll des Herstellers durchgeführt. Es wurde eine Größenauswahl von 300-600 bp vorgenommen, um die meisten Fragmente zu gewinnen. Die DNA-Fragmente wurden in 12 PCR-Zyklen mit dem Platinum Pfx Taq Polymerase Kit (Life® Technologies) und Illumina®-Adapter-spezifischen Primern amplifiziert. Die Bibliotheken wurden mit 0,8x AMPure XP-Beads (Beckmann Coulter) gereinigt und mit dem Agilent 2100 Bioanalyzer (Agilent® Technologies) und qPCR-Quantifizierung analysiert. Die Bibliotheken wurden unter Verwendung der 100-Basen-Lesechemie in einer Paired-End-Fließzelle auf dem Illumina® HiSeq2000 sequenziert.

Die Paired-End-Reads von Illumina wurden für Ovis auf das Schaf-Referenzgenom (OAR v3.1, GenBank-Assembly GCA_000298735.146) und für Capra auf das Ziegen-Referenzgenom (CHIR v1.0, GenBank-Assembly GCA_000317765.147) unter Verwendung von BWA-MEM48 gemappt. Die für jedes Individuum erstellte BAM-Datei wurde mit Picard SortSam sortiert und mit Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator und GATK IndelRealigner49 sowie SAMtools calmd50 verbessert.

Die Entdeckung von Varianten wurde mit drei verschiedenen Algorithmen durchgeführt: Samtools mpileup50, GATK UnifiedGenotyper51, und Freebayes52. Variant Sites wurden unabhängig für jede der sechs Gruppen identifiziert, wobei die Multi-Sample-Modi der Aufrufalgorithmen verwendet wurden: (i) 162 Proben von MOOA; (ii) 20 Proben von IROA; (iii) 14 Proben von IROO; (iv) 162 Proben von MOCH; (v) 20 Proben von IRCH; (vi) 19 Proben von IRCA. Für einige Gruppen waren im Rahmen des NextGen-Projekts (siehe oben) die WGS von mehr Individuen verfügbar. Die in der vorliegenden Studie verwendeten Proben wurden so ausgewählt, dass nach Möglichkeit ausgewogene Gruppen von 20 Individuen erhalten wurden. Für die IRCA- und IROO-Gruppen standen zu einem späteren Zeitpunkt zusätzliche Proben zur Verfügung, die für nachgeschaltete Analysen hinzugefügt wurden. Tiere mit niedriger Alignment- und Calling-Qualität wurden entfernt, um den endgültigen Datensatz zu erhalten (Supplementary Data 5).

In jeder Gruppe gab es zwei aufeinander folgende Runden der Filterung der Qualität der Variantenstellen. In der Filterstufe 1 wurden die Anrufe aus den drei Algorithmen zusammengeführt, während die Anrufe mit der niedrigsten Konfidenz herausgefiltert wurden. Ein Variantenort wurde akzeptiert, wenn er von mindestens zwei verschiedenen Algorithmen mit einer phred-Variantenqualität >30 aufgerufen wurde. Ein alternatives Allel an einer Stelle wurde akzeptiert, wenn es von einem der Algorithmen genannt wurde und die Genotypenzahl >0 war. In der Filterstufe 2 wurde die Rekalibrierung des Variant Quality Score durch GATK verwendet. Zunächst wurde ein Trainingsset der Variantenstellen mit der höchsten Konfidenz innerhalb der Gruppe generiert, bei dem (i) die Stelle von allen drei Varianten-Callern mit einer phred-skalierten Variantenqualität >100 genannt wird, (ii) die Stelle biallelisch ist, (iii) die Minor-Allel-Zahl mindestens 3 beträgt, wobei nur Proben mit einer phred-skalierten Genotypqualität >30 gezählt werden. Der Trainingssatz wurde zur Erstellung eines Gauß-Modells mit dem Tool GATK VariantRecalibrator unter Verwendung der folgenden Variantenannotationen von UnifiedGenotyper verwendet: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Ein Gaußsches Modell wurde auf den gesamten Datensatz angewandt und ein VQSLOD (log odds ratio of being a true variant) erzeugt. Standorte wurden gefiltert, wenn der VQSLOD <Cutoff-Wert war. Der Cutoff-Wert wurde für jede Gruppe wie folgt festgelegt: Minimaler VQSLOD = {der Medianwert des VQSLOD für die Varianten des Trainingssatzes}-3 × {die mediane absolute Abweichung des VQSLOD der Varianten des Trainingssatzes}. Das Verhältnis von Übergangs- zu Transversions-SNPs deutet darauf hin, dass das gewählte Cutoff-Kriterium das beste Gleichgewicht zwischen Selektivität und Sensitivität bietet.

SNPs-Call-Sets wurden für sechs Gruppen von Ovis- und Capra-Tieren generiert (d. h. iranische und marokkanische domestizierte Tiere sowie Wildtiere für jede Gattung). Da die in dieser Studie durchgeführten Analysen Vergleiche zwischen den Gruppen erforderten, erstellten wir Genotyp-Call-Sets an einem einheitlichen Satz von SNP-Stellen für alle Tiere aus jeder Gruppe. Für jede Gattung fügten wir die Variantenaufrufe aus ihren drei Gruppen zusammen und behielten nur biallelische Positionen ohne fehlende Daten bei. Die Genotypen wurden an jeder biallelischen SNP-Stelle für alle interessierenden Individuen von GATK UnifiedGenotyper unter Verwendung der Option GENOTYPE_GIVEN_ALLELES neu aufgerufen. In dieser Phase wurde die Liste der Individuen um die Tiere erweitert, die zu den Weltrassenpanels für Schafe und Ziegen (wpOA und wpCH) gehören, sowie um zusätzliche Wildproben, die in dieser Phase verfügbar wurden (4 O. orientalis und 4 C. aegagrus). Die Genotypen wurden innerhalb der Gruppen durch Beagle 453 verbessert und gestaffelt und dann herausgefiltert, wenn die Wahrscheinlichkeit eines Genotyps weniger als 0,95 betrug. Schließlich wurden die Stellen herausgefiltert, die in den verschiedenen Untergruppen von Individuen, die in dieser Studie verwendet wurden, monomorph waren (siehe unten).

Um die zwischen Ovis und Capra festgestellten Selektionssignale zu vergleichen, führten wir ein Cross-Alignment zwischen den beiden Referenzgenomen durch. Zunächst verwendeten wir die paarweise Alignment-Pipeline aus der Ensembl Release 69 Code Base54 , um die Referenzgenome von Schaf (OARv3.1) und Ziege (CHIR1.0) zu alignieren. Diese Pipeline verwendet LastZ55 für das Alignment auf DNA-Ebene, gefolgt von einer Nachbearbeitung, bei der die alignierten Blöcke entsprechend ihrer Position in beiden Genomen aneinandergereiht werden. Die LastZ-Pipeline zum paarweisen Alignment wird von Ensembl routinemäßig für alle unterstützten Arten ausgeführt, aber die Ziege ist noch nicht in Ensembl enthalten. Um eine Verzerrung zugunsten einer der beiden Arten zu vermeiden, haben wir zwei verschiedene interspezifische Alignments erstellt. Bei dem einen wurde das Schaf als Referenzgenom und die Ziege als Nicht-Referenz verwendet, bei dem anderen die Ziege als Referenzgenom und das Schaf als Nicht-Referenz. Der Unterschied besteht darin, dass genomische Regionen der Referenzspezies gezwungen sind, eindeutig auf einzelne Loci der Nicht-Referenzspezies zu kartieren, während genomische Regionen der Nicht-Referenzspezies auf mehrere Stellen der Referenzspezies kartiert werden können. Wir erhielten für Chromosomenabschnitte eines Referenzgenoms die Koordinaten auf dem Nicht-Referenzgenom. Schließlich verwendeten wir für die SNPs, die in einer Gattung entdeckt wurden, das Alignment des gesamten Genoms mit dem Referenzgenom der anderen Gattung, um die entsprechenden Positionen zu identifizieren (ergänzende Tabelle 6).

Genetische Struktur

Um die genetische Vielfalt innerhalb der Gruppen zu beschreiben, verwendeten wir VCFtools56 , um zusammenfassende Statistiken über die genetische Variation der 73 Individuen von Ovis zu berechnen (d. h., 13 IROO, 20 IROA, 20 MOOA und 20 wpOA) und 72 Individuen für Capra (d. h. 18 IRCA, 20 IRCH, 20 MOCH und 14 wpCH). Die gemessenen Statistiken waren die Gesamtzahl der polymorphen Varianten (S) für die Gesamtheit der Individuen in jeder Gattung und innerhalb jeder Gruppe, die gemittelte Nukleotiddiversität (π) innerhalb jeder Gruppe und der Inzuchtkoeffizient (F) für jedes Individuum. Innerhalb jeder Gattung wurden die Unterschiede zwischen der wilden Gruppe und jeder Hausgruppe mit einem einseitigen t-Test für individuelle Inzucht- und genetische Belastungswerte und einem zweiseitigen Mann-Whitney-Test für die Nukleotiddiversität pro Standort getestet.

Die Gesamtdivergenz zwischen den vier Gruppen innerhalb jeder Gattung (d. h., Wildform, iranische und marokkanische domestizierte Arten und Weltpanel) wurde anhand aller biallelischen SNPs und der durchschnittlichen gewichteten paarweisen Fst nach Weir und Cockerham57 geschätzt, wie in VCFtools56 implementiert. Die genetische Struktur zwischen den Gruppen wurde mit der Clustering-Methode sNMF26 bewertet, nachdem der Datensatz mit VCFtools bereinigt wurde, um SNPs mit einem Kopplungsungleichgewicht (r²) von mehr als 0,2 zu entfernen. Das Kopplungsungleichgewicht (r²) wurde zwischen SNP-Paaren innerhalb gleitender Fenster von 50 SNPs berechnet, wobei ein SNP pro Paar zufällig entfernt wurde, wenn r² größer als 0,2 war. Für jede sNMF-Analyse wurden fünf Durchläufe mit der gleichen Anzahl von Clustern (K) mit Werten von K zwischen 1 und 10 durchgeführt. Wir haben das Kriterium der Kreuzentropie verwendet, um die wahrscheinlichste Clusterlösung zu ermitteln. Es wurden jedoch auch alternative Partitionen für verschiedene Zahlen von K untersucht, um zu bewerten, wie die Individuen zwischen den Clustern aufgeteilt wurden.

Um zwischen gemeinsamer Abstammung und Vermischung zu unterscheiden, haben wir TreeMix27 ausgeführt, um gemeinsam Populationssplits und nachfolgende Vermischungsereignisse zu schätzen, wobei der für sNMF verwendete beschnittene Datensatz verwendet wurde. Wir führten TreeMix mit der Option -global aus, um unsere Maximum-Likelihood-Inferenzen zu verfeinern. Wir verwurzelten den TreeMix-Baum mit der Aufteilung zwischen wilden und heimischen Individuen. Die Blockgröße für das Jackknifing betrug -k 500 SNPs, was ungefähr 150 kb entspricht und die durchschnittlichen LD-Blöcke, die sowohl bei Schafen als auch bei Ziegen gefunden wurden, übersteigt. Wir erstellten einen Maximum-Likelihood-Baum ohne Migration und fügten dann Migrationsereignisse hinzu und untersuchten die inkrementelle Änderung der durch das Modell erklärten Varianz und die Restwerte zwischen den Individuen. Ziel war es, einen potenziell hohen Restwert oder eine Migrationskante zwischen wilden und domestizierten Individuen zu erkennen. Um die statistische Relevanz möglicher, von TreeMix identifizierter Vermischungsvektoren weiter zu erforschen (ergänzende Tabelle 3), berechneten wir den Drei-Populations-Test f328 als formalen Test auf genetische Introgression, indem wir das qp3Pop-Programm der ADMIXTOOLS-Suite58 für jede Kombination von Gruppen verwendeten. Für Capra wurde die wpCH-Gruppe in australische, französische und italienische Rassen aufgeteilt. Die Ergebnisse sind in den ergänzenden Daten 2 aufgeführt.

Demografische Inferenz

Für jede Gattung führten wir Analysen zur demografischen Inferenz der Vorfahren mit Hilfe des MSMC-Modells durch, das in der MSMC2-Software25 implementiert ist. MSMC basiert auf dem paarweisen sequentiellen Markovianischen Koaleszenzmodell59, verwendet jedoch Haplotypen aus phasenweisen Genomsequenzdaten als Input. Für jede Analyse wurden zwei Individuen aus einer Gruppe verwendet, also 4 Haplotypen. Jede Analyse wurde für einen anderen zufälligen Satz von zwei Individuen wiederholt, d. h. ein Replikat der Analyse pro Gruppe. Die Eingabe- und Ausgabedateien wurden mit den Python-Skripten erstellt und analysiert, die mit der MSMC-Software geliefert wurden und unter https://github.com/stschiff/msmc-tools zu finden sind. Die Analyseparameter wurden als Standardwerte beibehalten, mit Ausnahme der Mutationsrate, die auf 2,5×10-8 und die Generationslänge auf 2 Jahre festgelegt wurde. Um die Unsicherheit der Zeitschätzungen abzuschätzen, variierten wir diese Parameter (Mutationsrate von 2,5×10-8 und 1,0×10-8 in Kombination mit einer Generationslänge von 2 und 4 Jahren) und lieferten eine grobe Schätzung des Domestikationszeitraums (siehe ergänzende Abb. 2).

Genetische Belastung

Die genetische Belastung wurde auf zwei Arten geschätzt. Erstens durch Berechnung der genetischen Belastung für jedes Individuum als Summe der schädlichen Fitness-Effekte über alle Protein-kodierenden genomischen Positionen nach der Methode von Librado et al.60. Kurz gesagt, als Ersatz für evolutionäre Einschränkungen verwendeten wir die PhyloP-Scores aus dem 46-fachen Säugetier-Alignment (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Anhand dieses Alignments identifizierten wir proteinkodierende Stellen, die sich unter funktionalen Einschränkungen entwickeln (PhyloP-Score ≥1,5). Für jedes Ovis- oder Capra-Genom haben wir dann untersucht, ob diese Stellen mutiert waren. Wenn dies der Fall war, wurden die phyloP-Scores über alle mutierten Stellen summiert, so dass Mutationen an stark eingeschränkten Stellen proportional mehr zur Schätzung der Gesamtbelastung beitragen. Auf diese Weise erhielten wir eine Belastungsschätzung für jedes Schaf-/Ziegengenom. Um schließlich eine durchschnittliche Belastung pro Standort zu erhalten, wurde diese durch die Gesamtzahl der analysierten Positionen geteilt. Es ist erwähnenswert, dass wir uns auf homozygote Stellen beschränkt haben, um die Modellierung des Dominanzkoeffizienten von Mutationen an heterozygoten Stellen (z. B. rezessiv, intermediär, dominant) zu vermeiden. Zweitens verglichen wir die genetische schädliche Belastung in wilden und domestizierten Ovis-Gruppen mit einem Wilcoxon-Test, wobei die Alternativhypothese lautete, dass die domestizierten Tiere eine höhere Belastung aufweisen als ihre wilden Verwandten. p-Werte wurden für Mehrfachtests61 korrigiert, und wir wendeten einen Schwellenwert von angepassten p-Werten < 0,05 an. Wir führten eine Gen-Ontologie-Anreicherungsanalyse für den Satz von Genen durch, die eine signifikante Zunahme der genetischen Belastung zeigten, und verwendeten dazu WebGestalt62,63. Da die Referenzgenome für Gene nur unzureichend annotiert sind, stützten wir uns auf orthologe Einzelkopien zwischen unserer Spezies und Mensch und Maus. Gene vom X-Chromosom wurden aus dem Hintergrundsatz ausgeschlossen. Wir haben diese Analyse bei Capra aufgrund der höheren Inzucht, die in den Wildproben beobachtet wurde, nicht durchgeführt.

Erkennung von Selektionssignaturen

Für die Erkennung von Selektionssignaturen im Zusammenhang mit der Domestizierung haben wir alle biallelischen SNPs verwendet, die eine Minor-Allel-Häufigkeit von mehr als 0,10 in mindestens einer der drei getesteten Gruppen aufweisen (d. h. iranische und marokkanische Domestizierungsgruppen und die Wildgruppe für jede Gattung). Da wir davon ausgingen, dass Selektionsmerkmale, die mit dem Domestikationsprozess zusammenhängen, bei allen Haustieren vorhanden sind, verfolgten wir die folgende allgemeine Strategie: Wir testeten mit hapFLK29 (siehe ergänzende Anmerkung 5 und ergänzende Abbildungen 9, 10 und 11) für jede Gattung die wilde Gruppe gegen jede der traditionell verwalteten Haustiergruppen (d. h. die iranische und die marokkanische) und konzentrierten uns auf die gemeinsamen Regionen, die vermutlich unter Selektion stehen und in beiden Fällen entdeckt wurden. Die Stichprobengröße der Gruppen (n = 13-20) entsprach den Anforderungen der Methode29. Wir überprüften visuell, ob die mit hapFLK gefundenen konsistenten Selektionssignaturen auch in den entsprechenden Welt-Panels jeder Gattung vorhanden waren, nahmen diese Gruppen jedoch aufgrund ihrer Multi-Breed-Zusammensetzung nicht in den statistischen Test auf. Schließlich haben wir mit einem stratifizierten FDR-Ansatz nach gemeinsamen Selektionssignalen zwischen Ovis und Capra gesucht. Die Strategie ist in der ergänzenden Abb. 4 dargestellt.

Wir führten hapFLK-Tests durch, um die wilde Gruppe mit jeder der iranischen und marokkanischen Gruppen in jeder Gattung zu vergleichen. Die Verwandtschaftsmatrix wurde aus den Reynoldschen genetischen Abständen64 zwischen Gruppenpaaren berechnet, wobei eine zufällige Teilmenge von einem Prozent der Varianten verwendet wurde. Der abgeleitete Populationsbaum wurde mithilfe des Nachbarschaftsverbindungsalgorithmus erstellt. Für jeden SNP führten wir den hapFLK-Test durch, der haplotypische Informationen einbezieht, um die Aussagekraft zur Erkennung selektiver Sweeps zu erhöhen. Für jeden getesteten SNP berechnete die hapFLK-Statistik die Abweichung der haplotypischen Häufigkeiten in Bezug auf das neutrale Modell, das durch die Verwandtschaftsmatrix65 geschätzt wurde. Zur Nutzung von Informationen über Kopplungsungleichgewichte verwendet hapFLK das Mehrpunktmodell von Scheet und Stephens66 für Multilocus-Genotypen, das an unphasierte Daten angepasst werden kann. Eine der Hauptanwendungen dieses Modells ist die Durchführung von Phasenschätzungen (fastPHASE-Software66). In unserer Analyse wurde das Modell an unphasierten Daten trainiert, so dass unsere Analyse die Phasenunsicherheit berücksichtigt. Die Methode wurde verwendet, um lokale Haplotypen entlang der Chromosomen in einer bestimmten Anzahl von Clustern K zu gruppieren, die auf 25 festgelegt wurde, und zwar unter Verwendung eines Hidden-Markov-Modells.

Um die gemeinsamen Regionen zu identifizieren, die in den beiden traditionell verwalteten einheimischen Gruppen für jede Gattung vermutlich der Selektion unterliegen, haben wir die beiden vorherigen hapFLK-Analysen kombiniert. Für jede Analyse wurden die hapFLK-Werte an eine χ2-Verteilung angepasst, um p-Werte zu erhalten (Skript verfügbar unter https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Die Ergebnisse der beiden Kontraste zwischen der Wildtiergruppe und jeder der Haustiergruppen wurden mit der Stouffer-Methode67 kombiniert, um einzelne p-Werte für den Vergleich von Wild- und Haustieren zu erhalten. Schließlich wurde der FDR-Rahmen68 auf den gesamten SNP-Satz angewendet, um die kombinierten p-Werte in q-Werte umzuwandeln. SNPs mit q-Werten < 10-2 wurden beibehalten und in genomischen Regionen gruppiert, wenn sie weniger als 50 kb voneinander entfernt waren.

Um zu untersuchen, ob das Selektionssignal zwischen Ovis und Capra geteilt wurde, verwendeten wir zunächst das Cross Alignment der beiden Referenzgenome, um homologe Segmente zu identifizieren. Anschließend wendeten wir ein stratifiziertes FDR-Verfahren an69. Dieser Ansatz basiert auf der Tatsache, dass die Tests angesichts der Vorabinformationen in den genetischen Daten69 eine inhärente Stratifizierung aufweisen, da die zugrunde liegende Verteilung der wahren alternativen Hypothesen je nach der unterschiedlichen Dynamik der verschiedenen Genomregionen unterschiedlich sein kann, was zu unterschiedlichen Verteilungen der p-Werte führt. Dies erfordert, dass die FDR-bereinigten p-Werte (d. h. q-Werte) für die verschiedenen Schichten getrennt ermittelt werden. Wir suchten nach Konvergenzen in jeder Gattung, indem wir die Regionen, die homolog zu denen sind, die in der anderen Gattung entdeckt wurden (als gemeinsames Stratum bezeichnet), und den Rest des Genoms (als allgemeines Stratum bezeichnet) trennten. Wir extrahierten die p-Werte separat für jede der beiden definierten Schichten und berechneten dann die q-Werte mithilfe des FDR-Rahmens. Diese stratifizierten q-Werte waren die endgültigen Größen, die für die statistische Signifikanz (<10-2) in Betracht gezogen wurden, um SNPs unter Selektion zu erkennen und sie in die entsprechenden Genomregionen einzubinden.

Um auf konvergente Selektionssignaturen zu testen, die Wild- und Haustiere in beiden Gattungen unterscheiden, untersuchten wir die Beziehung zwischen der Signifikanzschwelle, die auf die q-Werte (die wir von 0,2 bis 0,002 variieren ließen) in einer Gattung angewandt wurde, und der geschätzten Wahrscheinlichkeit, dass ein SNP im gemeinsamen Stratum der anderen Gattung selektiert wird, unter Verwendung des Storey et al.70 Ansatzes. Ein Anstieg der abgeleiteten Wahrscheinlichkeit mit einer Verringerung des auf den q-Wert angewendeten Schwellenwerts (Erhöhung der Stringenz) deutet darauf hin, dass je signifikanter die Region in einer Gattung ist, desto wahrscheinlicher ist es, dass wir signifikante SNPs in der anderen Gattung finden.

Wir filterten die Selektionssignale heraus, die in den drei einheimischen Gruppen nicht konsistent waren. Für jede entdeckte Region verwendeten wir die phasierten Haplotypen von Individuen, die mit Hilfe von Neighbor-Joining-Bäumen auf der Grundlage der prozentualen Identität zwischen den Sequenzen geclustert wurden. Nur Regionen, die konsistente Signale zeigten, wurden beibehalten (ergänzende Abb. 5).

Um herauszufinden, ob die mit hapFLK entdeckten Selektionssignale auf eine Lockerung der Selektion oder eine positive Selektion bei den Haustieren hinweisen, schätzten wir den Unterschied in der Nukleotiddiversität (π) in jeder mutmaßlichen Selektionsregion zwischen der Wild- und der Haustiergruppe. Wir drückten diesen Unterschied als Δπ-Index aus, der für jede genomische Region als Differenz zwischen dem für die wilde Gruppe berechneten π und dem Durchschnitt von π für die iranische und marokkanische domestizierte Gruppe abzüglich der über das gesamte Genom berechneten Differenz in π zwischen diesen beiden Gruppen berechnet wurde:

Ein negativer Wert würde darauf hinweisen, dass die Nukleotidvielfalt in der wilden Gruppe im Vergleich zum Durchschnitt der beiden einheimischen Gruppen geringer ist, und würde als Anzeichen für eine nachlassende Selektion in den letztgenannten Gruppen, eine diversifizierende Selektion in den einheimischen Gruppen oder eine positive Selektion in den wilden Gruppen angesehen werden. Im Gegensatz dazu würde ein positiver Wert auf eine gerichtete positive oder stabilisierende Selektion hinweisen, die in den einheimischen Gruppen stattgefunden hat. Außerdem haben wir die Haplotyp-Clustering-Methode verwendet, um in jeder Region manuell zu überprüfen, ob der festgestellte Selektionsschub die durch den Δπ-Index gegebenen Hinweise bestätigt.

Wir haben die funktionalen Interpretationen wie folgt durchgeführt. Für jede Region, die selektiert wurde, betrachteten wir die Region plus 50 kb auf jeder Seite, um funktionelle Rollen zu identifizieren, und 5 kb stromaufwärts und stromabwärts von Genen, und wir bewerteten die Überlappung zwischen diesen Koordinaten, um die Gene von Interesse zu erhalten. Schließlich gingen wir davon aus, dass ein Gen mit einer bestimmten erkannten Region verwandt ist, wenn sich die Positionen der Region und des Gens überschneiden. Anschließend wurde ermittelt, welches Gen am ehesten von der Selektion betroffen war, indem das dem Top-Signal am nächsten liegende Gen, d. h. die Position mit dem niedrigsten q-Wert innerhalb der Region, berücksichtigt wurde. Die Gene wurden mit Hilfe von Uniprot (http://www.uniprot.org/) funktionell annotiert, indem ihre Beteiligung an 30 untergeordneten Begriffen (d. h. den direkten Ableitungen der Begriffe) der Kategorie „Biologischer Prozess“ (d. h. GO:0008150) berücksichtigt wurde. Wir haben alle GO-Terme, die jedem Gen entsprechen (ergänzende Daten 4), für 30 der 33 Kategorien abgerufen, da wir drei Terme nicht berücksichtigt haben, die nicht an Säugetierfunktionen beteiligt sind (d. h. GO:0006791 Schwefelverwertung, GO:0006794 Phosphorverwertung, GO:0015976 Kohlenstoffverwertung). Wir führten zwei χ2-Tests durch, um die Verteilung der Gene in den GO-Kategorien zu vergleichen, d.h. (i) Gene unter Selektion aus gattungsspezifischen Regionen im Vergleich zu denen aus homologen Regionen und (ii) alle Gene unter Selektion im Vergleich zu den 18.689 menschlichen Genen, die mit GO-Begriffen in Swiss-Prot assoziiert sind. Um die Funktionen von Genen im Kontext von Nutztieren zu interpretieren, haben wir auch die in der Literatur verfügbaren Informationen über ihre phänotypischen Auswirkungen abgerufen.

Um schließlich die SNPs innerhalb der zuvor entdeckten Regionen zu finden, die sich am stärksten zwischen Wild- und Haustiergruppen unterschieden, haben wir die FLK-Statistik verwendet. Wie bei hapFLK handelt es sich um die Abweichung der Allelhäufigkeiten der Einzelmarker vom neutralen Modell, das durch die Verwandtschaftsmatrix65 geschätzt wird. Es wurde dasselbe Verfahren verwendet, um die Ergebnisse der beiden Analysen an eine χ2-Verteilung anzupassen und die erhaltenen p-Werte zu kombinieren, wie es für den hapFLK-Test verwendet wurde. Die ungleichmäßige Verteilung der p-Werte schloss jedoch die Anwendung des FDR-Rahmens aus, und wir wählten SNPs innerhalb der mit hapFLK entdeckten Regionen mit p-Werten <10-4 aus. Für diese SNPs verwendeten wir die Variant Effect Predictor (VEP) Annotationen71 , die aus der Ensembl v74 Schaf OARv3.1 Genom Annotation für Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) und aus der Ziege CHIR1.0 Genom Annotation, die von der NCBI eukaryotischen Genom Annotation Pipeline für Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/) erstellt wurde. Die SNPs wurden als intergene, vor- und nachgeschaltete (einschließlich UTRs) sowie intronische und exonische Positionen klassifiziert. Die Unterschiede zwischen den Verteilungen der SNPs mit FLK p-values <10-4 und allen SNPs, die für die Erkennung von Selektionssignaturen verwendet wurden, wurden mit einem χ2-Test untersucht.

Datenverfügbarkeit

Die für die 73 Ovis- und 72 Capra-Proben, die in diesen Analysen verwendet wurden, generierten Sequenzen und Metadaten sind öffentlich verfügbar. Allgemeine Informationen und alle vcf-Dateien sind auf der Ensembl-Website zu finden (http://projects.ensembl.org/nextgen/). Alle Fastq-Dateien, Bam-Dateien und de novo-Zusammenstellungen von O. orientalis und C. aegagrus sind im European Nucleotide Archive (https://www.ebi.ac.uk/ena) unter dem Zugriffscode des Nextgen-Projekts (PRJEB7436) zu finden.