Bemonstering
Wij analyseerden (i) genoombrede DNA-gegevens geëxtraheerd uit tanden van de schedel (MCE1356, Fig. 2b, aanvullende Fig. S1) en van weefselmonsters van acht beluga’s en acht narwallen bemonsterd in Disko Bay, West-Groenland, en (ii) stabiele koolstof- en stikstofisotopensamenstellingen van botcollageen van MCE1356 en een referentiepanel van 18 beluga’s en 18 narwalschedels van West-Groenland (Fig. 1c). Weefselmonsters (huid) werden bewaard in 96% ethanol. De monsters werden verzameld door wetenschappers van het Groenlands Natuurinstituut in het kader van de algemene vergunning voor biologische bemonstering van de Inuit van de Groenlandse regering en naar Denemarken uitgevoerd met de CITES-vergunningen IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 en 116.2006. Informatie over de monsters is te vinden in de aanvullende tabel S2.
DNA-extractie, bibliotheekvoorbereiding en sequencing
Weefselmonsters
DNA van weefselmonsters werd geëxtraheerd met de Qiagen Blood and Tissue Kit volgens het protocol van de fabrikant. Het DNA werd gefragmenteerd met behulp van de Covaris M220 Focused-ultrasoonator om ~ 350-550 basispaar (bp) fragment lengtes te creëren. Bibliotheken werden gebouwd van de gefragmenteerde DNA-extracten met Illumina NeoPrep volgens de NeoPrep Library Prep System Guide met toepassing van standaardinstellingen. PCR amplificatie, kwantificering, en normalisatie werden alle uitgevoerd door de NeoPrep Library Prep System. De bibliotheken werden gescreend voor grootteverdeling op een Agilent 2100 Bioanalyzer en gepoold in equimolaire verhoudingen voor sequencing op een Illumina HiSeq 2500 met 80 bp SE technology.
Tanden / botmonsters
In tegenstelling tot de weefselmonsters, oude botten en tanden monsters hebben relatief lage DNA-concentraties, die verschillende extractie en bibliotheek te bouwen protocollen noodzakelijk maakt. Ongeveer 0,5 g botpoeder van vijf tanden en een botsplinter van specimen MCE1356 werd geboord met een hand-held Dremel. DNA werd uit het tand/botspoeder geëxtraheerd in een speciaal laboratorium voor het reinigen van oud DNA in het Natuurhistorisch Museum van Denemarken, Universiteit van Kopenhagen, met behulp van de extractiebuffer beschreven in11 met toevoeging van een predigestfase van 30 minuten12. In plaats van Zymo-Spin V-kolommen (Zymo Research) te gebruiken, werd de extractiebuffer geconcentreerd met Amicon Ultra 30 kDa Centrifugaal Filtereenheden en verder geconcentreerd en gereinigd met Qiagen Minelute-buisjes. De gezuiverde extracten werden vervolgens ingebouwd in Illumina bibliotheken volgens het protocol beschreven door13. We gebruikten qPCR om te controleren of de bibliotheek te bouwen succesvol was, om te selecteren welke bibliotheken sequentie, en het juiste aantal PCR-cycli nodig zijn om voldoende amplificatie van elke bibliotheek te berekenen zonder dat overamplificatie. In totaal werden vier bibliotheken geamplificeerd met unieke 6 bp indexen, en gescreend op endogene inhoud op het Illumina MiSeq platform met behulp van 250 bp SE sequencing. De beste bibliotheken werden vervolgens opnieuw gesequenced op de Illumina HiSeq 2500 platform met behulp van 80 bp SE technology.
Bioinformatic analyse
Alle mapping en DNA-beschadiging analyses werden uitgevoerd binnen de Paleomix pijplijn 1.2.1214. Leest werden getrimd met AdapterRemoval 2.2.015 met behulp van standaardinstellingen, behalve minimale leeslengte die werd ingesteld op 25 bp. Lezingen werden geïnspecteerd met FastQC en uitgelijnd met BWA16, waarbij het Backtrack algoritme werd toegepast en de start seed lengte werd uitgeschakeld. Als de gelezen op meerdere posities werden gemapt of als de mapping quality scores (MAPQ score van BWA) minder dan 30 waren, werden ze verwijderd met SAMtools17. Sequence duplicaten werden verwijderd met MarkDuplicates van Picard (beschikbaar via: http://broadinstitute.github.io/picard) en de uiteindelijke uitlijning werd opnieuw uitgelijnd rond indels met behulp van GATK18. De deaminatie van cytosine naar uracil in specimen MCE1356 werd beoordeeld met behulp van de output van mapDamage v2.0.619. MapDamage resultaten toonden geen duidelijk signaal van deaminatie in de sequenties (Supplementary Fig. S3).
Mitochondriale analyse
Om de maternale afstamming van MCE1356 te bepalen, werden DNA-sequencing lezingen in kaart gebracht met de beluga (GenBank toetreding: KY444734) en narwal (GenBank toetreding: NC_005279) mitochondriale referentie genomen en de gemiddelde dekking werd vergeleken. Lezingen van de acht beluga en acht narwal monsters die het referentiepanel vormden, werden in kaart gebracht aan hun respectievelijke mitochondriale referentie genomen. We construeerden mitochondriale consensus sequenties van MCE1356 en de 16 referentie panel monsters met regio’s bedekt door meer dan vijf gelezen met behulp van BEDtools20, SAMtools17 en GATK18. We creëerden twee sequentie-uitlijningen met ClustalW21, met standaardinstellingen, waarin de 16 referentiepanelmonsters waren opgenomen en ofwel de versie van MCE1356 die aan het beluga mitochondriale referentiegenoom was gekoppeld, ofwel de versie van MCE1356 die aan het narwal mitochondriale referentiegenoom was gekoppeld. We gebruikten de twee alignments om mediane-joining haplotype netwerken22 te construeren met behulp van Popart 1.723 (beschikbaar via: http://popart.otago.ac.nz), met uitsluiting van sites met indels of ontbrekende gegevens. Vervolgens werden zowel de post-mapping dekking als de twee haplotype netwerken gebruikt om de soort van MCE1356′s moederlijn te bepalen.
Nucleaire DNA-analyses
We brachten de DNA-sequencing leest van alle monsters in kaart met het orka-referentiegenoom (Orcinus orca) (GCA_000331955.2). Een hoge-kwaliteit beluga walvis genoom werd onlangs gepubliceerd24, maar mapping aan een van de twee potentiële ouderlijke soorten zou kunnen leiden tot vertekeningen in de analyses. Daarom hebben we de reads in kaart gebracht op het genoom van de orka, omdat het goed geassembleerd is en orka’s nog steeds relatief nauw verwant zijn aan beluga’s en narwallen (divergentieperiode van 12 MYA)25, maar toch ver genoeg om het risico van introgressie, die onze analyses zou compliceren, te minimaliseren.
Voor alle verdere filtering hebben we ANGSD v0.92326 gebruikt, een softwarepakket dat genotype waarschijnlijkheden gebruikt in plaats van opgeroepen genotypen, wat vooral nuttig is bij het analyseren van NGS-gegevens met een lage dekkingsgraad. We gebruikten de SAMtools17 methode geïmplementeerd in ANGSD om genotype waarschijnlijkheden te schatten, en leidden major en minor allelen direct uit de genotype waarschijnlijkheden af met behulp van een maximum likelihood aanpak zoals beschreven in27.
Om de in kaart gebrachte gegevens te visualiseren, hebben we de leesdiepteverdeling van alle individuen uitgezet, met uitsluiting van sites met meer dan twee allelen en sites met een Phred score lager dan 25. Het visualiseren van de gegevens van alle individuen samen stelde ons in staat om de gemiddelde leesdiepte van 4.14x te schatten en sites met een verhoogde leesdiepte te identificeren. Dergelijke sites waren meer waarschijnlijk afkomstig van paraloga en repetitieve regio’s van het genoom. De dataset werd visueel geïnspecteerd en verder gefilterd, waarbij alle sites met een leesdiepte groter dan negen (6,9% van de sites) werden uitgesloten. In ANGSD werden SNP’s genoemd op basis van hun allelfrequenties. De minor allelfrequentie (MAF) werd geschat uit de genotype waarschijnlijkheden en een likelihood ratio test werd gebruikt om te bepalen of de MAF verschillend was van nul. Indien de p-waarde van de likelihood ratio test <1e-4 was, werd de site als polymorf beschouwd en in de dataset opgenomen. Het toepassen van deze SNP-filter betekende dat geen sites met minder dan vier reads werden behouden, aangezien sites met minder reads niet zulke lage p-waarden konden krijgen. Deze filters werden toegepast op een dataset met alle 17 monsters, en een dataset zonder MCE1356.
Om te bepalen of MCE1356 een beluga/narwhal hybride was, hebben we de dataset van 17 individuen verder gefilterd, met uitsluiting van sites met geen enkele lezing in MCE1356. Op dit punt was de gemiddelde leesdiepte van MCE1356 slechts 1.15x, dus om er zeker van te zijn dat we geen multicopy loci analyseerden, sloten we sites uit die door meer dan één lees in MCE1356 werden gedekt.
Het was ons doel om de allelen gevonden in MCE1356 te vergelijken met de allelen in het referentiepanel, dus schatten we de waarschijnlijkheid van toewijzing van het allel gevonden in MCE1356 aan de verkeerde ouderlijke soort gegeven verschillende referentiepanel minimum unieke leesdiepten en allelfrequenties. Unieke leesdiepte werd gedefinieerd als het aantal lezingen dat een specifieke site dekt, waarbij alle lezingen van een uniek individu afkomstig zijn. Deze waarschijnlijkheid P werd berekend als in Vergelijking 1:
Waarbij f(ps) de allelfrequentie in de ouderlijke soort is, en urd(psp) de soortspecifieke unieke leesdiepte in het referentiepanel is.
De allelfrequentie van de ouderlijke soort die de hoogste waarschijnlijkheid f(ps-max) oplevert, kan worden beschreven als in vergelijking 2:
Door vergelijking 2 in vergelijking 1 in te voegen, werd de maximale waarschijnlijkheid P(max) van toewijzing van het allel gevonden in MCE1356 aan de verkeerde ouderlijke soort berekend als in vergelijking 3:
Resultaten toonden aan dat met een unieke leesdiepte van twee, drie en vier in elke ouderlijke soort, de maximale waarschijnlijkheid om het allel gevonden in MCE1356 aan de verkeerde ouderlijke soort toe te wijzen 0,148, 0,105 en 0,5 was.148, 0,105 en 0,082, respectievelijk. Deze maximumwaarden mogen niet als foutenpercentages worden geïnterpreteerd, aangezien zij alleen zouden worden verkregen als alle sites binnen de ouderlijke soorten variabel waren, en alle sites een MAF van precies (1 – (urd/urd + 1)) hadden. Bovendien, in de veronderstelling dat de MAF-verdelingen in beluga’s en narwallen gelijkaardig zijn, zou een foutieve toewijzing van allelen gevonden in MCE1356 beide soorten in gelijke mate beïnvloeden, en bijgevolg een minimale invloed hebben op de gevolgtrekkingen van hybridisatie. Een voordeel van het gebruik van de unieke leesdiepte in vergelijkingen 2 en 3 was dat het aantal individuen en het aantal lezingen gecombineerd werden in de schatting van de waarschijnlijkheid van toewijzing van het allel gevonden in MCE1356 aan de verkeerde ouderlijke soort. Vervolgens werden de leesdiepteverdeling, MAF (met vast major en minor allel), en het aantal individuen met gegevens in elke site afzonderlijk berekend voor elke ouderlijke soort.
Drie datasets werden geconstrueerd, die naast MCE1356 panels van ouderlijke soorten omvatten met minimale unieke leesdiepten van twee, drie en vier. Het aantal behouden sites in de drie datasets was respectievelijk 61.105, 11.764 en 360. Als compromis tussen het maximaliseren van het aantal sites en het minimaliseren van de verkeerde toewijzing van allelen gevonden in MCE1356, kozen wij ervoor om de dataset te gebruiken met één lees in MCE1356 en minimale unieke leesdiepten van drie in elke ouderlijke soort. Die dataset omvatte 11.764 sites, die in latere analyses werden gebruikt.
Samenvattende statistieken werden uitgevoerd op de dataset met 17 individuen, waaronder het aantal sites dat (i) gefixeerd was voor verschillende allelen in de beluga- en narwalsoortenpanels; (ii) polymorf in beluga’s, maar niet in narwals; (iii) polymorf in narwals, maar niet in beluga’s; (iv) polymorf in zowel beluga’s als narwals. De sites waarvan wordt aangenomen dat ze gefixeerd zijn tussen de panels van de twee ouderlijke soorten, zullen worden verrijkt met merkers die sterk gedifferentieerd zijn, d.w.z. grote allelfrequentieverschillen vertonen, tussen de twee ouderlijke soorten. Deze merkers, hoewel niet noodzakelijkerwijs gefixeerde verschillen tussen de twee ouderlijke soorten, zijn nog steeds zeer informatief voor voorouders in MCE1356.
We gebruikten de genotype waarschijnlijkheden van de dataset zonder MCE1356 om te verifiëren dat de beluga’s en narwallen in het referentiepanel niet zelf recentelijk gemixte individuen waren. We schatten hun individuele vermengingscoëfficiënten terwijl we twee populaties specificeerden (K = 2) met behulp van NgsAdmix28. Honderd runs werden uitgevoerd en het gemiddelde en de standaardafwijking van de vermengingscoëfficiënten werden gebruikt voor de verdere interpretatie. Om te bevestigen dat de filters niet eerder verborgen vermengde genetische profielen in het referentiepanel aan het licht hadden gebracht, werd deze analyse zowel vóór als na de toepassing van de unieke leesdieptefilters uitgevoerd.
Wij analyseerden de vermengingspercentages van MCE1356 met behulp van fastNGSadmix29, waarbij 1.000 bootstraps werden toegepast. fastNGSadmix gebruikt allelfrequenties van referentiepopulaties/soorten en de genotypewaarschijnlijkheid van een enkel individu om zijn vermengingspercentages te schatten. De software heeft bewezen robuust te zijn met NGS-gegevens met een dekking zo laag als 0,00015×29, en was daarom ideaal voor onze studie.
We hebben verder de hybride status van MCE1356 geschat door sites te onderzoeken die gefixeerd zijn voor verschillende allelen in het beluga-panel en het narwal-panel (9.178 sites), en de waargenomen proportie van gelezen die overeenkomen met het beluga-specifieke allel en het narwal-specifieke allel in MCE1356 te vergelijken met de verwachte proporties onder verschillende hybridisatiescenario’s. Om te bepalen hoe goed zeven verschillende hybridisatiescenario’s overeenkwamen met de waargenomen gegevens, berekenden we een Pearson’s Chi-kwadraat goodness of fit statistiek. De teststatistiek wordt berekend zoals in vergelijking 4:
waar Oi en Ei respectievelijk de waargenomen en verwachte tellingen zijn van de allelen die afkomstig zijn van de ouderlijke soort i. Onder de nulhypothese dat het gekozen scenario goed overeenkomt met de waargenomen gegevens, volgt de teststatistiek T een centrale χ2-verdeling. Het scenario dat het best overeenkomt met de waargenomen gegevens zou dus leiden tot de laagste teststatistiek.
Om de zeven verschillende hybridisatiescenario’s verder te onderzoeken, berekenden we de waarschijnlijkheid van de waargenomen allelen op plaatsen die gefixeerd waren voor verschillende allelen in de ouderlijke populaties. Meer bepaald berekenden we de waarschijnlijkheid van de geobserveerde allelen in de hybride MCE1356, onder een binomiaal model voor de overerving van de allelen van de ouderlijke soorten, verder uitgaande van onafhankelijkheid tussen merkers. Wij willen echter opmerken dat de schending van de onafhankelijkheidsaanname nog steeds tot onbevooroordeelde schattingen leidt. Aangenomen dat de hybride MCE1356 is samengesteld uit een deel b van beluga afstamming en (1-b) van narwal afstamming, kunnen we de waarschijnlijkheid van b schrijven als een product van de waarschijnlijkheid op k onafhankelijke sites, gegeven het aantal gelezen nib dat overeenkomt met het beluga allel op site i en nin dat overeenkomt met narwal afstamming, zoals in vergelijking 5.
We berekenden de log-waarschijnlijkheid van b over zijn hele bereik, van 0 tot 1, en vergeleken de log-waarschijnlijkheden voor de zeven verschillende hybridisatiescenario’s.
Seksebepaling
Om het geslacht van MCE1356 en de individuen in de beluga en narwhal referentiepanels te bepalen, onderzochten we de X-chromosoom tot autosomale dekkingsverhoudingen. Dit werd gedaan door het vergelijken van de dekking van scaffolds die vermoedelijk afkomstig zijn van het X-chromosoom, met die van scaffolds die vermoedelijk afkomstig zijn van de autosomen. We bepaalden welke scaffolds afkomstig zijn van het X-chromosoom door het uitlijnen van de orka genoom aan de koe (Bos taurus) X-chromosoom (CM008168.2) met SatsumaSynteny30. Bovendien, om biases die kunnen optreden als gevolg van homologie tussen sommige regio’s van de X-en Y-chromosomen te verwijderen, hebben we verder uitgelijnd de resulterende vermeende X-chromosoom steigers aan het menselijk Y-chromosoom (NC_000024.10) en verwijderd deze steigers van verdere analyse. Vervolgens berekenden we de dekking over de resterende X scaffolds en de vier grootste scaffolds niet uitlijnen naar het X-chromosoom (dat wil zeggen autosomale scaffolds), met behulp van de diepte-functie in SAMtools17. Om te compenseren voor eventuele variatie in de dekking over het genoom, selecteerden we willekeurig 10.000.000 sites, berekenden de gemiddelde dekking voor deze sites, en herhaalden deze stap 100 keer. Voor elk individu, 5% en 95% betrouwbaarheidsintervallen, evenals eerste en derde kwartielen, werden berekend en gebruikt voor verdere interpretatie.
Stabiele koolstof en stikstof isotoop analyse
Gepoederde botmonsters (~ 100 mg) van MCE1356, 18 beluga en 18 narwal schedels werden behandeld met 10 ml 2:1 chloroform/methanol (v / v) onder sonicatie gedurende 1 uur om lipiden te verwijderen. Na verwijdering van het oplosmiddel werden de monsters gedroogd onder normale atmosferische druk gedurende 24 uur. Vervolgens werden de monsters gedemineraliseerd in 10 ml 0,5 M HCl gedurende 4 uur, terwijl ze werden geschud door een rondschudapparaat. Na de demineralisatie werden de monsters tot neutraal gespoeld met type I water, en vervolgens gedurende 36 uur bij 75 °C verwarmd in 10-3 M HCl om het collageen op te lossen. Het in water oplosbare collageen werd vervolgens gevriesdroogd. Het collageenextract werd vervolgens behandeld met 10 ml 10:5:4 chloroform/methanol/water (v/v/v) onder sonicatie gedurende 1 uur om eventueel achtergebleven lipiden te verwijderen. Na centrifugatie werd de chloroform/methanol-laag verwijderd en werd de methanol uit de waterlaag verdampt bij 60 °C gedurende 24 uur. De monsters werden opnieuw gevriesdroogd en afgewogen in tinnen capsules voor elementaire en isotopische analyse. De stabiele isotopische en elementaire samenstellingen van koolstof en stikstof werden bepaald met een Nu Horizon (Nu Instruments, UK) massaspectrometer met continue stroom en isotopenverhouding aan de Trent University, Canada. Stabiele koolstof- en stikstofisotopensamenstellingen werden gekalibreerd ten opzichte van de VPDB- en AIR-schalen met behulp van USGS40 en USGS41a. De standaardonzekerheid werd vastgesteld op ±0,17‰ voor δ13C en ±0,22‰ voor δ15N31; bijkomende analytische details worden verstrekt in het aanvullend document S4. De statistische significantie van verschillen tussen beluga- en narwalisotopensamenstellingen werd beoordeeld met behulp van ongepaarde t-tests.