Sampling
Vi analyserede (i) DNA-data for hele genomet, som blev udtrukket fra tænder i kraniet (MCE1356, Fig. 2b, Supplerende fig. S1) og fra vævsprøver fra otte belugas og otte narhvaler, der blev udtaget i Diskobugten, Vestgrønland, og (ii) stabile kulstof- og kvælstofisotopersammensætninger af knoglekollagen fra MCE1356 og et referencepanel af 18 belugas og 18 narhvalskranier fra Vestgrønland (fig. 1c). Vævsprøver (hud) blev opbevaret i 96 % ethanol. Prøverne blev indsamlet af forskere fra Grønlands Naturinstitut i henhold til den generelle tilladelse til biologisk prøvetagning af inuitter fra Grønlands regering og eksporteret til Danmark i henhold til CITES-tilladelse IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 og 116.2006. Prøveoplysninger er detaljeret beskrevet i supplerende tabel S2.
DNA-ekstraktion, bibliotekspræparation og sekventering
Vævsprøver
DNA fra vævsprøver blev ekstraheret ved hjælp af Qiagen Blood and Tissue Kit i henhold til producentens protokol. DNA’et blev fragmenteret ved hjælp af Covaris M220 Focused-ultrasonicator for at skabe fragmentlængder på ~350-550 basepar (bp). Der blev opbygget biblioteker fra de fragmenterede DNA-ekstrakter ved hjælp af Illumina NeoPrep i henhold til NeoPrep Library Prep System Guide med standardindstillinger. PCR-amplifikation, kvantificering og normalisering blev alle udført af NeoPrep Library Prep System. Bibliotekerne blev screenet for størrelsesfordeling på en Agilent 2100 Bioanalyzer og puljet i ækimolære forhold før sekventering på en Illumina HiSeq 2500 med 80 bp SE-teknologi.
Tand/knogleprøver
I modsætning til vævsprøverne har gamle knogle- og tandprøver relativt lave DNA-koncentrationer, hvilket nødvendiggør forskellige ekstraktions- og biblioteksopbygningsprotokoller. Ca. 0,5 g knoglepulver fra fem tænder og et knogleskår fra prøve MCE1356 blev boret ved hjælp af en håndholdt Dremel. DNA blev ekstraheret fra tand- og knoglepulveret i et særligt laboratorium til rensning af oldtids-DNA på Danmarks Naturhistoriske Museum, Københavns Universitet, ved hjælp af den ekstraktionsbuffer, der er beskrevet i11 , med tilføjelse af en 30 minutters forfordigestningsfase12. I stedet for at anvende Zymo-Spin V-kolonner (Zymo Research) blev ekstraktionsbufferen koncentreret ved hjælp af Amicon Ultra 30 kDa centrifugalfilterenheder og yderligere koncentreret og renset ved hjælp af Qiagen Minelute-rør. De rensede ekstrakter blev derefter opbygget til Illumina-biblioteker efter den protokol, der er beskrevet i13. Vi brugte qPCR til at kontrollere, at opbygningen af biblioteket var vellykket, til at vælge, hvilke biblioteker der skulle sekventeres, og til at beregne det passende antal PCR-cyklusser, der var nødvendige for at amplificere hvert bibliotek tilstrækkeligt uden at forårsage overamplifikation. I alt blev fire biblioteker amplificeret med unikke 6 bp-indekser og screenet for endogent indhold på Illumina MiSeq-platformen ved hjælp af 250 bp SE-sekventering. De bedste biblioteker blev derefter re-seksekventeret på Illumina HiSeq 2500-platformen ved hjælp af 80 bp SE-teknologi.
Bioinformatisk analyse
Alle kortlægnings- og DNA-skadeanalyser blev udført inden for Paleomix-pipeline 1.2.1214. Læsninger blev trimmet med AdapterRemoval 2.2.015 ved hjælp af standardindstillinger undtagen mindste læselængde, som blev sat til 25 bp. Læsninger blev inspiceret ved hjælp af FastQC og justeret med BWA16 ved hjælp af Backtrack-algoritmen og deaktivering af startfrølængden. Hvis læsninger blev kortlagt til flere positioner eller havde kortlægningskvalitetsscorer (MAPQ-score fra BWA) på under 30, blev de fjernet ved hjælp af SAMtools17. Sekvensduplikater blev fjernet ved hjælp af MarkDuplicates fra Picard (tilgængelig fra: http://broadinstitute.github.io/picard), og den endelige alignment blev justeret omkring indels ved hjælp af GATK18. Deaminering af cytosin til uracil i prøve MCE1356 blev vurderet ved hjælp af output fra mapDamage v2.0.619. MapDamage-resultaterne viste ikke noget klart signal af deaminering i sekvenserne (Supplerende Fig. S3).
Mitokondrieanalyse
For at bestemme den maternelle afstamning af MCE1356 blev DNA-sekventeringslæsninger kortlagt til beluga (GenBank accession: KY444734) og narhval (GenBank accession: NC_005279) mitokondrielle referencegenomer, og den gennemsnitlige dækning blev sammenlignet. Læsninger fra de otte beluga- og otte narhvalprøver, der udgjorde referencepanelet, blev kortlagt til deres respektive mitokondrielle referencegenomer. Vi konstruerede mitokondrielle konsensussekvenser af MCE1356 og de 16 prøver fra referencepanelet med regioner, der var dækket af mere end fem læsninger, ved hjælp af BEDtools20, SAMtools17 og GATK18. Vi skabte to sekvensudligninger ved hjælp af ClustalW21 med standardindstillinger, som omfattede de 16 referencepanelprøver og enten den version af MCE1356, der er kortlagt til belugas mitokondrielle referencegenom, eller den version af MCE1356, der er kortlagt til narhvalens mitokondrielle referencegenom. Vi brugte de to alignments til at konstruere median-joining haplotypenetværk22 ved hjælp af Popart 1.723 (tilgængelig fra: http://popart.otago.ac.nz), idet vi udelukkede alle steder med indels eller manglende data. Efterfølgende blev både dækning efter kortlægningen og de to haplotypenetværk brugt til at bestemme arten af MCE1356′s moderlige slægt.
Nukleare DNA-analyser
Vi kortlagde DNA-sekventeringslæsningerne fra alle prøverne til spækhuggerhvalens (Orcinus orca) referencegenom (GCA_000331955.2). Der blev for nylig offentliggjort et beluga-hvalgenom af høj kvalitet24 , men kortlægning til en af de to potentielle forældrearter kunne skabe skævheder i analyserne. Derfor kortlagde vi læserne til spækhuggerhvalgenomet, da det er godt samlet, og spækhuggere stadig er relativt tæt beslægtet med belugas og narhvaler (divergenstid på 12 MYA)25, men alligevel langt nok fra hinanden til at minimere risikoen for introgression, der ville komplicere vores analyser.
Til al yderligere filtrering brugte vi ANGSD v0.92326, en softwarepakke, der bruger genotypesandsynligheder i stedet for kaldte genotyper, hvilket er særligt nyttigt ved analyse af NGS-data med lav dækning. Vi brugte SAMtools17-metoden, der er implementeret i ANGSD, til at estimere genotype-sandsynligheder, og udledte major- og minoralleler direkte fra genotype-sandsynlighederne ved hjælp af en maximum likelihood-tilgang som beskrevet i27.
For at visualisere de kortlagte data plottede vi læsedybdefordelingen fra alle individer, idet vi udelukkede steder med mere end to alleler og steder med en Phred-score på under 25. Ved at visualisere dataene fra alle individer kombineret kunne vi estimere den gennemsnitlige læsedybde på 4,14x og identificere steder med forhøjet læsedybde. Sådanne steder var mere tilbøjelige til at komme fra paraloger og repetitive regioner af genomet. Datasættet blev visuelt inspiceret og yderligere filtreret, idet alle steder med en læsedybde på mere end ni (6,9 % af stederne) blev udelukket. I ANGSD blev SNP’er kaldt på grundlag af deres allelfrekvenser. Den mindre allelfrekvens (MAF) blev estimeret ud fra genotypesandsynlighederne, og der blev anvendt en likelihood ratio-test for at afgøre, om MAF var forskellig fra nul. Hvis p-værdien fra likelihood ratio-testen var <1e-4, blev stedet betragtet som polymorft og bevaret i datasættet. Anvendelse af dette SNP-filter betød, at ingen steder med mindre end fire læsninger blev bevaret, da steder, der er dækket af færre læsninger, ikke kunne opnå så lave p-værdier. Disse filtre blev anvendt på et datasæt, der omfattede alle 17 prøver, og et datasæt uden MCE1356.
For at afgøre, om MCE1356 var en beluga/narhval-hybrid, filtrerede vi yderligere datasættet med 17 individer, idet vi udelukkede steder med ingen reads i MCE1356. På dette tidspunkt var den gennemsnitlige læsedybde i MCE1356 kun 1,15x, så for at sikre, at vi ikke analyserede multicopy loci, udelukkede vi steder, der var dækket af mere end én læsning i MCE1356.
Vores mål var at sammenligne de alleler, der blev fundet i MCE1356, med allelerne i referencepanelet, så vi estimerede sandsynligheden for at tildele den allel, der blev fundet i MCE1356, til den forkerte forældreart i betragtning af forskellige referencepanelers minimale unikke læsedybder og allelfrekvenser. Den unikke læsedybde blev defineret som antallet af læsninger, der dækker et specifikt sted, hvor alle læsninger kom fra et unikt individ. Denne sandsynlighed P blev beregnet som i ligning 1:
Hvor f(ps) er allelfrekvensen i forældrearten, og urd(psp) er den artsspecifikke unikke læsedybde i referencepanelet.
Den allelfrekvens for forældrearten, der giver den højeste sandsynlighed f(ps-max), kan beskrives som i ligning 2:
Ved indsættelse af ligning 2 i ligning 1, blev den maksimale sandsynlighed P(max) for at henføre den allel, der blev fundet i MCE1356, til den forkerte forældreart beregnet som i ligning 3:
Resultaterne viste, at med en unik læsedybde på to, tre og fire i hver forældreart, var den maksimale sandsynlighed for at tildele den allel, der blev fundet i MCE1356, til den forkerte forældreart 0.148, 0,105 og 0,082, henholdsvis. Disse maksimumsværdier bør ikke fortolkes som fejlprocenter, da de kun ville blive opnået, hvis alle steder var variable inden for forældrearterne, og alle steder havde en MAF på præcis (1 – (urd/urd + 1)). Hvis man desuden antager, at MAF-fordelingerne hos belugas og narhvaler var ens, ville fejlagtig tildeling af alleler fundet i MCE1356 påvirke begge arter lige meget og derfor have en minimal indflydelse på slutninger om hybridisering. En fordel ved at anvende den unikke læsedybde i ligning 2 og 3 var, at den kombinerede antallet af individer og antallet af læsninger i vurderingen af sandsynligheden for at tildele den allel, der blev fundet i MCE1356, til den forkerte forælderart. Efterfølgende blev read-dybdefordelingen, MAF (med fast major- og minorallel) og antallet af individer med data på hvert sted beregnet separat for hver forældreart.
Der blev konstrueret tre datasæt, som ud over MCE1356 omfattede forældreartspaneler med en unik read-dybde på mindst to, tre og fire. Antallet af steder, der blev bevaret i de tre datasæt, var henholdsvis 61 105, 11 764 og 360. Som et kompromis mellem maksimering af antallet af steder og minimering af den forkerte tildeling af alleler fundet i MCE1356 valgte vi at anvende datasættet med én læsning i MCE1356 og en minimum unik læsningsdybde på tre i hver af forældrenes arter. Dette datasæt indeholdt 11.764 steder, som blev anvendt i de efterfølgende analyser.
Summariske statistikker blev udført på datasættet med 17 individer, herunder antallet af steder, der var (i) fikseret for forskellige alleler i beluga- og narhvalartspanelerne; (ii) polymorfe hos belugas, men ikke hos narhvaler; (iii) polymorfe hos narhvaler, men ikke hos belugas; (iv) polymorfe hos både belugas og narhvaler. De steder, der skønnes at være faste mellem de to forældreartspaneler, vil være beriget med markører, der er stærkt differentierede, dvs. som har store forskelle i allelfrekvens mellem de to forældrearter. Disse markører er, selv om de ikke nødvendigvis har faste forskelle mellem de to forældrearter, stadig meget informative for afstamning i MCE1356.
Vi brugte genotype-sandsynlighederne fra datasættet uden MCE1356 til at verificere, at belugas og narhvaler i referencepanelet ikke selv var nyligt blandede individer. Vi estimerede deres individuelle blandingskoefficienter, mens vi specificerede to populationer (K = 2) ved hjælp af NgsAdmix28. Der blev udført 100 kørsler, og middelværdien og standardafvigelsen af blandingskoefficienterne blev anvendt til efterfølgende fortolkning. For at bekræfte, at filtrene ikke havde afsløret tidligere skjulte blandede genetiske profiler i referencepanelet, blev denne analyse udført både før og efter, at de unikke læsedybdefiltre blev anvendt.
Vi analyserede blandingsandelene for MCE1356 ved hjælp af fastNGSadmix29 og anvendte 1.000 bootstraps. fastNGSadmix anvender allelfrekvenser for referencepopulationer/arter og genotypesandsynlighederne for et enkelt individ til at estimere dets blandingsandele. Softwaren har vist sig at være robust over for NGS-data med en dækning så lav som 0,00015×29 og var derfor ideel til vores undersøgelse.
Vi estimerede yderligere hybridstatus for MCE1356 ved at undersøge steder, der er fikseret for forskellige alleler i beluga-panelet og narhval-panelet (9.178 steder), og ved at sammenligne den observerede andel af læsninger, der matcher den beluga-specifikke allel og den narhval-specifikke allel i MCE1356, med de forventede andele under forskellige hybridiseringsscenarier. For at afgøre, hvor godt syv forskellige hybridiseringsscenarier passede til de observerede data, beregnede vi en Pearsons Chi-kvadrat-statistik for tilpasningens godhed. Teststatistikken beregnes som i ligning 4:
hvor Oi og Ei er de observerede og forventede antal alleler, der stammer fra henholdsvis forældreart i. Under nulhypotesen, hvor det valgte scenarie stemmer godt overens med de observerede data, følger teststatistikken T en central χ2 -fordeling. Således vil det scenarie, der svarer bedst til de observerede data, føre til den laveste teststatistik.
For yderligere at undersøge de syv forskellige hybridiseringsscenarier beregnede vi sandsynligheden for de observerede alleler på steder, der var fikseret for forskellige alleler i forældrepopulationerne. Specifikt beregnede vi sandsynligheden for de observerede alleler i hybriden MCE1356 under en binomialmodel for nedarvning af allelerne fra forældrearterne, idet vi yderligere antog uafhængighed mellem markørerne. Vi vil dog gerne bemærke, at bruddet på uafhængighedsantagelsen stadig fører til ufordelagtige estimater. Hvis vi antager, at hybriden MCE1356 består af en andel b af beluga-forfædre og (1-b) af narhval-forfædre, kan vi skrive sandsynligheden for b som et produkt af sandsynligheden på k uafhængige steder, givet antallet af læsninger nib, der passer til beluga-allelen på sted i, og nin, der passer til narhval-forfædre, som i ligning 5.
Vi beregnede log-sandsynligheden for b over hele dens område, fra 0 til 1, og sammenlignede log-sandsynlighederne på tværs af de syv forskellige hybridiseringsscenarier.
Kønsbestemmelse
For at bestemme kønnet på MCE1356 og individerne i beluga- og narhval-referencepanelerne undersøgte vi forholdet mellem X-kromosomer og autosomale dækningsgrader. Dette blev gjort ved at sammenligne dækningen på tværs af stilladser, der formodes at stamme fra X-kromosomet, med dækningen på stilladser, der formodes at stamme fra autosomerne. Vi bestemte, hvilke stilladser der stammer fra X-kromosomet, ved at tilpasse spækhuggerhvalgenomet til koens (Bos taurus) X-kromosom (CM008168.2) med SatsumaSynteny30. For at fjerne skævheder, der kan opstå på grund af homologi mellem visse områder af X- og Y-kromosomerne, tilpassede vi desuden de resulterende formodede X-kromosom-stilladser yderligere til det menneskelige Y-kromosom (NC_000024.10) og fjernede disse stilladser fra den videre analyse. Derefter beregnede vi dækningen på tværs af de resterende X-stilladser og de fire største stilladser, der ikke er tilpasset til X-kromosomet (dvs. autosomale stilladser), ved hjælp af dybdefunktionen i SAMtools17. For at kompensere for eventuelle variationer i dækningen på tværs af genomet valgte vi tilfældigt 10 000 000 steder, beregnede den gennemsnitlige dækning for disse steder og gentog dette trin 100 gange. For hvert individ blev 5 % og 95 % konfidensintervaller samt første og tredje kvartil beregnet og anvendt til efterfølgende fortolkning.
Stabil kulstof- og kvælstofisotopanalyse
Pulveriserede knogleprøver (~100 mg) fra MCE1356, 18 beluga- og 18 narhvalskranier blev behandlet med 10 ml 2:1 kloroform/methanol (v/v) under sonicering i 1 time for at fjerne lipider. Efter fjernelse af opløsningsmidlet blev prøverne tørret under normalt atmosfærisk tryk i 24 timer. Prøverne blev derefter demineraliseret i 10 ml 0,5 M HCl i 4 timer under omrøring med en orbitalrystemaskine. Efter demineraliseringen blev prøverne skyllet til neutralitet med type I-vand og derefter opvarmet ved 75 °C i 36 timer i 10-3 M HCl for at opløse kollagenet. Det vandopløselige kollagen blev derefter frysetørret. Kollagenekstraktet blev derefter behandlet med 10 ml 10 ml 10:5:4 kloroform/methanol/vand (v/v/v/v) under sonicering i 1 time for at fjerne eventuelle resterende lipider. Efter centrifugering blev chloroform/methanol-laget fjernet, og methanolen blev fordampet fra vandlaget ved 60 °C i 24 timer. Prøverne blev igen frysetørret og afvejet i tinkapsler med henblik på grundstof- og isotopanalyse. De stabile isotop- og grundstofsammensætninger af kulstof og kvælstof blev bestemt ved hjælp af et Nu Horizon (Nu Instruments, UK) massespektrometer med kontinuerligt isotopforhold ved Trent University, Canada. Stabile kulstof- og kvælstofisotopsammensætninger blev kalibreret i forhold til VPDB- og AIR-skalaerne ved hjælp af USGS40 og USGS41a. Standardusikkerheden blev fastsat til ±0,17‰ for δ13C og ±0,22‰ for δ15N31; yderligere analytiske detaljer findes i supplerende dokument S4. Den statistiske signifikans af forskelle mellem belugas og narhvalers isotopsammensætninger blev vurderet ved hjælp af uparrede t-tests.