Výběr vzorků
Analyzovali jsme (i) celogenomová data DNA získaná ze zubů lebky (MCE1356, obr. 2b, doplňkový obr. S1) a ze vzorků tkání osmi běluh a osmi narvalů odebraných v zátoce Disko v západním Grónsku a (ii) stabilní izotopové složení uhlíku a dusíku kostního kolagenu z MCE1356 a referenčního panelu 18 lebek běluh a 18 lebek narvalů ze západního Grónska (obr. 1c). Vzorky tkání (kůže) byly uchovávány v 96% ethanolu. Vzorky byly odebrány vědci z Grónského institutu přírodních zdrojů na základě všeobecného povolení grónské vlády pro odběr biologických vzorků od Inuitů a vyvezeny do Dánska na základě povolení CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 a 116.2006. Informace o vzorcích jsou podrobně uvedeny v doplňkové tabulce S2.
Extrakce DNA, příprava knihovny a sekvenování
Vzorky tkání
DNA ze vzorků tkání byla extrahována pomocí soupravy Qiagen Blood and Tissue Kit podle protokolu výrobce. DNA byla fragmentována pomocí přístroje Covaris M220 Focused-ultrasonicator, aby se vytvořily fragmenty o délce ~350-550 párů bází (bp). Knihovny byly sestaveny z fragmentovaných extraktů DNA pomocí programu Illumina NeoPrep podle příručky systému přípravy knihoven NeoPrep s použitím výchozího nastavení. Amplifikace PCR, kvantifikace a normalizace byly provedeny pomocí systému NeoPrep Library Prep System. Knihovny byly zkontrolovány z hlediska distribuce velikosti na přístroji Agilent 2100 Bioanalyzer a před sekvenováním na přístroji Illumina HiSeq 2500 s technologií 80 bp SE byly spojeny v ekvimolárních poměrech.
Vzorky zubů/kostí
Na rozdíl od vzorků tkání mají vzorky starých kostí a zubů relativně nízkou koncentraci DNA, což vyžaduje odlišné protokoly extrakce a sestavení knihovny. Přibližně 0,5 g kostního prášku z pěti zubů a jednoho kostního úlomku ze vzorku MCE1356 bylo vyvrtáno pomocí ručního dremelu. DNA byla ze zubního/kostního prášku extrahována ve specializované laboratoři pro čištění starověké DNA v Dánském přírodovědném muzeu na Kodaňské univerzitě za použití extrakčního pufru popsaného v článku11 s přidáním 30minutové předtrávicí fáze12. Místo použití kolon Zymo-Spin V (Zymo Research) byl extrakční pufr zahuštěn pomocí odstředivých filtračních jednotek Amicon Ultra 30 kDa a dále zahuštěn a vyčištěn pomocí zkumavek Qiagen Minelute. Přečištěné extrakty byly poté sestaveny do knihoven Illumina podle protokolu popsaného13. Pomocí qPCR jsme zkontrolovali, zda sestavení knihovny proběhlo úspěšně, vybrali knihovny k sekvenování a vypočítali vhodný počet cyklů PCR potřebných k dostatečné amplifikaci každé knihovny, aniž by došlo k nadměrné amplifikaci. Celkem byly amplifikovány čtyři knihovny s unikátními 6bp indexy, které byly prověřeny na endogenní obsah na platformě Illumina MiSeq pomocí 250bp SE sekvenování. Nejlepší knihovny byly poté resekvenovány na platformě Illumina HiSeq 2500 pomocí technologie 80 bp SE.
Bioinformatická analýza
Všechny analýzy mapování a poškození DNA byly provedeny v rámci pipeline Paleomix 1.2.1214. Čtení byla ořezána pomocí AdapterRemoval 2.2.015 s použitím výchozího nastavení s výjimkou minimální délky čtení, která byla nastavena na 25 bp. Čtení byla zkontrolována pomocí FastQC a zarovnána pomocí BWA16 s použitím algoritmu Backtrack a vypnutím počáteční délky osiva. Pokud se čtení mapovala na více pozic nebo měla skóre kvality mapování (MAPQ skóre z BWA) nižší než 30, byla odstraněna pomocí SAMtools17. Duplikáty sekvencí byly odstraněny pomocí MarkDuplicates z programu Picard (dostupné z: http://broadinstitute.github.io/picard) a konečné zarovnání bylo znovu zarovnáno kolem indelů pomocí GATK18. Deaminace cytosinu na uracil u vzorku MCE1356 byla posouzena pomocí výstupu z mapDamage v2.0.619. Výsledky mapDamage neprokázaly v sekvencích žádný jasný signál deaminace (doplňkový obr. S3).
Mitochondriální analýza
Pro určení mateřské linie MCE1356 byla čtení sekvenování DNA mapována na referenční mitochondriální genomy belugy (GenBank accession: KY444734) a narvala (GenBank accession: NC_005279) a bylo porovnáno průměrné pokrytí. Čtení z osmi vzorků belugy a osmi vzorků narvalů, které tvořily referenční panel, byla mapována na jejich příslušné mitochondriální referenční genomy. Pomocí nástrojů BEDtools20, SAMtools17 a GATK18 jsme sestavili mitochondriální konsenzuální sekvence MCE1356 a 16 vzorků referenčního panelu s oblastmi pokrytými více než pěti čteními. Pomocí programu ClustalW21 jsme vytvořili dvě zarovnání sekvencí s použitím výchozího nastavení, která zahrnovala 16 referenčních panelových vzorků a buď verzi MCE1356 mapovanou na referenční mitochondriální genom belugy, nebo verzi MCE1356 mapovanou na referenční mitochondriální genom narvala. Obě zarovnání jsme použili ke konstrukci mediánových spojovacích haplotypových sítí22 pomocí programu Popart 1.723 (dostupný z: http://popart.otago.ac.nz), přičemž jsme vyloučili všechna místa s indely nebo chybějícími údaji. Následně jsme k určení druhu mateřské linie MCE1356 použili jak pokrytí po mapování, tak obě haplotypové sítě.
Analýzy jaderné DNA
Mapovali jsme čtení sekvenování DNA ze všech vzorků na referenční genom kosatky dravé (Orcinus orca) (GCA_000331955.2). Nedávno byl zveřejněn vysoce kvalitní genom velryby běluhy24 , ale mapování na jeden ze dvou potenciálních rodičovských druhů by mohlo způsobit zkreslení analýz. Proto jsme čtení mapovali na genom kosatky, protože je dobře sestavený a kosatky jsou stále relativně blízce příbuzné s běluhami a narvaly (doba divergence 12 MYA)25, ale zároveň dostatečně vzdálené, aby se minimalizovalo riziko introgrese, která by komplikovala naše analýzy.
Pro veškeré další filtrování jsme použili ANGSD v0.92326, softwarový balík, který používá pravděpodobnosti genotypů místo volaných genotypů, což je užitečné zejména při analýze dat NGS s malým pokrytím. K odhadu genotypových pravděpodobností jsme použili metodu SAMtools17 implementovanou v ANGSD a hlavní a vedlejší alely jsme odvodili přímo z genotypových pravděpodobností pomocí přístupu maximální věrohodnosti, jak je popsáno v27.
Pro vizualizaci mapovaných dat jsme vykreslili rozložení hloubky čtení ze všech jedinců, přičemž jsme vyloučili místa s více než dvěma alelami a místa se skóre Phred nižším než 25 bodů. Vizualizace dat od všech jedinců dohromady nám umožnila odhadnout průměrnou hloubku čtení 4,14x a identifikovat místa se zvýšenou hloubkou čtení. Taková místa s větší pravděpodobností pocházela z paralogů a repetitivních oblastí genomu. Soubor dat byl vizuálně zkontrolován a dále filtrován, přičemž byla vyloučena všechna místa s hloubkou čtení větší než devět (6,9 % míst). V ANGSD byly SNP pojmenovány na základě jejich alelových frekvencí. Frekvence minoritních alel (MAF) byla odhadnuta z pravděpodobností genotypů a k určení, zda se MAF liší od nuly, byl použit test poměru pravděpodobnosti. Pokud byla p hodnota z testu poměru pravděpodobnosti <1e-4, bylo místo považováno za polymorfní a bylo ponecháno v souboru dat. Použití tohoto filtru SNP znamenalo, že nebyly ponechány žádné lokality s méně než čtyřmi čteními, protože lokality pokryté menším počtem čtení nemohly získat tak nízké hodnoty p. Tyto filtry byly použity na soubor dat zahrnující všech 17 vzorků a na soubor dat bez MCE1356.
Abychom zjistili, zda MCE1356 byl křížencem belugy a nosorožce, dále jsme filtrovali soubor dat 17 jedinců a vyloučili místa bez čtení v MCE1356. V tomto okamžiku byla průměrná hloubka čtení MCE1356 pouze 1,15x, takže abychom se ujistili, že neanalyzujeme vícekopiové lokusy, vyloučili jsme místa pokrytá více než jedním čtením v MCE1356.
Naším cílem bylo porovnat alely nalezené v MCE1356 s alelami v referenčním panelu, takže jsme odhadli pravděpodobnost přiřazení alely nalezené v MCE1356 nesprávnému rodičovskému druhu vzhledem k různým minimálním jedinečným hloubkám čtení a frekvencím alel v referenčním panelu. Unikátní hloubka čtení byla definována jako počet čtení pokrývajících konkrétní místo, kde všechna čtení pocházejí od unikátního jedince. Tato pravděpodobnost P byla vypočtena podle rovnice 1:
Kde f(ps) je frekvence alel u rodičovského druhu a urd(psp) je hloubka jedinečných čtení specifická pro daný druh v referenčním panelu.
Frekvenci alel rodičovského druhu, která dává nejvyšší pravděpodobnost f(ps-max), lze popsat podle rovnice 2:
Vložením rovnice 2 do rovnice 1, byla vypočtena maximální pravděpodobnost P(max) přiřazení alely nalezené v MCE1356 nesprávnému rodičovskému druhu podle rovnice 3:
Výsledky ukázaly, že při jedinečné hloubce čtení dva, tři a čtyři u každého rodičovského druhu byla maximální pravděpodobnost přiřazení alely nalezené v MCE1356 nesprávnému rodičovskému druhu 0.148, 0,105 a 0,082. Tyto maximální hodnoty by neměly být interpretovány jako chybovost, protože by jich bylo dosaženo pouze v případě, že by všechna místa byla v rámci rodičovských druhů variabilní a všechna místa by měla MAF přesně (1 – (urd/urd + 1)). Navíc za předpokladu, že by rozložení MAF u běluh a narvalů bylo podobné, chybné přiřazení alel nalezených v MCE1356 by se týkalo obou druhů stejně, a proto by mělo minimální vliv na závěry o hybridizaci. Výhodou použití unikátní hloubky čtení v rovnicích 2 a 3 bylo, že se při odhadu pravděpodobnosti přiřazení alely nalezené v MCE1356 nesprávnému rodičovskému druhu kombinoval počet jedinců a počet čtení. Následně bylo pro každý rodičovský druh zvlášť vypočteno rozložení hloubky čtení, MAF (s fixní hlavní a vedlejší alelou) a počet jedinců s daty v každém místě.
Byly zkonstruovány tři soubory dat, které kromě MCE1356 zahrnovaly panely rodičovských druhů s minimální unikátní hloubkou čtení dvě, tři a čtyři. Počet lokalit zachovaných v těchto třech souborech dat byl 61 105, 11 764 a 360. V těchto třech souborech dat byly zachyceny i další lokality. Jako kompromis mezi maximalizací počtu míst a minimalizací chybného přiřazení alel nalezených v MCE1356 jsme se rozhodli použít datovou sadu s jedním čtením v MCE1356 a minimální unikátní hloubkou čtení tři v každém rodičovském druhu. Tento soubor dat obsahoval 11 764 míst, která byla použita v následných analýzách.
Souhrnná statistika byla provedena na souboru dat se 17 jedinci, včetně počtu míst, která byla (i) fixována pro různé alely v panelech druhů běluh a narval; (ii) polymorfní u běluh, ale ne u narval; (iii) polymorfní u narval, ale ne u běluh; (iv) polymorfní u běluh i narval. Místa, která se odhadují jako fixní mezi oběma rodičovskými druhy, budou obohacena o markery, které jsou vysoce diferencované, tj. mají velké rozdíly ve frekvenci alel mezi oběma rodičovskými druhy. Tyto markery, ačkoli nemusí nutně znamenat fixní rozdíly mezi oběma rodičovskými druhy, jsou přesto vysoce informativní pro původ v MCE1356.
Použili jsme pravděpodobnosti genotypů ze souboru dat bez MCE1356, abychom ověřili, že běluhy a narvalové v referenčním panelu nebyli sami nedávno smíšenými jedinci. Odhadli jsme jejich individuální koeficienty příměsi při zadání dvou populací (K = 2) pomocí NgsAdmix28. Bylo provedeno sto běhů a pro následnou interpretaci byly použity průměr a směrodatná odchylka koeficientů příměsí. Abychom potvrdili, že filtry neodhalily dříve skryté příměsové genetické profily v referenčním panelu, byla tato analýza provedena jak před, tak i po použití filtrů jedinečné hloubky čtení.
Analyzovali jsme příměsové poměry MCE1356 pomocí fastNGSadmix29 s použitím 1000 bootstrapů. fastNGSadmix používá frekvence alel referenčních populací/druhů a pravděpodobnosti genotypů jednoho jedince k odhadu jeho příměsových poměrů. Tento software se ukázal jako robustní u dat NGS s pokrytím až 0,00015×29, a proto byl pro naši studii ideální.
Dále jsme odhadovali hybridní status MCE1356 zkoumáním míst fixovaných pro různé alely v panelu beluga a panelu narval (9 178 míst) a porovnáváním pozorovaného podílu čtení odpovídajících alele specifické pro belugu a alele specifické pro narvala v MCE1356 s očekávanými podíly při různých scénářích hybridizace. Abychom určili, jak dobře sedm různých hybridizačních scénářů odpovídá pozorovaným údajům, vypočítali jsme Pearsonův chí-kvadrát statistiku dobré shody. Testovací statistika se vypočítá podle rovnice 4:
kde Oi a Ei jsou pozorované a očekávané počty alel odvozených od rodičovského druhu i, v tomto pořadí. Za nulové hypotézy, kdy zvolený scénář dobře odpovídá pozorovaným údajům, se testovací statistika T řídí centrálním rozdělením χ2. Tedy scénář, který nejlépe odpovídá pozorovaným datům, by vedl k nejnižší testové statistice.
Pro další zkoumání sedmi různých scénářů hybridizace jsme vypočítali pravděpodobnost pozorovaných alel na místech, která byla fixována pro různé alely v rodičovských populacích. Konkrétně jsme vypočítali pravděpodobnost pozorovaných alel v hybridu MCE1356 podle binomického modelu pro dědičnost alel z rodičovských druhů, dále za předpokladu nezávislosti mezi markery. Rádi bychom však poznamenali, že porušení předpokladu nezávislosti stále vede k nezkresleným odhadům. Za předpokladu, že hybrid MCE1356 je složen z podílu b předků belugy a (1-b) předků narvala, můžeme pravděpodobnost b zapsat jako součin pravděpodobnosti na k nezávislých místech, vzhledem k počtu čtení nib, která odpovídají alele belugy v místě i, a nin, která odpovídají předkům narvala, jak je uvedeno v rovnici 5.
Vypočítali jsme logaritmickou pravděpodobnost b v celém jeho rozsahu, od 0 do 1 a porovnali logaritmické pravděpodobnosti v sedmi různých scénářích hybridizace.
Determinace pohlaví
Pro určení pohlaví MCE1356 a jedinců v referenčních panelech belugy a narvala jsme zkoumali poměr pokrytí chromozomu X a autozomu. To bylo provedeno porovnáním pokrytí napříč scaffoldy domněle pocházejícími z chromozomu X a scaffoldy domněle pocházejícími z autozomů. Určili jsme, které scaffoldy pocházejí z chromozomu X, zarovnáním genomu kosatky s chromozomem X krávy (Bos taurus) (CM008168.2) pomocí SatsumaSynteny30. Navíc, abychom odstranili zkreslení, ke kterému může dojít v důsledku homologie mezi některými oblastmi chromozomů X a Y, jsme výsledné domnělé scaffoldy chromozomu X dále zarovnali k lidskému chromozomu Y (NC_000024.10) a tyto scaffoldy jsme z další analýzy odstranili. Poté jsme vypočítali pokrytí zbývajících scaffoldů X a čtyř největších scaffoldů nezarovnaných k chromozomu X (tj. autozomálních scaffoldů) pomocí funkce depth v nástroji SAMtools17. Abychom kompenzovali případné rozdíly v pokrytí napříč genomem, náhodně jsme vybrali 10 000 000 míst, vypočítali průměrné pokrytí pro tato místa a tento krok 100krát zopakovali. Pro každého jedince byly vypočteny 5% a 95% intervaly spolehlivosti a také první a třetí kvartil, které byly použity pro následnou interpretaci.
Analýza stabilních izotopů uhlíku a dusíku
Vzorky práškových kostí (~100 mg) z lebek MCE1356, 18 belug a 18 narvalů byly ošetřeny 10 ml chloroformu/methanolu v poměru 2:1 (v/v) za sonikace po dobu 1 hodiny, aby se odstranily lipidy. Po odstranění rozpouštědla byly vzorky sušeny za normálního atmosférického tlaku po dobu 24 h. Poté byly vzorky demineralizovány v 10 ml 0,5 M HCl po dobu 4 h za současného míchání orbitální třepačkou. Po demineralizaci byly vzorky opláchnuty vodou typu I do neutrální polohy a poté zahřívány při 75 °C po dobu 36 h v 10-3 M HCl, aby došlo k rozpouštění kolagenu. Ve vodě rozpustný kolagen byl poté vysušen mrazem. Kolagenový extrakt byl poté ošetřen 10 ml chloroformu/methanolu/vody v poměru 10:5:4 (v/v/v) za sonikace po dobu 1 h, aby se odstranily veškeré zbytkové lipidy. Po odstředění byla vrstva chloroformu/methanolu odstraněna a methanol byl odpařen z vodní vrstvy při 60 °C po dobu 24 h. Vzorky byly opět vysušeny mrazem a zváženy do cínových kapslí pro prvkovou a izotopovou analýzu. Stabilní izotopové a prvkové složení uhlíku a dusíku bylo stanoveno pomocí kontinuálního průtokového hmotnostního spektrometru Nu Horizon (Nu Instruments, UK) na Trent University v Kanadě. Stabilní izotopové složení uhlíku a dusíku bylo kalibrováno vzhledem ke stupnicím VPDB a AIR pomocí USGS40 a USGS41a. Standardní nejistota byla stanovena na ±0,17 ‰ pro δ13C a ±0,22 ‰ pro δ15N31; další analytické podrobnosti jsou uvedeny v doplňkovém dokumentu S4. Statistická významnost rozdílů mezi izotopovým složením běluhy a narvala byla posouzena pomocí nepárového t-testu.
.