Sampling
Vi analyserade (i) DNA-data för hela arvsmassan som extraherats från tänder i skallen (MCE1356, Fig. 2b, Supplementary Fig. S1) och från vävnadsprover från åtta belugor och åtta narvaler som provtagits i Diskobukten, västra Grönland, och (ii) stabila kol- och kväveisotopsammansättningar av benkollagen från MCE1356 och en referenspanel med 18 belugor och 18 narvalskallar från västra Grönland (Fig. 1c). Vävnadsprover (hud) förvarades i 96 % etanol. Proverna samlades in av forskare från Grönlands naturresursinstitut enligt Grönlands regerings allmänna tillstånd för biologisk provtagning av inuiterna och exporterades till Danmark enligt CITES-tillstånden IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 och 116.2006. Provinformation finns i detalj i kompletterande tabell S2.
DNA-extraktion, biblioteksberedning och sekvensering
Vävnadsprover
DNA från vävnadsprover extraherades med hjälp av Qiagen Blood and Tissue Kit enligt tillverkarens protokoll. DNA fragmenterades med hjälp av Covaris M220 Focused-ultrasonicator för att skapa fragmentlängder på ~350-550 baspar (bp). Bibliotek byggdes från de fragmenterade DNA-extrakten med hjälp av Illumina NeoPrep enligt NeoPrep Library Prep System Guide med standardinställningar. PCR-amplifiering, kvantifiering och normalisering utfördes av NeoPrep Library Prep System. Biblioteken screenades för storleksfördelning på en Agilent 2100 Bioanalyzer och sammanfördes i ekvimolära förhållanden före sekvensering på en Illumina HiSeq 2500 med 80 bp SE-teknik.
Tand-/benprover
Till skillnad från vävnadsproverna har gamla ben- och tandprover relativt låga DNA-koncentrationer, vilket kräver olika extraktions- och biblioteksbyggnadsprotokoll. Ungefär 0,5 g benpulver från fem tänder och en bensplitter från exemplar MCE1356 borrades med hjälp av en handhållen Dremel. DNA extraherades från tand- och benpulvret i ett särskilt laboratorium för rening av forntida DNA vid Danmarks naturhistoriska museum, Köpenhamns universitet, med hjälp av den extraktionsbuffert som beskrivs i11 med tillägg av ett 30 minuter långt förfördelningsskede12. I stället för att använda Zymo-Spin V-kolonner (Zymo Research) koncentrerades extraktionsbufferten med Amicon Ultra 30 kDa centrifugalfilterenheter och koncentrerades och rengjordes ytterligare med Qiagen Minelute-rör. De renade extrakten byggdes sedan in i Illumina-bibliotek enligt det protokoll som beskrivs i13. Vi använde qPCR för att kontrollera att biblioteksuppbyggnaden var framgångsrik, för att välja ut vilka bibliotek som skulle sekvenseras och för att beräkna det lämpliga antalet PCR-cykler som krävs för att amplifiera varje bibliotek tillräckligt utan att orsaka överamplifiering. Totalt amplifierades fyra bibliotek med unika index på 6 bp och undersöktes för endogent innehåll på Illumina MiSeq-plattformen med hjälp av 250 bp SE-sekvensering. De bästa biblioteken sekvenserades sedan på nytt på Illumina HiSeq 2500-plattformen med hjälp av 80 bp SE-teknik.
Bioinformatisk analys
Alla kartläggningar och analyser av DNA-skador utfördes inom Paleomix pipeline 1.2.1214. Läsningar trimmades med AdapterRemoval 2.2.015 med standardinställningar utom minsta läsningslängd som sattes till 25 bp. Läsningarna inspekterades med FastQC och anpassades med BWA16 genom att tillämpa Backtrack-algoritmen och inaktivera startlängden för frön. Om läsningar mappades till flera positioner eller hade mappningskvalitetspoäng (MAPQ-poäng från BWA) som var lägre än 30, togs de bort med hjälp av SAMtools17. Sekvensdubbletter togs bort med hjälp av MarkDuplicates från Picard (tillgänglig från: http://broadinstitute.github.io/picard) och den slutliga anpassningen omjusterades runt indels med hjälp av GATK18. Deaminering av cytosin till uracil i prov MCE1356 bedömdes med hjälp av resultatet från mapDamage v2.0.619. MapDamage-resultaten visade inte på någon tydlig signal av deaminering i sekvenserna (Supplementary Fig. S3).
Mitokondriell analys
För att fastställa den maternella härstamningen av MCE1356 kartlades DNA-sekvenseringsavläsningar till belugas (GenBank accession: KY444734) och narvalens (GenBank accession: NC_005279) mitokondriella referensgenom och den genomsnittliga täckningen jämfördes. Läsningar från de åtta beluga- och åtta narvalprover som ingick i referenspanelen kartlades till deras respektive mitokondriella referensgenom. Vi konstruerade mitokondriella konsensussekvenser av MCE1356 och de 16 referenspanelproverna med regioner som täcktes av mer än fem läsningar med hjälp av BEDtools20, SAMtools17 och GATK18. Vi skapade två sekvensanpassningar med ClustalW21 , med standardinställningar, som inkluderade de 16 referenspanelproverna och antingen den version av MCE1356 som kartlagts till belugas mitokondriella referensgenom eller den version av MCE1356 som kartlagts till narvalens mitokondriella referensgenom. Vi använde de två anpassningarna för att konstruera median-joining haplotypnätverk22 med hjälp av Popart 1.723 (tillgänglig från: http://popart.otago.ac.nz), och uteslöt alla platser med indels eller saknade data. Därefter användes både täckningen efter kartläggningen och de två haplotypnätverken för att bestämma arten i MCE1356′s moderliga släktlinje.
Kärn-DNA-analyser
Vi kartlade DNA-sekvenseringsavläsningarna från alla prover till späckhuggarnas (Orcinus orca) referensgenom (GCA_000331955.2). Ett högkvalitativt belugavalsgenom publicerades nyligen24 , men mappning till en av de två potentiella föräldraarterna skulle kunna skapa bias i analyserna. Därför mappade vi läsningarna till späckhuggarens genom, eftersom det är väl sammansatt och späckhuggare fortfarande är relativt nära besläktade med belugas och narvalar (divergenstid på 12 MYA)25, men ändå tillräckligt avlägsna för att minimera risken för introgression som skulle komplicera våra analyser.
För all ytterligare filtrering använde vi ANGSD v0.92326, ett programvarupaket som använder genotypsannolikheter i stället för kallade genotyper, vilket är särskilt användbart när man analyserar NGS-data med låg täckning. Vi använde SAMtools17-metoden som är implementerad i ANGSD för att uppskatta genotypsannolikheter, och härledde major- och minoralleler direkt från genotypsannolikheterna med hjälp av en maximal sannolikhetsmetod enligt beskrivningen i27.
För att visualisera de mappade uppgifterna plottade vi fördelningen av läsdjupet från alla individer, med uteslutande av platser med mer än två alleler och platser med en Phred-poäng under 25. Genom att visualisera data från alla individer tillsammans kunde vi uppskatta det genomsnittliga läsdjupet på 4,14x och identifiera platser med förhöjt läsdjup. Det var mer sannolikt att sådana platser kom från paraloger och repetitiva områden i genomet. Datamängden inspekterades visuellt och filtrerades ytterligare, varvid alla platser med ett läsdjup som var större än nio (6,9 % av platserna) exkluderades. I ANGSD kallades SNPs utifrån deras allelfrekvenser. Den mindre allelfrekvensen (MAF) uppskattades från genotypens sannolikheter och ett likelihood ratio-test användes för att avgöra om MAF var annorlunda än noll. Om p-värdet från likelihood ratio-testet var <1e-4 betraktades platsen som polymorf och behölls i datasetet. Tillämpning av detta SNP-filter innebar att inga platser med färre än fyra reads behölls, eftersom platser som täcks av färre reads inte kunde få så låga p-värden. Dessa filter tillämpades på en datauppsättning med alla 17 prover och en datauppsättning utan MCE1356.
För att avgöra om MCE1356 var en beluga/narval-hybrid filtrerade vi ytterligare datauppsättningen med 17 individer och exkluderade platser med inga reads i MCE1356. Vid denna tidpunkt var det genomsnittliga avläsningsdjupet i MCE1356 endast 1,15x, så för att säkerställa att vi inte analyserade multikopiska loci exkluderade vi platser som täcktes av mer än en avläsning i MCE1356.
Vårt mål var att jämföra de alleler som hittades i MCE1356 med allelerna i referenspanelen, så vi uppskattade sannolikheten för att tilldela allelen som hittades i MCE1356 till fel föräldraart med tanke på olika referenspanelers minimala unika avläsningsdjup och allelfrekvenser. Unikt avläsningsdjup definierades som antalet avläsningar som täcker en specifik plats, där alla avläsningar kom från en unik individ. Denna sannolikhet P beräknades enligt ekvation 1:
Varvid f(ps) är allelfrekvensen i den föräldrade arten och urd(psp) är det artspecifika unika läsdjupet i referenspanelen.
Den allelfrekvens hos föräldraarten som ger den högsta sannolikheten f(ps-max) kan beskrivas enligt ekvation 2:
Där ekvation 2 infogas i ekvation 1, beräknades den maximala sannolikheten P(max) för att den allel som hittades i MCE1356 skulle hänföras till fel föräldraart enligt ekvation 3:
Resultaten visade att med ett unikt läsdjup på två, tre och fyra i varje föräldraart, den maximala sannolikheten för att tilldela den allel som hittades i MCE1356 till fel föräldraart var 0.148, 0,105 respektive 0,082. Dessa maximala värden bör inte tolkas som felprocent, eftersom de endast skulle erhållas om alla platser var variabla inom föräldraarterna och alla platser hade en MAF på exakt (1 – (urd/urd + 1)). Om man dessutom antar att MAF-fördelningarna hos belugor och narvaler var likartade, skulle felaktig tilldelning av alleler som hittats i MCE1356 påverka båda arterna lika mycket och därför ha en minimal inverkan på slutsatserna om hybridisering. En fördel med att använda det unika läsdjupet i ekvationerna 2 och 3 var att det kombinerade antalet individer och antalet läsuppgifter i uppskattningen av sannolikheten för att tilldela den allel som hittades i MCE1356 till fel föräldraart. Därefter beräknades fördelningen av avläsningsdjup, MAF (med fast major- och minorallel) och antalet individer med data på varje plats separat för varje föräldraart.
Tre dataset konstruerades, som förutom MCE1356 innehöll föräldraartspaneler med minsta unika avläsningsdjup på två, tre och fyra. Antalet platser som bevarades i de tre datamängderna var 61 105, 11 764 respektive 360. Som en kompromiss mellan att maximera antalet platser och minimera den felaktiga tilldelningen av alleler som hittades i MCE1356, valde vi att använda datasetet med en avläsning i MCE1356 och minsta unika avläsningsdjup på tre i varje föräldraart. Det datasetet innehöll 11 764 platser, som användes i efterföljande analyser.
Samlingsstatistik utfördes på datasetetet med 17 individer, inklusive antal platser som var (i) fixerade för olika alleler i panelerna för beluga- och narvalarterna; (ii) polymorfa hos belugor, men inte hos narvalar; (iii) polymorfa hos narvalar, men inte hos belugor; (iv) polymorfa hos både belugor och narvalar. De platser som uppskattas vara fasta mellan de två föräldraarterna kommer att berikas med markörer som är starkt differentierade, dvs. har stora skillnader i allelfrekvens, mellan de två föräldraarterna. Dessa markörer, även om de inte nödvändigtvis är fasta skillnader mellan de två föräldraarterna, är fortfarande mycket informativa för härstamning i MCE1356.
Vi använde genotypsannolikheterna från datasetet utan MCE1356 för att verifiera att belugorna och narvalarna i referenspanelen inte själva var nyligen blandade individer. Vi uppskattade deras individuella blandningskoefficienter samtidigt som vi specificerade två populationer (K = 2) med hjälp av NgsAdmix28. Hundra körningar utfördes och medelvärde och standardavvikelse för blandningskoefficienterna användes för efterföljande tolkning. För att bekräfta att filtren inte hade avslöjat tidigare dolda blandade genetiska profiler i referenspanelen utfördes denna analys både före och efter det att filtren för unikt läsdjup tillämpades.
Vi analyserade blandningskvoterna för MCE1356 med hjälp av fastNGSadmix29 och tillämpade 1 000 bootstraps. fastNGSadmix använder allelfrekvenser för referenspopulationer/arter och genotypsannolikheterna för en enskild individ för att skatta dess blandningskvoter. Programvaran har visat sig vara robust med NGS-data med så låg täckning som 0,00015×29 och var därför idealisk för vår studie.
Vi skattade vidare hybridstatusen hos MCE1356 genom att undersöka platser som är fixerade för olika alleler i beluga-panelen och narval-panelen (9 178 platser), och genom att jämföra den observerade andelen läsningar som matchar den belugaspecifika allelen och den narvalspecifika allelen i MCE1356 med de förväntade proportionerna under olika hybridiseringsscenarier. För att avgöra hur väl sju olika hybridiseringsscenarier stämde överens med de observerade uppgifterna beräknade vi en Pearsons Chi-square goodness of fit-statistik. Teststatistiken beräknas enligt ekvation 4:
där Oi och Ei är det observerade och förväntade antalet alleler som härrör från föräldraarten i, respektive. Under nollhypotesen där det valda scenariot stämmer väl överens med de observerade uppgifterna följer teststatistiken T en central χ2-fördelning. Således skulle det scenario som bäst motsvarar de observerade uppgifterna leda till den lägsta teststatistiken.
För att ytterligare undersöka de sju olika hybridiseringsscenarierna beräknade vi sannolikheten för de observerade allelerna på platser som var fixerade för olika alleler i föräldrapopulationerna. Specifikt beräknade vi sannolikheten för de observerade allelerna i hybriden MCE1356, enligt en binomialmodell för nedärvningen av allelerna från föräldraarterna, med ytterligare antagande om oberoende mellan markörerna. Vi vill dock notera att brottet mot antagandet om oberoende fortfarande leder till opartiska skattningar. Om vi antar att hybriden MCE1356 består av en andel b av belugaförstamning och (1-b) av narvalförstamning, kan vi skriva sannolikheten för b som en produkt av sannolikheten på k oberoende platser, givet antalet läsningar nib som matchar belugaallelen på plats i och nin som matchar narvalförstamning, enligt ekvation 5.
Vi beräknade log-likelihood för b över hela dess intervall, från 0 till 1, och jämförde log-vannolikheterna mellan de sju olika hybridiseringsscenarierna.
Könsbestämning
För att bestämma könet hos MCE1356 och individerna i referenspanelerna för beluga och narval undersökte vi förhållandet mellan X-kromosom och autosomal täckning. Detta gjordes genom att jämföra täckningen över scaffolds som antas härröra från X-kromosomen med täckningen över scaffolds som antas härröra från autosomerna. Vi fastställde vilka scaffolds som kommer från X-kromosomen genom att anpassa späckhuggarens genom till koens (Bos taurus) X-kromosom (CM008168.2) med SatsumaSynteny30. För att avlägsna de förvrängningar som kan uppstå på grund av homologin mellan vissa områden av X- och Y-kromosomerna, anpassade vi dessutom de resulterande putativa X-kromosom-ställningar till människans Y-kromosom (NC_000024.10) och avlägsnade dessa ställningar från vidare analys. Vi beräknade sedan täckningen över de återstående X-ställningar och de fyra största ställningarna som inte anpassats till X-kromosomen (dvs. autosomala ställningar) med hjälp av djupfunktionen i SAMtools17. För att kompensera för eventuella variationer i täckningen över hela genomet valde vi slumpmässigt ut 10 000 000 platser, beräknade den genomsnittliga täckningen för dessa platser och upprepade detta steg 100 gånger. För varje individ beräknades 5 % och 95 % konfidensintervall samt första och tredje kvartilen och användes för efterföljande tolkning.
Stabil kol- och kväveisotopanalys
Pulveriserade benprover (~100 mg) från MCE1356, 18 belugakranier och 18 narvalskranier behandlades med 10 ml 2:1 kloroform/metanol (v/v) under sonikation i 1 timme för att avlägsna lipider. Efter att lösningsmedlet avlägsnats torkades proverna under normalt atmosfärstryck i 24 timmar. Proverna avmineraliserades sedan i 10 ml 0,5 M HCl i 4 timmar under omrörning med en omrörare. Efter demineraliseringen sköljdes proverna till neutralitet med vatten av typ I och värmdes sedan vid 75 °C i 36 timmar i 10-3 M HCl för att lösa upp kollagenet. Det vattenlösliga kollagenet frystorkades sedan. Kollagenextraktet behandlades sedan med 10 ml 10:5:4 kloroform/metanol/vatten (v/v/v/v) under sonikation i 1 timme för att avlägsna eventuella kvarvarande lipider. Efter centrifugering avlägsnades kloroform/metanol-skiktet och metanolen avdunstades från vattenskiktet vid 60 °C i 24 timmar. Proverna frystorkades återigen och vägdes in i tennkapslar för elementär- och isotopanalys. Stabila isotop- och grundämneskompositioner av kol och kväve bestämdes med hjälp av en Nu Horizon (Nu Instruments, UK) masspektrometer med kontinuerligt flöde för isotopförhållande vid Trent University, Kanada. Stabila kol- och kväveisotopsammansättningar kalibrerades i förhållande till VPDB- och AIR-skalorna med hjälp av USGS40 och USGS41a. Standardosäkerheten fastställdes till ±0,17‰ för δ13C och ±0,22‰ för δ15N31; ytterligare analytiska detaljer finns i tilläggsdokument S4. Den statistiska betydelsen av skillnaderna mellan belugas och narvalars isotopkompositioner bedömdes med hjälp av opartade t-test.