Probenahme
Wir analysierten (i) genomweite DNA-Daten, die aus Zähnen des Schädels (MCE1356, Fig. 2b, ergänzende Abb. S1) und aus Gewebeproben von acht Belugas und acht Narwalen, die in der Diskobucht, Westgrönland, entnommen wurden, und (ii) stabile Kohlenstoff- und Stickstoff-Isotopenzusammensetzungen von Knochenkollagen aus MCE1356 und einer Referenzgruppe von 18 Belugas und 18 Narwalen aus Westgrönland (Abb. 1c). Die Gewebeproben (Haut) wurden in 96%igem Ethanol gelagert. Die Proben wurden von Wissenschaftlern des Greenland Institute of Natural Resources im Rahmen der allgemeinen Genehmigung der grönländischen Regierung für biologische Probenahmen bei den Inuit gesammelt und mit den CITES-Genehmigungen IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 und 116.2006 nach Dänemark exportiert. Die Probeninformationen sind in der ergänzenden Tabelle S2 aufgeführt.
DNA-Extraktion, Bibliotheksvorbereitung und Sequenzierung
Gewebeproben
DNA aus Gewebeproben wurde mit dem Qiagen Blood and Tissue Kit gemäß dem Herstellerprotokoll extrahiert. Die DNA wurde mit dem Covaris M220 Focused-Ultraschallgerät fragmentiert, um Fragmente mit einer Länge von ~350-550 Basenpaaren (bp) zu erzeugen. Die Bibliotheken wurden aus den fragmentierten DNA-Extrakten mit Illumina NeoPrep unter Verwendung der Standardeinstellungen gemäß dem NeoPrep Library Prep System Guide erstellt. PCR-Amplifikation, Quantifizierung und Normalisierung wurden alle mit dem NeoPrep Library Prep System durchgeführt. Die Bibliotheken wurden auf einem Agilent 2100 Bioanalyzer auf ihre Größenverteilung hin überprüft und in äquimolaren Verhältnissen gepoolt, bevor sie auf einem Illumina HiSeq 2500 mit 80 bp SE-Technologie sequenziert wurden.
Zahn-/Knochenproben
Im Gegensatz zu den Gewebeproben weisen alte Knochen- und Zahnproben relativ geringe DNA-Konzentrationen auf, was andere Extraktions- und Bibliotheksaufbauprotokolle erfordert. Etwa 0,5 g Knochenpulver von fünf Zähnen und ein Knochensplitter der Probe MCE1356 wurden mit einem handgeführten Dremel gebohrt. Die DNA wurde aus dem Zahn-/Knochenpulver in einem speziellen Labor für antike DNA am Naturhistorischen Museum Dänemarks, Universität Kopenhagen, unter Verwendung des in11 beschriebenen Extraktionspuffers und unter Hinzufügung einer 30-minütigen Vorverdauungsphase12 extrahiert. Anstatt Zymo-Spin V Säulen (Zymo Research) zu verwenden, wurde der Extraktionspuffer mit Amicon Ultra 30 kDa Zentrifugalfiltereinheiten konzentriert und mit Qiagen Minelute Röhrchen weiter konzentriert und gereinigt. Die gereinigten Extrakte wurden dann nach dem in13 beschriebenen Protokoll in Illumina-Bibliotheken eingebaut. Wir haben qPCR verwendet, um zu überprüfen, ob der Bibliotheksaufbau erfolgreich war, um auszuwählen, welche Bibliotheken sequenziert werden sollen, und um die entsprechende Anzahl von PCR-Zyklen zu berechnen, die erforderlich sind, um jede Bibliothek ausreichend zu amplifizieren, ohne eine Überamplifikation zu verursachen. Insgesamt wurden vier Bibliotheken mit einzigartigen 6 bp-Indizes amplifiziert und auf der Illumina MiSeq-Plattform mit 250 bp-SE-Sequenzierung auf endogene Inhalte untersucht. Die besten Bibliotheken wurden dann auf der Illumina HiSeq 2500-Plattform unter Verwendung der 80 bp SE-Technologie erneut sequenziert.
Bioinformatische Analyse
Alle Kartierungs- und DNA-Schadensanalysen wurden mit der Paleomix-Pipeline 1.2.1214 durchgeführt. Die Reads wurden mit AdapterRemoval 2.2.015 unter Verwendung der Standardeinstellungen getrimmt, mit Ausnahme der minimalen Leselänge, die auf 25 bp gesetzt wurde. Die Reads wurden mit FastQC inspiziert und mit BWA16 unter Anwendung des Backtrack-Algorithmus und Deaktivierung der Start-Seed-Länge ausgerichtet. Wenn Reads auf mehrere Positionen kartiert wurden oder einen Mapping-Qualitätsscore (MAPQ-Score von BWA) von weniger als 30 hatten, wurden sie mit SAMtools17 entfernt. Sequenzduplikate wurden mit MarkDuplicates von Picard (verfügbar unter: http://broadinstitute.github.io/picard) entfernt und das endgültige Alignment wurde mit GATK18 um Indels herum neu ausgerichtet. Die Deaminierung von Cytosin zu Uracil in der Probe MCE1356 wurde anhand der Ausgabe von mapDamage v2.0.619 bewertet. Die Ergebnisse von mapDamage zeigten kein klares Signal für eine Deaminierung in den Sequenzen (ergänzende Abb. S3).
Mitochondriale Analyse
Um die mütterliche Abstammung von MCE1356 zu bestimmen, wurden die DNA-Sequenzierungs-Reads auf die mitochondrialen Referenzgenome von Beluga (GenBank-Zugang: KY444734) und Narwal (GenBank-Zugang: NC_005279) abgebildet und die durchschnittliche Abdeckung verglichen. Die Reads der acht Beluga- und acht Narwalproben, die das Referenzpanel bildeten, wurden auf die jeweiligen mitochondrialen Referenzgenome abgebildet. Mit BEDtools20, SAMtools17 und GATK18 wurden mitochondriale Konsensussequenzen von MCE1356 und den 16 Proben des Referenzpanels erstellt, deren Regionen von mehr als fünf Reads abgedeckt wurden. Wir erstellten zwei Sequenzalignments mit ClustalW21 unter Verwendung der Standardeinstellungen, die die 16 Proben des Referenzpanels und entweder die Version von MCE1356, die dem mitochondrialen Referenzgenom des Belugas zugeordnet ist, oder die Version von MCE1356, die dem mitochondrialen Referenzgenom des Narwals zugeordnet ist, enthielten. Aus den beiden Alignments wurden mit Hilfe von Popart 1.723 (verfügbar unter: http://popart.otago.ac.nz) Median-Join-Haplotypen-Netzwerke22 konstruiert, wobei Standorte mit Indels oder fehlenden Daten ausgeschlossen wurden. Anschließend wurden sowohl die Post-Mapping-Abdeckung als auch die beiden Haplotyp-Netzwerke verwendet, um die Spezies von MCE1356′s mütterlicher Abstammung zu bestimmen.
Nukleare DNA-Analysen
Wir haben die DNA-Sequenzierungs-Reads von allen Proben auf das Referenzgenom des Schwertwals (Orcinus orca) (GCA_000331955.2) gemappt. Vor kurzem wurde ein hochwertiges Beluga-Genom veröffentlicht24, aber die Zuordnung zu einer der beiden potenziellen Elternarten könnte zu Verzerrungen bei den Analysen führen. Daher kartierten wir die Reads auf das Genom des Schwertwals, da es gut assembliert ist und Schwertwale immer noch relativ eng mit Belugas und Narwalen verwandt sind (Divergenzzeit von 12 MYA)25, jedoch weit genug entfernt, um das Risiko einer Introgression zu minimieren, die unsere Analysen erschweren würde.
Für alle weiteren Filterungen verwendeten wir ANGSD v0.92326, ein Softwarepaket, das Genotyp-Likelihoods anstelle von aufgerufenen Genotypen verwendet, was besonders bei der Analyse von NGS-Daten mit geringer Abdeckung nützlich ist. Wir verwendeten die in ANGSD implementierte SAMtools17 -Methode zur Schätzung der Genotyp-Likelihoods und leiteten die Haupt- und Nebenallele direkt aus den Genotyp-Likelihoods ab, indem wir einen Maximum-Likelihood-Ansatz wie in27 beschrieben verwendeten.
Um die gemappten Daten zu visualisieren, zeichneten wir die Lesetiefenverteilung aller Individuen auf, wobei Sites mit mehr als zwei Allelen und Sites mit einem Phred-Score unter 25 ausgeschlossen wurden. Die Visualisierung der Daten aller Individuen zusammen ermöglichte es uns, die mittlere Lesetiefe von 4,14x zu schätzen und Stellen mit erhöhter Lesetiefe zu identifizieren. Solche Stellen stammten mit größerer Wahrscheinlichkeit von Paralogen und repetitiven Regionen des Genoms. Der Datensatz wurde visuell inspiziert und weiter gefiltert, wobei alle Stellen mit einer Lesetiefe von mehr als neun (6,9 % der Stellen) ausgeschlossen wurden. In ANGSD wurden die SNPs auf der Grundlage ihrer Allelhäufigkeit bestimmt. Die Minor-Allel-Häufigkeit (MAF) wurde aus den Genotyp-Likelihoods geschätzt, und ein Likelihood-Ratio-Test wurde verwendet, um festzustellen, ob die MAF von Null verschieden war. Wenn der p-Wert des Likelihood-Ratio-Tests <1e-4 war, wurde der Standort als polymorph betrachtet und im Datensatz beibehalten. Die Anwendung dieses SNP-Filters bedeutete, dass keine Sites mit weniger als vier Reads beibehalten wurden, da Sites, die von weniger Reads abgedeckt werden, keine so niedrigen p-Werte erhalten können. Diese Filter wurden auf einen Datensatz mit allen 17 Proben und einen Datensatz ohne MCE1356 angewandt.
Um festzustellen, ob es sich bei MCE1356 um einen Beluga/Narwal-Hybriden handelt, wurde der Datensatz mit 17 Individuen weiter gefiltert, wobei Stellen ohne Reads in MCE1356 ausgeschlossen wurden. Zu diesem Zeitpunkt betrug die mittlere Lesetiefe von MCE1356 nur das 1,15-fache. Um sicherzustellen, dass wir keine Multikopie-Loci analysierten, schlossen wir Stellen aus, die von mehr als einem Read in MCE1356 abgedeckt wurden.
Unser Ziel war es, die in MCE1356 gefundenen Allele mit den Allelen im Referenzpanel zu vergleichen, daher schätzten wir die Wahrscheinlichkeit ab, das in MCE1356 gefundene Allel der falschen Elternart zuzuordnen, wenn man die minimalen eindeutigen Lesetiefen und Allelfrequenzen des Referenzpanels berücksichtigt. Die eindeutige Lesetiefe wurde als die Anzahl der Reads definiert, die eine bestimmte Stelle abdecken, wobei alle Reads von einem einzigen Individuum stammen. Diese Wahrscheinlichkeit P wurde wie in Gleichung 1 berechnet:
Wobei f(ps) die Allelfrequenz in der Elternspezies und urd(psp) die speziesspezifische eindeutige Lesetiefe im Referenzpanel ist.
Die Allelfrequenz der Elternspezies, die die höchste Wahrscheinlichkeit f(ps-max) ergibt, kann wie in Gleichung 2 beschrieben werden:
Durch Einsetzen von Gleichung 2 in Gleichung 1, wurde die maximale Wahrscheinlichkeit P(max), das in MCE1356 gefundene Allel der falschen Elternart zuzuordnen, wie in Gleichung 3 berechnet:
Die Ergebnisse zeigten, dass bei einer einzigartigen Lesetiefe von zwei, drei und vier in jeder Elternart, die maximale Wahrscheinlichkeit, das in MCE1356 gefundene Allel der falschen Elternart zuzuordnen, 0.148, 0,105 bzw. 0,082. Diese Maximalwerte sollten nicht als Fehlerquoten interpretiert werden, da sie sich nur ergeben würden, wenn alle Standorte innerhalb der Elternart variabel wären und alle Standorte einen MAF von genau (1 – (urd/urd + 1)) hätten. Unter der Annahme, dass die MAF-Verteilungen bei Belugas und Narwalen ähnlich sind, würde eine fehlerhafte Zuordnung von Allelen, die in MCE1356 gefunden wurden, beide Arten gleichermaßen betreffen und daher einen minimalen Einfluss auf die Rückschlüsse auf die Hybridisierung haben. Ein Vorteil der Verwendung der eindeutigen Lesetiefe in den Gleichungen 2 und 3 bestand darin, dass die Anzahl der Individuen und die Anzahl der Lesungen bei der Schätzung der Wahrscheinlichkeit, das in MCE1356 gefundene Allel der falschen Elternart zuzuordnen, kombiniert wurden. Anschließend wurden die Verteilung der Lesetiefe, der MAF (mit festem Haupt- und Nebenallel) und die Anzahl der Individuen mit Daten an jedem Standort für jede Elternart getrennt berechnet.
Es wurden drei Datensätze erstellt, die neben MCE1356 Elternarten-Panels mit einer minimalen Lesetiefe von zwei, drei und vier enthielten. Die Anzahl der in den drei Datensätzen enthaltenen Stellen betrug 61 105, 11 764 bzw. 360. Als Kompromiss zwischen der Maximierung der Anzahl der Sites und der Minimierung der falschen Zuordnung der in MCE1356 gefundenen Allele entschieden wir uns für den Datensatz mit einem Read in MCE1356 und einer minimalen eindeutigen Lesetiefe von drei in jeder Elternart. Dieser Datensatz umfasste 11.764 Stellen, die in den nachfolgenden Analysen verwendet wurden.
Zusammenfassende Statistiken wurden für den Datensatz mit 17 Individuen durchgeführt, einschließlich der Anzahl der Stellen, die (i) für verschiedene Allele in den Beluga- und Narwal-Arten-Panels fixiert waren; (ii) polymorph bei Belugas, aber nicht bei Narwalen; (iii) polymorph bei Narwalen, aber nicht bei Belugas; (iv) polymorph sowohl bei Belugas als auch bei Narwalen. Die Stellen, bei denen davon ausgegangen wird, dass sie zwischen den beiden Elternarten fixiert sind, werden mit Markern angereichert, die zwischen den beiden Elternarten stark differenziert sind, d. h. große Unterschiede in der Allelfrequenz aufweisen. Diese Marker sind zwar nicht notwendigerweise fixe Unterschiede zwischen den beiden Elternarten, aber dennoch sehr informativ für die Abstammung in MCE1356.
Wir haben die Genotyp-Likelihoods aus dem Datensatz ohne MCE1356 verwendet, um zu überprüfen, ob die Belugas und Narwale im Referenzpanel nicht selbst kürzlich gemischte Individuen waren. Wir schätzten ihre individuellen Vermischungskoeffizienten unter Angabe von zwei Populationen (K = 2) mit NgsAdmix28. Es wurden einhundert Durchläufe durchgeführt, und Mittelwert und Standardabweichung der Vermischungskoeffizienten wurden für die anschließende Auswertung verwendet. Um zu bestätigen, dass die Filter keine zuvor verborgenen genetischen Mischungsprofile im Referenzpanel aufgedeckt hatten, wurde diese Analyse sowohl vor als auch nach Anwendung der Filter für die eindeutige Lesetiefe durchgeführt.
Wir analysierten die Vermischungsanteile von MCE1356 mithilfe von fastNGSadmix29 unter Anwendung von 1.000 Bootstraps. fastNGSadmix verwendet Allelfrequenzen von Referenzpopulationen/-arten und die Genotypwahrscheinlichkeiten eines einzelnen Individuums, um dessen Vermischungsanteile zu schätzen. Die Software hat sich bei NGS-Daten mit einer Abdeckung von nur 0,00015×29 als robust erwiesen und war daher ideal für unsere Studie.
Wir schätzten den Hybridstatus von MCE1356 weiter ab, indem wir Standorte untersuchten, die für verschiedene Allele im Beluga-Panel und im Narwal-Panel (9.178 Standorte) fixiert waren, und den beobachteten Anteil von Reads, die mit dem belugaspezifischen Allel und dem narwalspezifischen Allel in MCE1356 übereinstimmten, mit den erwarteten Anteilen unter verschiedenen Hybridisierungsszenarien verglichen. Um festzustellen, wie gut die sieben verschiedenen Hybridisierungsszenarien mit den beobachteten Daten übereinstimmen, haben wir eine Pearson-Chi-Quadrat-Statistik für die Anpassungsgüte berechnet. Die Teststatistik wird wie in Gleichung 4 berechnet:
wobei Oi und Ei die beobachtete bzw. erwartete Anzahl von Allelen sind, die von der Elternart i stammen. Unter der Nullhypothese, dass das gewählte Szenario gut mit den beobachteten Daten übereinstimmt, folgt die Teststatistik T einer zentralen χ2-Verteilung. Somit würde das Szenario, das am besten mit den beobachteten Daten übereinstimmt, zur niedrigsten Teststatistik führen.
Um die sieben verschiedenen Hybridisierungsszenarien weiter zu untersuchen, berechneten wir die Wahrscheinlichkeit der beobachteten Allele an Stellen, die für verschiedene Allele in den Elternpopulationen fixiert waren. Insbesondere berechneten wir die Wahrscheinlichkeit der beobachteten Allele in der Hybride MCE1356 unter der Annahme eines Binomialmodells für die Vererbung der Allele von den Elternarten, wobei wir weiterhin von der Unabhängigkeit zwischen den Markern ausgingen. Wir möchten jedoch anmerken, dass die Verletzung der Unabhängigkeitsannahme immer noch zu unverzerrten Schätzungen führt. Unter der Annahme, dass der Hybrid MCE1356 aus einem Anteil b der Beluga-Abstammung und (1-b) der Narwal-Abstammung besteht, können wir die Wahrscheinlichkeit von b als Produkt der Wahrscheinlichkeit an k unabhängigen Stellen schreiben, angesichts der Anzahl der Reads nib, die mit dem Beluga-Allel an Stelle i übereinstimmen, und nin, die mit der Narwal-Abstammung übereinstimmen, wie in Gleichung 5.
Wir berechneten die Log-Wahrscheinlichkeit von b über seinen gesamten Bereich, von 0 bis 1, und verglichen die Log-Wahrscheinlichkeiten über die sieben verschiedenen Hybridisierungsszenarien.
Geschlechtsbestimmung
Um das Geschlecht von MCE1356 und den Individuen in den Beluga- und Narwal-Referenzpanels zu bestimmen, untersuchten wir das Verhältnis zwischen X-Chromosom und autosomaler Abdeckung. Dazu verglichen wir die Abdeckung von Gerüsten, die vermutlich vom X-Chromosom stammen, mit der von Gerüsten, die vermutlich von den Autosomen stammen. Wir bestimmten, welche Gerüste vom X-Chromosom stammen, indem wir das Genom des Killerwals mit SatsumaSynteny30 an das X-Chromosom der Kuh (Bos taurus) (CM008168.2) anglichen. Um Verzerrungen zu beseitigen, die aufgrund der Homologie zwischen einigen Regionen des X- und des Y-Chromosoms auftreten können, haben wir die resultierenden mutmaßlichen X-Chromosom-Gerüste weiter an das menschliche Y-Chromosom (NC_000024.10) angeglichen und diese Gerüste aus der weiteren Analyse entfernt. Anschließend berechneten wir die Abdeckung der verbleibenden X-Gerüste und der vier größten Gerüste, die nicht am X-Chromosom ausgerichtet waren (d. h. autosomale Gerüste), mithilfe der Tiefenfunktion in SAMtools17. Um Schwankungen in der Abdeckung im gesamten Genom auszugleichen, wurden 10.000.000 Stellen zufällig ausgewählt, die durchschnittliche Abdeckung für diese Stellen berechnet und dieser Schritt 100 Mal wiederholt. Für jedes Individuum wurden 5 %- und 95 %-Konfidenzintervalle sowie erste und dritte Quartile berechnet und für die anschließende Auswertung verwendet.
Stabile Kohlenstoff- und Stickstoffisotopenanalyse
Knochenpulverproben (~100 mg) von MCE1356, 18 Beluga- und 18 Narwalschädeln wurden mit 10 ml 2:1 Chloroform/Methanol (v/v) unter Beschallung 1 Stunde lang behandelt, um Lipide zu entfernen. Nach dem Entfernen des Lösungsmittels wurden die Proben 24 Stunden lang unter normalem atmosphärischem Druck getrocknet. Anschließend wurden die Proben 4 Stunden lang in 10 ml 0,5 M HCl entmineralisiert, während sie auf einem Orbitalschüttler geschüttelt wurden. Nach der Entmineralisierung wurden die Proben bis zur Neutralität mit Typ-I-Wasser gespült und anschließend 36 Stunden lang bei 75 °C in 10-3 M HCl erhitzt, um das Kollagen zu lösen. Das wasserlösliche Kollagen wurde dann gefriergetrocknet. Der Kollagenextrakt wurde dann mit 10 ml 10:5:4 Chloroform/Methanol/Wasser (v/v/v) unter Beschallung 1 Stunde lang behandelt, um alle restlichen Lipide zu entfernen. Nach dem Zentrifugieren wurde die Chloroform/Methanol-Schicht entfernt, und das Methanol wurde 24 Stunden lang bei 60 °C aus der Wasserschicht verdampft. Die stabile Isotopen- und Elementzusammensetzung von Kohlenstoff und Stickstoff wurde mit einem Nu Horizon (Nu Instruments, UK) Massenspektrometer mit kontinuierlichem Isotopenverhältnis an der Universität Trent, Kanada, bestimmt. Die stabilen Kohlenstoff- und Stickstoff-Isotopenzusammensetzungen wurden mit Hilfe der USGS40 und USGS41a relativ zu den VPDB- und AIR-Skalen kalibriert. Die Standardunsicherheit wurde mit ±0,17‰ für δ13C und ±0,22‰ für δ15N31 bestimmt; zusätzliche analytische Details sind im ergänzenden Dokument S4 enthalten. Die statistische Signifikanz der Unterschiede zwischen den Isotopenzusammensetzungen von Beluga und Narwal wurde mit ungepaarten t-Tests ermittelt.