Eșantionare
Am analizat (i) datele ADN la nivel genomic extrase din dinții craniului (MCE1356, Fig. 2b, Fig. Suplimentară S1) și din eșantioane de țesut de la opt beluga și opt narvali prelevate în Golful Disko, în vestul Groenlandei, și (ii) compozițiile izotopice stabile de carbon și azot ale colagenului osos din MCE1356 și dintr-un panel de referință de 18 cranii de beluga și 18 cranii de narval din vestul Groenlandei (Fig. 1c). Probele de țesut (piele) au fost depozitate în etanol 96%. Eșantioanele au fost recoltate de oamenii de știință de la Institutul de Resurse Naturale din Groenlanda în temeiul permisului general de prelevare de probe biologice de la Inuit de la guvernul Groenlandei și au fost exportate în Danemarca în baza permiselor CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 și 116.2006. Informațiile despre eșantioane sunt detaliate în tabelul suplimentar S2.
Extragerea ADN-ului, pregătirea bibliotecilor și secvențierea
Eșantioane de țesut
ADN-ul din eșantioanele de țesut a fost extras cu ajutorul kitului Qiagen Blood and Tissue Kit în conformitate cu protocolul producătorului. ADN-ul a fost fragmentat cu ajutorul aparatului Covaris M220 Focused-ultrasonicator pentru a crea fragmente cu o lungime de ~350-550 de perechi de baze (bp). Bibliotecile au fost construite din extractele de ADN fragmentate folosind Illumina NeoPrep, în conformitate cu Ghidul sistemului NeoPrep Library Prep, aplicând setările implicite. Amplificarea PCR, cuantificarea și normalizarea au fost toate efectuate de NeoPrep Library Prep System. Bibliotecile au fost verificate pentru distribuția dimensională pe un bioanalizator Agilent 2100 și au fost grupate în proporții echimolare înainte de secvențierea pe un Illumina HiSeq 2500 cu tehnologie SE de 80 bp.
Eșantioane de dinți/oase
Spre deosebire de probele de țesut, probele de oase și dinți bătrâni au concentrații relativ scăzute de ADN, ceea ce necesită protocoale diferite de extracție și construire a bibliotecilor. Aproximativ 0,5 g de pulbere de os de la cinci dinți și un ciob de os de la specimenul MCE1356 au fost forate cu ajutorul unui Dremel portabil. ADN-ul a fost extras din pulberea de dinți/oase într-un laborator dedicat curățării ADN-ului antic din cadrul Muzeului de Istorie Naturală din Danemarca, Universitatea din Copenhaga, utilizând tamponul de extracție descris în11 cu adăugarea unei etape de predigerare de 30 de minute12. În loc să se utilizeze coloanele Zymo-Spin V (Zymo Research), tamponul de extracție a fost concentrat cu ajutorul unităților de filtrare centrifugale Amicon Ultra 30 kDa și concentrat și curățat în continuare cu ajutorul tuburilor Qiagen Minelute. Extractele purificate au fost apoi încorporate în bibliotecile Illumina în conformitate cu protocolul descris de13. Am utilizat qPCR pentru a verifica dacă biblioteca a fost construită cu succes, pentru a selecta bibliotecile care vor fi secvențiate și pentru a calcula numărul adecvat de cicluri PCR necesare pentru a amplifica suficient fiecare bibliotecă fără a provoca supraamplificare. În total, patru biblioteci au fost amplificate cu indici unici de 6 bp și au fost analizate pentru conținut endogen pe platforma Illumina MiSeq, utilizând secvențierea SE de 250 bp. Cele mai bune biblioteci au fost apoi resecvențiate pe platforma Illumina HiSeq 2500 utilizând tehnologia 80 bp SE.
Analiză bioinformatică
Toate analizele de cartografiere și de deteriorare a ADN-ului au fost efectuate în cadrul pipeline-ului Paleomix 1.2.1214. Lecturile au fost tăiate cu AdapterRemoval 2.2.015 utilizând setările implicite, cu excepția lungimii minime de citire, care a fost setată la 25 bp. Citirile au fost inspectate cu FastQC și aliniate cu BWA16 aplicând algoritmul Backtrack și dezactivând lungimea de pornire a semințelor. În cazul în care citirile se corelau cu mai multe poziții sau aveau scoruri de calitate a cartografierii (scor MAPQ din BWA) mai mici de 30, acestea au fost eliminate cu ajutorul SAMtools17. Dublurile de secvențe au fost eliminate folosind MarkDuplicates din Picard (disponibil la: http://broadinstitute.github.io/picard), iar alinierea finală a fost realiniată în jurul indelilor folosind GATK18. Deaminarea citozinei în uracil în specimenul MCE1356 a fost evaluată utilizând rezultatul din mapDamage v2.0.619. Rezultatele mapDamage nu au arătat niciun semnal clar de dezaminare în secvențe (Fig. Suplimentară S3).
Analiză mitocondrială
Pentru a determina descendența maternă a MCE1356, citirile de secvențiere a ADN-ului au fost mapate pe genomul mitocondrial de referință al belugii (accesare GenBank: KY444734) și al narvalului (accesare GenBank: NC_005279) și a fost comparată acoperirea medie. Lecturile din cele opt eșantioane de beluga și opt eșantioane de narval care au alcătuit panoul de referință au fost puse în corespondență cu genomurile lor mitocondriale de referință respective. Am construit secvențe consensuale mitocondriale ale MCE1356 și ale celor 16 eșantioane din panoul de referință cu regiuni acoperite de mai mult de cinci citiri, utilizând BEDtools20, SAMtools17 și GATK18. Am creat două alinieri de secvențe cu ajutorul ClustalW21, aplicând setările implicite, care au inclus cele 16 eșantioane din panoul de referință și fie versiunea MCE1356 cartografiată la genomul mitocondrial de referință al belugii, fie versiunea MCE1356 cartografiată la genomul mitocondrial de referință al narvalului. Am folosit cele două alinieri pentru a construi rețele de haplotipuri cu joncțiune mediană22 folosind Popart 1.723 (disponibil la: http://popart.otago.ac.nz), excluzând orice situs cu indels sau date lipsă. Ulterior, atât acoperirea post-mapping, cât și cele două rețele de haplotipuri au fost utilizate pentru a determina specia liniei materne a MCE1356′s.
Analize ADN nuclear
Am cartografiat citirile de secvențiere a ADN-ului din toate probele la genomul de referință al balenei ucigașe (Orcinus orca) (GCA_000331955.2). Un genom de înaltă calitate al balenei beluga a fost publicat recent24, dar cartografierea la una dintre cele două specii parentale potențiale ar putea crea distorsiuni în analize. Prin urmare, am cartografiat citirile la genomul balenei ucigașe, deoarece acesta este bine asamblat și balenele ucigașe sunt încă relativ strâns înrudite cu beluga și narvalii (timp de divergență de 12 MYA)25, dar suficient de îndepărtate pentru a minimiza riscul de introgresie care ar complica analizele noastre.
Pentru toate filtrările ulterioare am utilizat ANGSD v0.92326, un pachet software care utilizează probabilitățile genotipurilor în loc de genotipurile numite, care este deosebit de util atunci când se analizează date NGS cu acoperire redusă. Am utilizat metoda SAMtools17 implementată în ANGSD pentru a estima probabilitățile genotipurilor și am dedus alelele majore și minore direct din probabilitățile genotipurilor utilizând o abordare de maximă probabilitate, așa cum este descrisă în27.
Pentru a vizualiza datele cartografiate, am trasat distribuția adâncimii citite de la toți indivizii, excluzând siturile cu mai mult de două alele și siturile cu un scor Phred sub 25. Vizualizarea datelor de la toți indivizii combinați ne-a permis să estimăm adâncimea medie de citire de 4,14x și să identificăm site-urile cu adâncime de citire ridicată. Astfel de situri au fost mai susceptibile de a proveni din paralogi și regiuni repetitive ale genomului. Setul de date a fost inspectat vizual și filtrat în continuare, excluzând toate site-urile cu o adâncime de citire mai mare de nouă (6,9 % din site-uri). În ANGSD, SNP-urile au fost numite pe baza frecvențelor alelelor lor. Frecvența alelei minore (MAF) a fost estimată din probabilitățile genotipurilor și a fost utilizat un test al raportului de verosimilitate pentru a determina dacă MAF era diferită de zero. În cazul în care valoarea p din testul raportului de verosimilitate a fost <1e-4, situl a fost considerat polimorf și a fost păstrat în setul de date. Aplicarea acestui filtru SNP a însemnat că niciun sit cu mai puțin de patru citiri nu a fost reținut, deoarece siturile acoperite de mai puține citiri nu puteau obține valori p atât de scăzute. Aceste filtre au fost aplicate unui set de date care includea toate cele 17 eșantioane, precum și unui set de date fără MCE1356.
Pentru a determina dacă MCE1356 era un hibrid beluga/narval, am filtrat în continuare setul de date de 17 indivizi, excluzând site-urile fără lecturi în MCE1356. În acest moment, adâncimea medie de citire a MCE1356 era de numai 1,15x, astfel încât, pentru a ne asigura că nu analizăm loci multicopie, am exclus siturile acoperite de mai mult de o citire în MCE1356.
Scopul nostru a fost să comparăm alelele găsite în MCE1356 cu alelele din panoul de referință, astfel încât am estimat probabilitatea de a atribui alela găsită în MCE1356 la specia parentală greșită, având în vedere diferitele adâncimi de citire și frecvențe alele unice minime ale panoului de referință. Adâncimea unică de citire a fost definită ca fiind numărul de citiri care acoperă un anumit sit, în cazul în care toate citirile provin de la un singur individ. Această probabilitate P a fost calculată conform ecuației 1:
Unde f(ps) este frecvența alelei în specia parentală, iar urd(psp) este adâncimea de citire unică specifică speciei în panoul de referință.
Frecvența alelei din specia parentală care dă cea mai mare probabilitate f(ps-max) ar putea fi descrisă ca în ecuația 2:
Prin introducerea Ecuației 2 în ecuația 1, probabilitatea maximă P(max) de a atribui alela găsită la MCE1356 la specia parentală greșită a fost calculată ca în ecuația 3:
Rezultatele au arătat că, cu o adâncime de citire unică de două, trei și patru în fiecare specie parentală, probabilitatea maximă de a atribui alela găsită în MCE1356 la specia parentală greșită a fost 0.148, 0,105 și, respectiv, 0,082. Aceste valori maxime nu trebuie interpretate ca fiind rate de eroare, deoarece acestea ar fi obținute numai dacă toate siturile ar fi variabile în cadrul speciei parentale și toate siturile ar avea un MAF de exact (1 – (urd/urd + 1)). În plus, presupunând că distribuțiile MAF la beluga și narvalul sunt similare, atribuirea eronată a alelelor găsite în MCE1356 ar afecta ambele specii în mod egal și, prin urmare, ar avea o influență minimă asupra inferențelor de hibridare. Un avantaj al utilizării profunzimii unice de citire în ecuațiile 2 și 3 a fost acela că a combinat numărul de indivizi și numărul de citiri în estimarea probabilității de atribuire a alelei găsite în MCE1356 la specia parentală greșită. Ulterior, distribuția adâncimii de citire, MAF (cu alela majoră și minoră fixă) și numărul de indivizi cu date în fiecare sit au fost calculate separat pentru fiecare specie parentală.
Au fost construite trei seturi de date, care, pe lângă MCE1356, au inclus panouri de specii parentale cu adâncimi de citire unice minime de două, trei și patru. Numărul de situri reținute în cele trei seturi de date a fost de 61.105, 11.764 și, respectiv, 360. Ca un compromis între maximizarea numărului de situri și minimizarea atribuirii greșite a alelelor găsite în MCE1356, am ales să folosim setul de date cu o singură citire în MCE1356 și adâncimi minime unice de citire de trei în fiecare specie parentală. Acest set de date a inclus 11.764 de situri, care au fost utilizate în analizele ulterioare.
Statisticile sumare au fost efectuate pe setul de date cu 17 indivizi, inclusiv numărul de situri care au fost (i) fixate pentru diferite alele în panourile de specii de beluga și narval; (ii) polimorfice la beluga, dar nu și la narval; (iii) polimorfice la narval, dar nu și la beluga; (iv) polimorfice atât la beluga, cât și la narval. Locurile care se estimează a fi fixe între cele două paneluri de specii parentale vor fi îmbogățite pentru markerii care sunt foarte diferențiați, adică au diferențe mari de frecvență alelelor, între cele două specii parentale. Acești markeri, deși nu sunt neapărat diferențe fixe între cele două specii parentale, sunt totuși foarte informativi pentru strămoșii din MCE1356.
Am folosit probabilitățile genotipurilor din setul de date fără MCE1356 pentru a verifica dacă belugii și narvalii din panoul de referință nu au fost ei înșiși indivizi amestecați recent. Am estimat coeficienții lor individuali de amestec în timp ce am specificat două populații (K = 2) folosind NgsAdmix28. Au fost efectuate o sută de rulări, iar media și deviația standard a coeficienților de amestec au fost utilizate pentru interpretarea ulterioară. Pentru a confirma faptul că filtrele nu au dezvăluit profiluri genetice de amestec ascunse anterior în panoul de referință, această analiză a fost efectuată atât înainte, cât și după aplicarea filtrelor de adâncime de citire unice.
Am analizat proporțiile de amestec ale MCE1356 utilizând fastNGSadmix29, aplicând 1 000 de bootstrap-uri. fastNGSadmix utilizează frecvențele alelelor populațiilor/speciilor de referință și probabilitățile genotipice ale unui singur individ pentru a estima proporțiile de amestec ale acestuia. Software-ul s-a dovedit robust cu date NGS cu o acoperire de până la 0,00015×29 și, prin urmare, a fost ideal pentru studiul nostru.
Am estimat în continuare statutul hibrid al MCE1356 prin investigarea siturilor fixate pentru diferite alele în panoul de beluga și panoul de narval (9 178 de situri) și prin compararea proporției observate de citiri care corespund alelei specifice belugii și alelei specifice narvalului din MCE1356 cu proporțiile așteptate în diferite scenarii de hibridare. Pentru a determina în ce măsură cele șapte scenarii de hibridare diferite se potrivesc cu datele observate, am calculat o statistică Pearson’s Chi-square goodness of fit. Statistica testului este calculată ca în Ecuația 4:
unde Oi și Ei sunt numerele observate și, respectiv, așteptate de alele derivate din specia parentală i. În ipoteza nulă în care scenariul ales corespunde bine datelor observate, statistica de test T urmează o distribuție centrală χ2. Astfel, scenariul care corespunde cel mai bine cu datele observate ar conduce la cea mai mică statistică de testare.
Pentru a investiga în continuare cele șapte scenarii diferite de hibridare, am calculat probabilitatea alelelor observate în situri care au fost fixate pentru diferite alele în populațiile parentale. Mai exact, am calculat probabilitatea alelelor observate în hibridul MCE1356, în cadrul unui model binomial pentru moștenirea alelelor de la speciile parentale, presupunând în continuare independența între markeri. Cu toate acestea, am dori să observăm că încălcarea ipotezei de independență conduce în continuare la estimări nepărtinitoare. Presupunând că hibridul MCE1356 este compus dintr-o proporție b de strămoși beluga și (1-b) de strămoși narval, putem scrie probabilitatea lui b ca produs al probabilității la k situsuri independente, având în vedere numărul de citiri nib care corespund alelei beluga la situl i și nin care corespund strămoșilor narval, ca în ecuația 5.
Am calculat probabilitatea logaritmică a lui b pe tot intervalul său, de la 0 la 1, și am comparat probabilitățile logaritmice pentru cele șapte scenarii diferite de hibridare.
Determinarea sexului
Pentru a determina sexul lui MCE1356 și al indivizilor din panourile de referință beluga și narval, am investigat raportul de acoperire a cromozomului X față de autosomal. Acest lucru a fost realizat prin compararea acoperirii între scheletele care provin probabil din cromozomul X, cu cea a scheletelor care provin probabil din autosomi. Am determinat schelele care provin din cromozomul X prin alinierea genomului balenei ucigașe la cromozomul X al vacii (Bos taurus) (CM008168.2) cu SatsumaSynteny30. În plus, pentru a elimina prejudecățile care ar putea apărea din cauza omologiei dintre unele regiuni ale cromozomilor X și Y, am aliniat în continuare scheletele presupuse ale cromozomului X rezultate la cromozomul Y uman (NC_000024.10) și am eliminat aceste schele de la analiza ulterioară. Am calculat apoi acoperirea pentru schelele X rămase și pentru cele patru schele mai mari care nu se aliniază la cromozomul X (adică schelele autosomale), utilizând funcția de adâncime din SAMtools17. Pentru a compensa orice variație a acoperirii în cadrul genomului, am selectat la întâmplare 10.000.000 de situri, am calculat acoperirea medie pentru aceste situri și am repetat acest pas de 100 de ori. Pentru fiecare individ, s-au calculat intervalele de încredere de 5% și 95%, precum și primul și al treilea cuartil, care au fost utilizate pentru interpretarea ulterioară.
Analiză a izotopilor stabili de carbon și azot
Eșantioanele de os pulverizat (~100 mg) din craniile MCE1356, 18 beluga și 18 narval au fost tratate cu 10 ml de cloroform/metanol 2:1 (v/v) sub sonicare timp de 1 h pentru a elimina lipidele. După îndepărtarea solventului, probele au fost uscate la presiune atmosferică normală timp de 24 h. Probele au fost apoi demineralizate în 10 ml de HCl 0,5 M timp de 4 h, fiind agitate cu un agitator orbital. După demineralizare, probele au fost clătite până la neutralitate cu apă de tip I, iar apoi au fost încălzite la 75 °C timp de 36 h în HCl 10-3 M pentru a solubiliza colagenul. Colagenul solubil în apă a fost apoi liofilizat. Extractul de colagen a fost apoi tratat cu 10 ml de 10:5:4 cloroform/metanol/apă (v/v/v/v) sub sonicare timp de 1 h pentru a elimina orice lipide reziduale. După centrifugare, stratul de cloroform/metanol a fost îndepărtat, iar metanolul a fost evaporat din stratul de apă la 60 °C timp de 24 h. Probele au fost din nou liofilizate și cântărite în capsule de staniu pentru analiza elementară și izotopică. Compozițiile izotopice stabile și elementare ale carbonului și azotului au fost determinate cu ajutorul unui spectrometru de masă de raport izotopic în flux continuu Nu Horizon (Nu Instruments, Regatul Unit) la Universitatea Trent, Canada. Compozițiile izotopice stabile ale carbonului și azotului au fost calibrate în raport cu scările VPDB și AIR folosind USGS40 și USGS41a. Incertitudinea standard a fost determinată ca fiind de ±0,17‰ pentru δ13C și ±0,22‰ pentru δ15N31; detalii analitice suplimentare sunt furnizate în documentul suplimentar S4. Semnificația statistică a diferențelor dintre compozițiile izotopice ale belugilor și narvalilor a fost evaluată cu ajutorul testelor t nepereche.
.