Pobieranie próbek
Przeanalizowaliśmy (i) dane genomowego DNA wyodrębnionego z zębów czaszki (MCE1356, Fig. 2b, Supplementary Fig. S1) oraz z próbek tkanek ośmiu belugów i ośmiu narwali pobranych w Zatoce Disko, Zachodnia Grenlandia, oraz (ii) stabilne składy izotopowe węgla i azotu kolagenu kostnego z MCE1356 i panelu referencyjnego 18 czaszek belugów i 18 czaszek narwali z Zachodniej Grenlandii (Fig. 1c). Próbki tkanek (skóra) były przechowywane w 96% etanolu. Próbki zostały pobrane przez naukowców z Grenlandzkiego Instytutu Zasobów Naturalnych na podstawie ogólnego zezwolenia na pobieranie próbek biologicznych od Eskimosów wydanego przez rząd Grenlandii i wywiezione do Danii na podstawie zezwolenia CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 i 116.2006. Informacje o próbkach są wyszczególnione w Supplementary Table S2.
Ekstrakcja DNA, przygotowanie biblioteki i sekwencjonowanie
Próbki tkanek
DNA z próbek tkanek ekstrahowano przy użyciu Qiagen Blood and Tissue Kit zgodnie z protokołem producenta. DNA zostało pofragmentowane przy użyciu Covaris M220 Focused-ultrasonicator w celu utworzenia fragmentów o długości ~350-550 par zasad (bp). Biblioteki budowano z fragmentowanych ekstraktów DNA przy użyciu programu Illumina NeoPrep zgodnie z przewodnikiem NeoPrep Library Prep System Guide z zastosowaniem ustawień domyślnych. Amplifikacja PCR, kwantyfikacja i normalizacja zostały przeprowadzone przez system przygotowania bibliotek NeoPrep. Biblioteki przesiewano pod kątem rozkładu wielkości na bioanalizatorze Agilent 2100 i łączono w równych proporcjach przed sekwencjonowaniem na urządzeniu Illumina HiSeq 2500 z technologią 80 bp SE.
Próbki zębów/kości
W przeciwieństwie do próbek tkanek, próbki starej kości i zębów mają stosunkowo niskie stężenie DNA, co wymaga innych protokołów ekstrakcji i tworzenia bibliotek. Około 0,5 g proszku kostnego z pięciu zębów i jeden odłamek kości z okazu MCE1356 został wywiercony przy użyciu ręcznego narzędzia Dremel. DNA zostało wyekstrahowane z proszku zębowego/kostnego w specjalnym laboratorium czyszczenia starożytnego DNA w Natural History Museum of Denmark, University of Copenhagen, przy użyciu buforu ekstrakcyjnego opisanego w11 z dodatkiem 30-minutowego etapu wstępnego trawienia12. Zamiast używać kolumn Zymo-Spin V (Zymo Research), bufor ekstrakcyjny został zagęszczony przy użyciu jednostek filtracyjnych Amicon Ultra 30 kDa, a następnie dalej zagęszczony i oczyszczony przy użyciu probówek Qiagen Minelute. Oczyszczone ekstrakty zostały następnie wbudowane w biblioteki Illumina zgodnie z protokołem opisanym przez13. Użyliśmy qPCR do sprawdzenia, czy budowanie biblioteki zakończyło się sukcesem, do wybrania bibliotek do sekwencjonowania i do obliczenia odpowiedniej liczby cykli PCR wymaganych do wystarczającej amplifikacji każdej biblioteki bez powodowania nadmiernej amplifikacji. W sumie, cztery biblioteki zostały amplifikowane z unikalnymi indeksami 6 bp i przesiewane pod kątem zawartości endogennej na platformie Illumina MiSeq przy użyciu sekwencjonowania 250 bp SE. Najlepsze biblioteki zostały następnie poddane ponownej sekwencjonowaniu na platformie Illumina HiSeq 2500 przy użyciu technologii 80 bp SE.
Analiza bioinformatyczna
Wszystkie analizy mapowania i uszkodzeń DNA zostały wykonane w ramach potoku Paleomix 1.2.1214. Odczyty zostały przycięte za pomocą AdapterRemoval 2.2.015 przy użyciu ustawień domyślnych, z wyjątkiem minimalnej długości odczytu, która została ustawiona na 25 bp. Odczyty były kontrolowane przy użyciu FastQC i wyrównywane przy użyciu BWA16 z zastosowaniem algorytmu Backtrack i wyłączeniem początkowej długości ziarna. Jeśli odczyty mapowały do wielu pozycji lub miały wynik jakości mapowania (MAPQ score z BWA) mniejszy niż 30, usuwano je za pomocą SAMtools17. Duplikaty sekwencji zostały usunięte przy użyciu MarkDuplicates z Picarda (dostępne pod adresem: http://broadinstitute.github.io/picard), a ostateczne wyrównanie wokół indeli wykonano przy użyciu GATK18. Deaminacja cytozyny do uracylu w próbce MCE1356 została oceniona przy użyciu danych wyjściowych z programu mapDamage v2.0.619. Wyniki MapDamage nie wykazały żadnego wyraźnego sygnału deaminacji w sekwencjach (Supplementary Fig. S3).
Analiza mitochondrialna
Aby określić linię matczyną MCE1356, odczyty sekwencjonowania DNA zmapowano do mitochondrialnych genomów referencyjnych belugi (GenBank accession: KY444734) i narwala (GenBank accession: NC_005279) i porównano średnie pokrycie. Odczyty z ośmiu próbek belugi i ośmiu próbek narwala, które stanowiły panel referencyjny, zostały zmapowane do ich odpowiednich mitochondrialnych genomów referencyjnych. Skonstruowaliśmy mitochondrialne sekwencje konsensusowe dla MCE1356 i 16 próbek panelu referencyjnego z regionami pokrytymi przez więcej niż pięć odczytów przy użyciu BEDtools20, SAMtools17 i GATK18. Utworzyliśmy dwa wyrównania sekwencji przy użyciu ClustalW21, stosując ustawienia domyślne, które obejmowały 16 próbek panelu referencyjnego i albo wersję MCE1356 zmapowaną do mitochondrialnego genomu referencyjnego belugi, albo wersję MCE1356 zmapowaną do mitochondrialnego genomu referencyjnego narwala. Użyliśmy tych dwóch wyrównań do skonstruowania medianowych sieci haplotypów22 przy użyciu Popart 1.723 (dostępny w: http://popart.otago.ac.nz), wykluczając wszelkie miejsca z indelami lub brakującymi danymi. Następnie, zarówno pokrycie po mapowaniu, jak i dwie sieci haplotypów zostały użyte do określenia gatunku linii matczynej MCE1356.
Analizy jądrowego DNA
Zmapowaliśmy odczyty sekwencjonowania DNA ze wszystkich próbek do genomu referencyjnego wieloryba zabójcy (Orcinus orca) (GCA_000331955.2). Wysokiej jakości genom wieloryba beluga został niedawno opublikowany24, ale mapowanie do jednego z dwóch potencjalnych gatunków rodzicielskich mogłoby spowodować błędy w analizach. Dlatego zmapowaliśmy odczyty do genomu wieloryba zabójcy, ponieważ jest on dobrze złożony, a wieloryby zabójcy są nadal stosunkowo blisko spokrewnione z belugami i narwalami (czas dywergencji 12 MYA)25, ale wystarczająco odległe, aby zminimalizować ryzyko introgresji, które skomplikowałoby nasze analizy.
Do wszystkich dalszych filtrów użyliśmy ANGSD v0.92326, pakietu oprogramowania, który wykorzystuje prawdopodobieństwa genotypów zamiast nazwanych genotypów, co jest szczególnie przydatne podczas analizy danych NGS o niskim pokryciu. Użyliśmy metody SAMtools17 zaimplementowanej w ANGSD do oszacowania prawdopodobieństw genotypów i wywnioskowaliśmy główne i pomniejsze allele bezpośrednio z prawdopodobieństw genotypów przy użyciu podejścia maksymalnej wiarygodności, jak opisano w27.
Aby zwizualizować zmapowane dane, wykreśliliśmy rozkład głębokości odczytu ze wszystkich osobników, wyłączając miejsca z więcej niż dwoma allelami i miejscami z wynikiem Phred poniżej 25. Wizualizacja danych ze wszystkich osobników łącznie umożliwiła nam oszacowanie średniej głębokości odczytu 4,14x i zidentyfikowanie miejsc o podwyższonej głębokości odczytu. Miejsca takie pochodziły najprawdopodobniej z paralogów i powtarzających się regionów genomu. Zestaw danych został poddany wizualnej inspekcji i dalszej filtracji, wykluczając wszystkie miejsca o głębokości odczytu większej niż dziewięć (6,9% miejsc). W ANGSD, SNP były nazywane na podstawie częstości alleli. Częstość alleli mniejszościowych (MAF) została oszacowana na podstawie prawdopodobieństwa genotypu, a test współczynnika prawdopodobieństwa został użyty do określenia, czy MAF jest różna od zera. Jeśli wartość p z testu stosunku prawdopodobieństwa wynosiła <1e-4, miejsce uznawano za polimorficzne i zachowywano w zbiorze danych. Zastosowanie tego filtra SNP oznaczało, że żadne miejsca z mniej niż czterema odczytami nie zostały zachowane, ponieważ miejsca objęte mniejszą liczbą odczytów nie mogły uzyskać tak niskich wartości p. Filtry te zastosowano do zbioru danych obejmującego wszystkie 17 próbek oraz do zbioru danych bez MCE1356.
Aby ustalić, czy MCE1356 był hybrydą beluga/narwal, dokonaliśmy dalszej filtracji zbioru danych 17 osobników, wykluczając miejsca bez odczytów w MCE1356. W tym momencie średnia głębokość odczytu MCE1356 była tylko 1,15x, więc aby upewnić się, że nie analizujemy loci wielokopii, wykluczyliśmy miejsca objęte więcej niż jednym odczytem w MCE1356.
Naszym celem było porównanie alleli znalezionych w MCE1356 z allelami w panelu referencyjnym, więc oszacowaliśmy prawdopodobieństwo przypisania allelu znalezionego w MCE1356 do niewłaściwego gatunku rodzicielskiego, biorąc pod uwagę różne minimalne unikalne głębokości odczytu panelu referencyjnego i częstotliwości alleli. Unikalna głębokość odczytu została zdefiniowana jako liczba odczytów obejmujących określone miejsce, gdzie wszystkie odczyty pochodzą od unikalnego osobnika. Prawdopodobieństwo P zostało obliczone jak w równaniu 1:
Gdzie f(ps) to częstotliwość alleli w gatunku rodzicielskim, a urd(psp) to specyficzna dla gatunku unikalna głębokość odczytu w panelu referencyjnym.
Częstotliwość alleli gatunku rodzicielskiego dająca najwyższe prawdopodobieństwo f(ps-max) może być opisana jak w równaniu 2:
Wstawiając równanie 2 do równania 1, obliczono maksymalne prawdopodobieństwo P(max) przypisania allelu znalezionego w MCE1356 do niewłaściwego gatunku rodzicielskiego, jak w równaniu 3:
Results revealed that with a unique read depth of two, trzy i cztery w każdym gatunku rodzicielskim, maksymalne prawdopodobieństwo przypisania allelu znalezionego w MCE1356 do niewłaściwego gatunku rodzicielskiego wynosiło 0.148, 0,105 i 0,082, odpowiednio. Te maksymalne wartości nie powinny być interpretowane jako wskaźniki błędu, ponieważ zostałyby one uzyskane tylko wtedy, gdyby wszystkie stanowiska były zmienne w obrębie gatunków rodzicielskich, a wszystkie stanowiska miały MAF równy dokładnie (1 – (urd/urd + 1)). Ponadto, zakładając, że rozkłady MAF u belugów i narwali były podobne, błędne przypisanie alleli znalezionych w MCE1356 dotknęłoby oba gatunki w równym stopniu, a zatem miałoby minimalny wpływ na wnioskowanie o hybrydyzacji. Zaletą zastosowania unikalnej głębokości odczytu w równaniach 2 i 3 było połączenie liczby osobników i liczby odczytów w oszacowaniu prawdopodobieństwa przypisania allelu znalezionego w MCE1356 do niewłaściwego gatunku rodzicielskiego. Następnie rozkład głębokości odczytu, MAF (z ustalonym głównym i mniejszym allelem) oraz liczbę osobników z danymi w każdym miejscu obliczono oddzielnie dla każdego gatunku rodzicielskiego.
Zbudowano trzy zestawy danych, które oprócz MCE1356 obejmowały panele gatunków rodzicielskich z minimalną unikalną głębokością odczytu dwóch, trzech i czterech. Liczba miejsc zachowanych w tych trzech zestawach danych wynosiła odpowiednio 61,105, 11,764 i 360. Jako kompromis pomiędzy maksymalizacją liczby miejsc a minimalizacją błędnego przypisania alleli znalezionych w MCE1356, zdecydowaliśmy się użyć zestawu danych z jednym odczytem w MCE1356 i minimalną unikalną głębokością odczytu wynoszącą trzy w każdym gatunku rodzicielskim. Ten zbiór danych zawierał 11 764 miejsc, które zostały wykorzystane w kolejnych analizach.
Podsumowujące statystyki zostały wykonane na zbiorze danych z 17 osobnikami, w tym liczba miejsc, które były (i) ustalone dla różnych alleli w panelach gatunków beluga i narwala; (ii) polimorficzne w belugach, ale nie w narwalach; (iii) polimorficzne w narwalach, ale nie w belugach; (iv) polimorficzne zarówno w belugach, jak i narwalach. Miejsca, które szacuje się jako stałe między dwoma panelami gatunków rodzicielskich, będą wzbogacone o markery, które są wysoce zróżnicowane, tj. mają duże różnice w częstotliwości alleli między dwoma gatunkami rodzicielskimi. Te markery, chociaż niekoniecznie stałe różnice między dwoma gatunkami rodzicielskimi, są nadal wysoce informatywne dla rodowodu w MCE1356.
Użyliśmy prawdopodobieństw genotypów ze zbioru danych bez MCE1356, aby sprawdzić, czy belugi i narwale w panelu referencyjnym nie były same niedawno zmieszane osobniki. Oszacowaliśmy ich indywidualne współczynniki domieszkowania, określając dwie populacje (K = 2) przy użyciu NgsAdmix28. Wykonano sto przebiegów, a średnia i odchylenie standardowe współczynników domieszkowania posłużyły do późniejszej interpretacji. Aby potwierdzić, że filtry nie ujawniły wcześniej ukrytych profili genetycznych domieszek w panelu referencyjnym, analiza ta została przeprowadzona zarówno przed, jak i po zastosowaniu filtrów unikalnej głębokości odczytu.
Zanalizowaliśmy proporcje domieszek w MCE1356 przy użyciu fastNGSadmix29, stosując 1000 bootstrapów. fastNGSadmix wykorzystuje częstotliwości alleli populacji/gatunków referencyjnych i prawdopodobieństwa genotypu pojedynczego osobnika, aby oszacować jego proporcje domieszek. Oprogramowanie to okazało się odporne na dane NGS o pokryciu tak niskim jak 0,00015×29, a zatem było idealne do naszego badania.
Oszacowaliśmy dalej status hybrydowy MCE1356 badając miejsca ustalone dla różnych alleli w panelu belugi i panelu narwala (9 178 miejsc) i porównując zaobserwowany odsetek odczytów pasujących do alleli specyficznych dla belugi i alleli specyficznych dla narwala w MCE1356 do oczekiwanych proporcji w różnych scenariuszach hybrydyzacji. Aby określić, jak dobrze siedem różnych scenariuszy hybrydyzacji pasowało do obserwowanych danych, obliczyliśmy statystykę dobroci dopasowania Pearsona Chi-square. Statystyka testu jest obliczana jak w równaniu 4:
gdzie Oi i Ei są odpowiednio obserwowaną i oczekiwaną liczbą alleli pochodzących od gatunku rodzicielskiego i. Przy hipotezie zerowej, w której wybrany scenariusz dobrze odpowiada zaobserwowanym danym, statystyka testowa T ma rozkład centralny χ2. Tak więc scenariusz, który najlepiej odpowiada obserwowanym danym, prowadziłby do najniższej statystyki testowej.
Aby dokładniej zbadać siedem różnych scenariuszy hybrydyzacji, obliczyliśmy prawdopodobieństwo obserwowanych alleli w miejscach, które były ustalone dla różnych alleli w populacjach rodzicielskich. Konkretnie, obliczyliśmy prawdopodobieństwo obserwowanych alleli w hybrydzie MCE1356, w ramach modelu dwumianowego dla dziedziczenia alleli z gatunków rodzicielskich, dalej zakładając niezależność między markerami. Chcielibyśmy jednak zauważyć, że naruszenie założenia o niezależności nadal prowadzi do nieobiektywnych oszacowań. Zakładając, że hybryda MCE1356 składa się z proporcji b przodków belugi i (1-b) przodków narwala, możemy zapisać prawdopodobieństwo b jako iloczyn prawdopodobieństwa w k niezależnych miejscach, biorąc pod uwagę liczbę odczytów nib, które pasują do allelu belugi w miejscu i i nin, które pasują do przodków narwala, jak w równaniu 5.
Obliczyliśmy log-likelihood dla b w całym jego zakresie, od 0 do 1, i porównaliśmy log-likelihoody w siedmiu różnych scenariuszach hybrydyzacji.
Określanie płci
Aby określić płeć MCE1356 oraz osobników z paneli referencyjnych belugi i narwala, zbadaliśmy stosunek pokrycia chromosomu X do autosomalnego. Dokonano tego poprzez porównanie pokrycia rusztowań pochodzących przypuszczalnie z chromosomu X, z pokryciem rusztowań pochodzących przypuszczalnie z autosomów. Określiliśmy, które rusztowania pochodzą z chromosomu X poprzez wyrównanie genomu wieloryba zabójcy do chromosomu X krowy (Bos taurus) (CM008168.2) za pomocą SatsumaSynteny30. Ponadto, aby usunąć błędy, które mogą wystąpić z powodu homologii między niektórymi regionami chromosomów X i Y, wyrównaliśmy otrzymane rusztowania chromosomu X do chromosomu Y człowieka (NC_000024.10) i usunęliśmy te rusztowania z dalszej analizy. Następnie obliczyliśmy pokrycie dla pozostałych rusztowań chromosomu X oraz czterech największych rusztowań nie wyrównujących się do chromosomu X (tj. rusztowań autosomalnych), używając funkcji depth w SAMtools17. Aby skompensować wszelkie różnice w pokryciu w całym genomie, wybraliśmy losowo 10 000 000 miejsc, obliczyliśmy średnie pokrycie dla tych miejsc i powtórzyliśmy ten krok 100 razy. Dla każdego osobnika obliczono 5% i 95% przedziały ufności, a także pierwszy i trzeci kwartyl, które wykorzystano do późniejszej interpretacji.
Analiza stabilnego węgla i izotopu azotu
Sproszkowane próbki kości (~100 mg) z czaszek MCE1356, 18 beluga i 18 narwala poddano działaniu 10 ml 2:1 chloroformu/metanolu (v/v) pod sonikacją przez 1 h w celu usunięcia lipidów. Po usunięciu rozpuszczalnika próbki suszono pod normalnym ciśnieniem atmosferycznym przez 24 h. Następnie próbki poddano demineralizacji w 10 ml 0,5 M HCl przez 4 h, mieszając je na wytrząsarce orbitalnej. Po demineralizacji, próbki płukano do neutralności wodą typu I, a następnie ogrzewano w temperaturze 75 °C przez 36 h w 10-3 M HCl w celu rozpuszczenia kolagenu. Rozpuszczalny w wodzie kolagen był następnie liofilizowany. Ekstrakt kolagenu był następnie traktowany 10 ml 10:5:4 chloroformu/metanolu/wody (v/v/v) pod sonikacją przez 1 h w celu usunięcia wszelkich pozostałości lipidów. Po odwirowaniu usunięto warstwę chloroformu/metanolu, a metanol odparowano z warstwy wodnej w temperaturze 60 °C przez 24 h. Próbki ponownie liofilizowano i odważono w cynowych kapsułkach do analizy pierwiastkowej i izotopowej. Stabilny izotop węgla i azotu oraz skład pierwiastkowy oznaczono przy użyciu spektrometru masowego Nu Horizon (Nu Instruments, UK) o ciągłym przepływie izotopów na Uniwersytecie Trent w Kanadzie. Stabilne izotopowe składy węgla i azotu zostały skalibrowane względem skali VPDB i AIR przy użyciu USGS40 i USGS41a. Niepewność standardową określono na ±0,17‰ dla δ13C i ±0,22‰ dla δ15N31; dodatkowe szczegóły analityczne podano w Dokumencie Uzupełniającym S4. Statystyczną istotność różnic między składem izotopowym belugi i narwala oceniano za pomocą nieparowanych testów t.
.