Näytteenotto
Analysoimme (i) genominlaajuista DNA-dataa, joka oli poimittu kallon hampaista (MCE1356, Fig. 2b, Supplementary Fig. S1) ja kudosnäytteistä, jotka otettiin kahdeksan belugan ja kahdeksan narvivalaan näytteistä Diskonlahdelta, Länsi-Grönlannista, ja (ii) luun kollageenin stabiilit hiili- ja typpi-isotooppikoostumukset MCE1356:sta ja vertailupaneelista, joka koostui 18 belugan ja 18 narvivalaan kallosta Länsi-Grönlannista (Fig. 1c). Kudosnäytteet (iho) säilytettiin 96-prosenttisessa etanolissa. Näytteet keräsivät Grönlannin luonnonvarainstituutin tutkijat Grönlannin hallituksen inuiittien biologista näytteenottoa koskevan yleisen luvan nojalla, ja ne vietiin Tanskaan CITES-lupien IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 ja 116.2006 nojalla. Näytetiedot on esitetty yksityiskohtaisesti lisätaulukossa S2.
DNA:n uutto, kirjaston valmistelu ja sekvensointi
Kudosnäytteet
Kudosnäytteiden DNA uutettiin Qiagenin Blood and Tissue Kitillä valmistajan protokollan mukaisesti. DNA fragmentoitiin Covaris M220 Focused-ultraäänilaitteella ~350-550 emäsparin (bp) pituisten fragmenttien tuottamiseksi. Kirjastot muodostettiin fragmentoiduista DNA-uutteista Illumina NeoPrep -ohjelmalla NeoPrep Library Prep System -oppaan mukaisesti käyttäen oletusasetuksia. PCR-monistaminen, kvantifiointi ja normalisointi suoritettiin NeoPrep Library Prep System -järjestelmällä. Kirjastot seulottiin kokojakauman osalta Agilent 2100 -bioanalysaattorilla ja yhdistettiin ekvimolaarisissa suhteissa ennen sekvensointia Illumina HiSeq 2500 -laitteella 80 bp:n SE-tekniikalla.
Hammas-/luunäytteet
Kudosnäytteistä poiketen vanhoissa luu- ja hampaidenäytteissä on suhteellisen vähän DNA:ta, mikä edellyttää erilaista uutto- ja kirjastojen rakentamisprotokollaa. Noin 0,5 g luujauhetta viidestä hampaasta ja yksi luun sirpale näytteestä MCE1356 porattiin käsikäyttöisellä Dremelillä. DNA uutettiin hammas-/luujauheesta Kööpenhaminan yliopiston Tanskan luonnonhistoriallisessa museossa sijaitsevassa muinais-DNA:n puhdistuslaboratoriossa käyttäen artikkelissa11 kuvattua uuttopuskuria, johon lisättiin 30 minuutin esilöytövaihe12. Zymo-Spin V -kolonnien (Zymo Research) sijasta uuttopuskuri väkevöitiin Amicon Ultra 30 kDa -sentrifugisuodatinyksiköillä ja väkevöitiin ja puhdistettiin edelleen Qiagen Minelute -putkissa. Puhdistetuista uutteista muodostettiin sitten Illumina-kirjastoja13 kuvatun protokollan mukaisesti. Käytimme qPCR:ää tarkistaaksemme, että kirjaston rakentaminen onnistui, valitaksemme, mitkä kirjastot sekvensoitiin, ja laskeaksemme sopivan määrän PCR-syklejä, joita tarvittiin kunkin kirjaston riittävään monistamiseen aiheuttamatta yliampulointia. Yhteensä neljä kirjastoa monistettiin yksilöllisillä 6 bp:n indekseillä ja seulottiin endogeenisen sisällön varalta Illumina MiSeq -alustalla käyttäen 250 bp:n SE-sekvensointia. Parhaat kirjastot sekvensoitiin sitten uudelleen Illumina HiSeq 2500 -alustalla käyttäen 80 bp SE-tekniikkaa.
Bioinformatiikan analyysi
Kaikki kartoitus- ja DNA-vaurioanalyysit suoritettiin Paleomix pipeline 1.2.1214 -ohjelmalla. Lukemat leikattiin AdapterRemoval 2.2.015:llä käyttäen oletusasetuksia lukuun ottamatta lukujen vähimmäispituutta, joka asetettiin 25 bp:ksi. Lukemat tarkastettiin FastQC:llä ja linjattiin BWA16:lla käyttäen Backtrack-algoritmia ja poistamalla aloitussiemenen pituus käytöstä. Jos lukemat karttuivat useisiin positioihin tai niiden kartoituksen laatupisteet (BWA:n MAPQ-pisteet) olivat alle 30, ne poistettiin SAMtools-ohjelmalla17. Sekvenssiduplikaatit poistettiin Picardin MarkDuplicates-ohjelmalla (saatavilla osoitteesta: http://broadinstitute.github.io/picard), ja lopullinen kohdistus kohdistettiin uudelleen indeleiden ympärille GATK:lla18. Sytosiinin deaminoituminen urasiiliksi näytteessä MCE1356 arvioitiin käyttämällä mapDamage v2.0.619 -ohjelman tulosta. MapDamage-tulokset eivät osoittaneet selvää deaminaatiosignaalia sekvensseissä (Supplementary Fig. S3).
Mitokondrioiden analyysi
MCE1356:n emälinjan määrittämiseksi DNA:n sekvensointilukemat kartoitettiin belugan (GenBank-tietue: KY444734) ja norsunvalaan (GenBank-tietue: NC_005279) mitokondrioiden referenssigenomiin, ja niiden keskimääräistä kattavuutta verrattiin. Vertailupaneelin muodostaneiden kahdeksan beluga- ja kahdeksan narvivalasnäytteen lukemat kartoitettiin vastaaviin mitokondriaalisiin vertailugenomeihin. Muodostimme mitokondriaaliset konsensussekvenssit MCE1356:sta ja 16:sta referenssipaneelin näytteestä, joiden alueet kattoivat yli viisi lukua, käyttämällä BEDtools20-, SAMtools17- ja GATK18-ohjelmia. Loimme kaksi sekvenssikohdistusta käyttäen ClustalW:tä21 käyttäen oletusasetuksia, jotka sisälsivät 16 referenssipaneelin näytettä ja joko sen version MCE1356:sta, joka on kartoitettu belugan mitokondriaaliseen referenssigenomiin, tai sen version MCE1356:sta, joka on kartoitettu narrivalaan mitokondriaaliseen referenssigenomiin. Käytimme näitä kahta linjausta rakentaaksemme mediaaniyhdistävät haplotyyppiverkostot22 käyttäen Popart 1.723 -ohjelmaa (saatavana osoitteesta: http://popart.otago.ac.nz) ja jättäen pois kaikki kohdat, joissa oli indeleitä tai puuttuvia tietoja. Tämän jälkeen sekä kartoituksen jälkeistä kattavuutta että kahta haplotyyppiverkostoa käytettiin MCE1356′:n emälinjan lajin määrittämiseen.
Ydin-DNA-analyysit
Kartoitimme kaikista näytteistä saadut DNA:n sekvensointilukemat tappajavalaan (Orcinus orca) referenssigenomiin (GCA_000331955.2). Laadukas valaan genomi julkaistiin hiljattain24 , mutta kartoitus jompaankumpaan mahdolliseen vanhempiin lajeihin saattaisi aiheuttaa vääristymiä analyyseihin. Siksi kartoitimme lukemat tappajavalaan genomiin, koska se on hyvin koottu ja koska tappajavalaat ovat edelleen suhteellisen läheistä sukua belugoille ja narvivalaaneille (divergenssiaika 12 MYA)25 , mutta silti riittävän kaukana toisistaan, jotta analyysejämme hankaloittavan introgression riski olisi mahdollisimman pieni.
Kaikessa muussa suodatuksessa käytimme ANGSD v0.92326 -ohjelmistopakettia, joka käyttää genotyyppien todennäköisyyksiä kutsuttujen genotyyppien sijaan, mikä on erityisen hyödyllistä silloin, kun analyysissä käytetään matalan kattavuuden omaavia NGS-tietoja. Käytimme ANGSD:ssä toteutettua SAMtools17 -menetelmää genotyyppien todennäköisyyksien arvioimiseksi ja päättelimme pää- ja sivu-alleelit suoraan genotyyppien todennäköisyyksistä käyttäen maksimiluotettavuuslähestymistapaa, joka on kuvattu teoksessa27.
Kartoitetun datan visualisoimiseksi piirsimme lukusyvyysjakauman kaikilta yksilöiltä sulkemalla pois paikat, joissa oli enemmän kuin kaksi alleelia, ja paikat, joiden Phred-pistemäärä oli alle 25. Visualisoimalla kaikkien yksilöiden tiedot yhdessä pystyimme arvioimaan keskimääräisen lukusyvyyden 4,14-kertaiseksi ja tunnistamaan kohteet, joissa lukusyvyys oli koholla. Tällaiset kohteet olivat todennäköisemmin peräisin paralogeista ja genomin toistuvista alueista. Tietokokonaisuutta tarkasteltiin visuaalisesti ja suodatettiin edelleen, jolloin jätettiin pois kaikki kohteet, joiden lukusyvyys oli yli yhdeksän (6,9 % kohteista). ANGSD:ssä SNP:t kutsuttiin niiden alleelifrekvenssien perusteella. Vähäinen alleelifrekvenssi (minor allele frequency, MAF) arvioitiin genotyyppien todennäköisyyksistä, ja todennäköisyyssuhdetestiä käytettiin sen määrittämiseksi, poikkesiko MAF nollasta. Jos likelihood ratio -testin p-arvo oli <1e-4, kohdetta pidettiin polymorfisena ja se säilytettiin aineistossa. Tämän SNP-suodattimen soveltaminen merkitsi sitä, että yhtään kohdetta, jossa oli vähemmän kuin neljä lukua, ei säilytetty, koska vähemmän lukuja sisältävien kohteiden p-arvot eivät voi olla näin alhaiset. Näitä suodattimia sovellettiin aineistoon, joka sisälsi kaikki 17 näytettä, ja aineistoon, jossa ei ollut MCE1356:ta.
Määrittääksemme, oliko MCE1356 belugan ja norsunvalaan hybridi, suodatimme edelleen 17 yksilön aineiston, jättäen pois sivustot, joissa ei ollut lukemia MCE1356:ssa. Tässä vaiheessa MCE1356:n keskimääräinen lukusyvyys oli vain 1,15-kertainen, joten varmistaaksemme, ettemme analysoi monikopiointilokuksia, jätimme pois paikat, jotka kattavat useamman kuin yhden lukeman MCE1356:ssa.
Tarkoituksemme oli verrata MCE1356:sta löydettyjä alleeleja referenssipaneelissa esiintyviin alleeleihin, joten arvioimme todennäköisyyden, jolla MCE1356:sta löydetty alleeli olisi osoitettu väärälle vanhemmalle lajille, kun otetaan huomioon referenssipaneelissa esiintyvien lukusyvyydeltään yksilöllisten lukusyvyydeltään minimaalisten lukusyvyydeltään yksilöllisten lukujen erilaiset minimiarvot ja alleelien esiintymistiheydet. Yksilöllinen lukusyvyys määriteltiin tietyn paikan kattavien lukujen määräksi, kun kaikki lukemat olivat peräisin yksilöltä. Tämä todennäköisyys P laskettiin yhtälön 1 mukaisesti:
Jossa f(ps) on alleelifrekvenssi vanhempien lajeissa ja urd(psp) lajikohtainen ainutkertainen lukusyvyys referenssipanelissa.
Perhelajin alleelifrekvenssi, joka antaa suurimman todennäköisyyden f(ps-max), voidaan kuvata yhtälön 2 mukaisesti:
Sisällytetään yhtälö 2 yhtälöön 1, laskettiin suurin todennäköisyys P(max) sille, että MCE1356:sta löydetty alleeli osoitetaan väärälle vanhempien lajille, yhtälön 3 mukaisesti:
Tulokset osoittivat, että ainutlaatuisella lukusyvyydellä kaksi, kolme ja neljä kussakin kantalajissa, suurin todennäköisyys määrittää MCE1356:sta löydetty alleeli väärään kantalajiin oli 0.148, 0,105 ja 0,082. Näitä enimmäisarvoja ei pitäisi tulkita virhemääriksi, sillä ne saataisiin vain, jos kaikki paikat olisivat vaihtelevia vanhempien lajien sisällä ja kaikkien paikkojen MAF olisi täsmälleen (1 – (urd/urd + 1)). Lisäksi olettaen, että MAF-jakaumat belugoilla ja norsunvalailla olisivat samanlaiset, MCE1356:ssa havaittujen alleelien virheellinen osoittaminen vaikuttaisi molempiin lajeihin yhtä paljon ja vaikuttaisi näin ollen vain vähän hybridisaatiota koskeviin johtopäätöksiin. Yksilöllisen lukusyvyyden käyttämisestä yhtälöissä 2 ja 3 oli se hyöty, että se yhdisti yksilöiden määrän ja lukujen määrän arvioitaessa todennäköisyyttä, jolla MCE1356:sta löydetty alleeli osoitetaan väärälle vanhempien lajille. Tämän jälkeen lukusyvyysjakauma, MAF (kiinteällä pää- ja sivu-alleelilla) ja niiden yksilöiden lukumäärä, joilla oli tietoja kussakin paikassa, laskettiin erikseen kullekin vanhempaislajille.
Konstruoitiin kolme tietokokonaisuutta, jotka sisälsivät MCE1356:n lisäksi vanhempaislajipaneeleita, joiden minimiyksilölliset lukusyvyydet olivat kaksi, kolme ja neljä. Näissä kolmessa tietokokonaisuudessa säilytettyjen paikkojen määrä oli 61 105, 11 764 ja 360. Kompromissina sivustojen lukumäärän maksimoimisen ja MCE1356:ssa havaittujen alleelien virheellisen osoittamisen minimoimisen välillä päätimme käyttää tietokokonaisuutta, jossa MCE1356:n lukusyvyys on yksi ja joka sisältää vähintään kolmen lukusyvyyden lukusyvyyden kussakin vanhempien lajissa. Tämä tietokokonaisuus sisälsi 11 764 paikkaa, joita käytettiin myöhemmissä analyyseissä.
Yhteenvetotilastot tehtiin tietokokonaisuudelle, jossa oli 17 yksilöä, mukaan lukien niiden paikkojen lukumäärä, jotka olivat (i) kiinteitä eri alleeleille beluga- ja narvivalaslajipaneeleissa; (ii) polymorfisia belugoissa, mutta eivät narvivalaslajeissa; (iii) polymorfisia narvivalaslajipaneeleissa, mutta eivät belugoissa; (iv) polymorfisia sekä belugoissa että narvivalaslajeissa. Kohteet, joiden arvioidaan olevan kiinteitä kahden kantalajipaneelin välillä, rikastuvat markkereiden osalta, jotka eroavat suuresti toisistaan eli joilla on suuria alleelifrekvenssieroja näiden kahden kantalajin välillä. Vaikka nämä markkerit eivät välttämättä olekaan kiinteitä eroja kahden vanhemman lajin välillä, ne ovat silti erittäin informatiivisia MCE1356:n syntyperän suhteen.
Käytimme genotyyppien todennäköisyyksiä aineistosta ilman MCE1356:ta varmistaaksemme, että referenssipaneelissa olevat belugat ja narvivalaat eivät itse olleet hiljattain sekoittuneita yksilöitä. Arvioimme niiden yksilölliset sekoittumiskertoimet, kun määrittelimme kaksi populaatiota (K = 2) käyttäen NgsAdmix28 -ohjelmaa. Suoritettiin sata ajoa, ja sekoittumiskertoimien keskiarvoa ja keskihajontaa käytettiin myöhempää tulkintaa varten. Varmistaaksemme, että suodattimet eivät olleet paljastaneet aiemmin piilossa olleita sekoittuneita geneettisiä profiileja referenssipaneelissa, tämä analyysi suoritettiin sekä ennen että jälkeen uniikkien lukusyvyyssuodattimien soveltamisen.
Analysoimme MCE1356:n sekoittumisprosentit käyttämällä fastNGSadmix29 -ohjelmaa soveltaen 1000 bootstrap-analyysiä. fastNGSadmix-ohjelma käyttää referenssipopulaatioiden/lajien alleelifrekvenssejä ja yksittäisen yksilön genotyyppien todennäköisyyttä estimoidakseen sekoittumisprosentit. Ohjelmisto on osoittautunut kestäväksi NGS-datalla, jonka kattavuus on niinkin alhainen kuin 0,00015×29, ja se oli siksi ihanteellinen tutkimukseemme.
Arvioimme edelleen MCE1356:n hybridistatusta tutkimalla eri alleeleille kiinnitettyjä paikkoja belugapaneelissa ja narrivalaspaneelissa (9 178 paikkaa) ja vertaamalla havaittua belugaspesifisen alleelin ja narrivalaspesifisen alleelin kanssa täsmäävän lukeman osuutta MCE1356:ssa odotettavissa oleviin osuuksiin erilaisissa hybridisaatioskenaarioissa. Määrittääksemme, kuinka hyvin seitsemän erilaista hybridisaatioskenaariota sopi havaittuihin tietoihin, laskimme Pearsonin Khiin neliö -tilaston sopivuuden hyvyydestä. Testistatistiikka lasketaan yhtälön 4 mukaisesti:
joissa Oi ja Ei ovat vanhempien lajista i peräisin olevien alleelien havaitut ja odotetut määrät. Nollahypoteesissa, jossa valittu skenaario vastaa hyvin havaittuja tietoja, testistatistiikka T noudattaa keskeistä χ2-jakaumaa. Näin ollen skenaario, joka vastaa parhaiten havaittuja tietoja, johtaisi pienimpään testistatistiikkaan.
Seitsemän eri hybridisaatioskenaarion tarkemmaksi tutkimiseksi laskimme havaittujen alleelien todennäköisyyden paikoissa, jotka oli kiinnitetty eri alleeleille vanhempien populaatioissa. Tarkemmin sanottuna laskimme hybridissä MCE1356 havaittujen alleelien todennäköisyyden vanhempien lajeista peräisin olevien alleelien periytymisen binomiaalimallilla olettaen lisäksi, että markkerien välillä vallitsee riippumattomuus. Haluamme kuitenkin huomauttaa, että riippumattomuusolettaman rikkominen johtaa edelleen harhattomiin estimaatteihin. Jos oletetaan, että hybridi MCE1356 koostuu osuudesta b, joka on peräisin belugan polveutumisesta, ja osuudesta (1-b), joka on peräisin narvaalin polveutumisesta, voimme kirjoittaa b:n todennäköisyyden k:n riippumattoman paikan todennäköisyyden tulona, kun otetaan huomioon niiden lukemien lukumäärä nib, jotka täsmäävät belugan alleeliin paikassa i, ja niiden lukemien lukumäärä nin, jotka täsmäävät narvaalin polveutumiseen, kuten yhtälössä 5 on esitetty.
Laskimme b:n log-likelihoodin koko sen alueella, 0:sta 1:een, ja vertasimme log-likelihoodeja seitsemässä eri hybridisaatioskenaariossa.
Sukupuolen määrittäminen
Myös MCE1356:n sekä beluga- ja narvivalasvertailupaneelien yksilöiden sukupuolen määrittämiseksi tutkittiin X-kromosomin ja autosomaalisen kattavuuden suhdetta. Tämä tehtiin vertaamalla X-kromosomista oletettavasti peräisin olevien telineiden kattavuutta autosomeista oletettavasti peräisin olevien telineiden kattavuuteen. Määritimme, mitkä telineet ovat peräisin X-kromosomista, kohdistamalla tappajavalaan genomi lehmän (Bos taurus) X-kromosomiin (CM008168.2) SatsumaSynteny30-ohjelmalla. Lisäksi poistaaksemme vääristymät, jotka voivat johtua X- ja Y-kromosomien joidenkin alueiden homologiasta, kohdistimme tuloksena saadut oletetut X-kromosomin telineet ihmisen Y-kromosomiin (NC_000024.10) ja poistimme nämä telineet jatkoanalyysistä. Tämän jälkeen laskimme kattavuuden jäljelle jääneille X-telineille ja neljälle suurimmalle telineelle, jotka eivät olleet linjassa X-kromosomin kanssa (eli autosomaalisille telineille), käyttämällä SAMtools-ohjelman syvyysfunktiota17 . Kompensoidaksemme peittävyyden vaihtelua koko genomissa valitsimme satunnaisesti 10 000 000 paikkaa, laskimme näiden paikkojen keskimääräisen peittävyyden ja toistimme tämän vaiheen 100 kertaa. Kullekin yksilölle laskettiin 5 %:n ja 95 %:n luottamusvälit sekä ensimmäinen ja kolmas kvartiili, joita käytettiin myöhempää tulkintaa varten.
Vakaan hiilen ja typen isotooppianalyysi
Jauhemaiset luunäytteet (~100 mg) MCE1356:n, 18 belugan ja 18 narvivalaan kalloista käsiteltiin 10 ml:lla kloroformia/metanolia (v/v) (2:1) ultraäänihajotuksen alaisena 1 tunnin ajan lipidien poistamiseksi. Liuottimen poistamisen jälkeen näytteet kuivattiin normaalissa ilmanpaineessa 24 tunnin ajan. Tämän jälkeen näytteitä demineralisoitiin 10 ml:ssa 0,5 M HCl:ää 4 tunnin ajan, kun niitä sekoitettiin kiertorattaalla. Demineralisoinnin jälkeen näytteet huuhdeltiin neutraaliksi tyypin I vedellä, minkä jälkeen niitä lämmitettiin 75 °C:ssa 36 tuntia 10-3 M HCl:ssä kollageenin liuottamiseksi. Vesiliukoinen kollageeni pakastettiin sen jälkeen. Tämän jälkeen kollageeniuutetta käsiteltiin 10 ml:lla 10:5:4 kloroformia/metanolia/ vettä (v/v/v) sonikoimalla 1 tunnin ajan jäljellä olevien lipidien poistamiseksi. Sentrifugoinnin jälkeen kloroformi/metanolikerros poistettiin ja metanoli haihdutettiin pois vesikerroksesta 60 °C:ssa 24 tunnin ajan. Näytteet pakastettiin jälleen ja punnittiin tinakapseleihin alkuaine- ja isotooppianalyysiä varten. Hiilen ja typen stabiilit isotooppi- ja alkuainekoostumukset määritettiin jatkuvatoimisella Nu Horizon (Nu Instruments, Yhdistynyt kuningaskunta) -isotooppisuhdemassaspektrometrillä Trentin yliopistossa Kanadassa. Hiilen ja typen stabiilit isotooppikoostumukset kalibroitiin suhteessa VPDB- ja AIR-asteikkoihin käyttäen USGS40- ja USGS41a-mittareita. Standardiepävarmuudeksi määritettiin ±0,17‰ δ13C:n osalta ja ±0,22‰ δ15N31:n osalta; analyyttiset lisätiedot on esitetty lisäasiakirjassa S4. Belugan ja norsunvalaan isotooppikoostumusten välisten erojen tilastollinen merkitsevyys arvioitiin käyttämällä parittomia t-testejä.