Campionamento
Abbiamo analizzato (i) dati sul DNA genomico estratti da denti del cranio (MCE1356, Fig. 2b, Fig. S1 supplementare) e da campioni di tessuto di otto beluga e otto narvali campionati nella baia di Disko, Groenlandia occidentale, e (ii) composizioni isotopiche stabili di carbonio e azoto del collagene osseo da MCE1356 e un pannello di riferimento di 18 beluga e 18 teschi narvalo dalla Groenlandia occidentale (Fig. 1c). I campioni di tessuto (pelle) sono stati conservati in etanolo al 96%. I campioni sono stati raccolti da scienziati del Greenland Institute of Natural Resources sotto il permesso generale per il campionamento biologico degli Inuit del governo della Groenlandia ed esportati in Danimarca sotto il permesso CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 e 116.2006. Le informazioni sui campioni sono dettagliate nella tabella supplementare S2.
Estrazione del DNA, preparazione della libreria e sequenziamento
Campioni di tessuto
Il DNA dai campioni di tessuto è stato estratto utilizzando il kit Qiagen Blood and Tissue seguendo il protocollo del produttore. Il DNA è stato frammentato utilizzando il Covaris M220 Focused-ultrasonicator per creare ~ 350-550 paia di basi (bp) lunghezze del frammento. Librerie sono stati costruiti dagli estratti di DNA frammentato utilizzando Illumina NeoPrep seguendo la Guida NeoPrep Library Prep sistema applicando le impostazioni predefinite. Amplificazione PCR, quantificazione e normalizzazione sono stati tutti effettuati dal NeoPrep Library Prep System. Le librerie sono stati controllati per la distribuzione delle dimensioni su un Agilent 2100 Bioanalyzer e raggruppati in rapporti equimolari prima di sequenziamento su un Illumina HiSeq 2500 con 80 bp SE tecnologia.
Dente / osso campioni
A differenza dei campioni di tessuto, osso vecchio e campioni di denti hanno relativamente basse concentrazioni di DNA, che richiede l’estrazione diversa e protocolli biblioteca costruire. Circa 0,5 g di polvere d’osso da cinque denti e un frammento d’osso dal campione MCE1356 sono stati forati usando un Dremel manuale. Il DNA è stato estratto dalla polvere di dente/ossa in un laboratorio dedicato alla pulizia del DNA antico presso il Museo di Storia Naturale della Danimarca, Università di Copenhagen, utilizzando il buffer di estrazione descritto in11 con l’aggiunta di una fase di predigestione di 30 minuti12. Invece di utilizzare le colonne Zymo-Spin V (Zymo Research), il buffer di estrazione è stato concentrato utilizzando le unità di filtraggio centrifugo Amicon Ultra 30 kDa e ulteriormente concentrato e pulito utilizzando le provette Qiagen Minelute. Gli estratti purificati sono stati poi costruiti in Illumina librerie seguendo il protocollo descritto da13. Abbiamo usato qPCR per verificare che la biblioteca costruire è stato successo, per selezionare quali librerie di sequenza, e per calcolare il numero appropriato di cicli di PCR necessari per amplificare sufficientemente ogni biblioteca senza causare sovracampionamento. In totale, quattro librerie sono state amplificate con indici unici di 6 bp, e vagliate per il contenuto endogeno sulla piattaforma Illumina MiSeq usando un sequenziamento SE di 250 bp. Le migliori librerie sono state poi ri-sequenziate sulla piattaforma Illumina HiSeq 2500 utilizzando la tecnologia 80 bp SE.
Analisi bioinformatica
Tutte le analisi di mappatura e danni al DNA sono state eseguite all’interno della pipeline Paleomix 1.2.1214. Le letture sono state tagliate con AdapterRemoval 2.2.015 utilizzando le impostazioni predefinite, tranne la lunghezza minima di lettura che è stata impostata a 25 bp. Le letture sono state ispezionate utilizzando FastQC e allineate con BWA16 applicando l’algoritmo Backtrack e disabilitando la lunghezza iniziale del seme. Se le letture mappate in posizioni multiple o con punteggi di qualità di mappatura (MAPQ score da BWA) inferiori a 30, sono state rimosse utilizzando SAMtools17. I duplicati di sequenza sono stati rimossi usando MarkDuplicates di Picard (disponibile da: http://broadinstitute.github.io/picard) e l’allineamento finale è stato riallineato intorno agli indel usando GATK18. La deaminazione della citosina a uracile nel campione MCE1356 è stata valutata utilizzando l’output di mapDamage v2.0.619. I risultati di MapDamage non hanno mostrato alcun chiaro segnale di deaminazione nelle sequenze (Fig. S3 supplementare).
Analisi mitocondriale
Per determinare la discendenza materna di MCE1356, le letture di sequenziamento del DNA sono state mappate ai genomi mitocondriali di riferimento beluga (accesso GenBank: KY444734) e narvalo (accesso GenBank: NC_005279) e la copertura media è stata confrontata. Le letture degli otto campioni di beluga e degli otto campioni di narvalo che comprendevano il pannello di riferimento sono state mappate ai rispettivi genomi mitocondriali di riferimento. Abbiamo costruito le sequenze di consenso mitocondriali di MCE1356 e dei 16 campioni del pannello di riferimento con regioni coperte da più di cinque letture utilizzando BEDtools20, SAMtools17 e GATK18. Abbiamo creato due allineamenti di sequenza utilizzando ClustalW21, applicando le impostazioni predefinite, che includevano i 16 campioni del pannello di riferimento e la versione di MCE1356 mappata al genoma mitocondriale di riferimento beluga o la versione di MCE1356 mappata al genoma mitocondriale narvalo di riferimento. Abbiamo usato i due allineamenti per costruire reti di aplotipi con giunzione mediana22 usando Popart 1.723 (disponibile da: http://popart.otago.ac.nz), escludendo tutti i siti con indel o dati mancanti. Successivamente, sia la copertura post-mapping che le due reti di aplotipi sono state utilizzate per determinare la specie del lignaggio materno di MCE1356.
Analisi del DNA nucleare
Abbiamo mappato le letture di sequenziamento del DNA di tutti i campioni al genoma di riferimento della balena assassina (Orcinus orca) (GCA_000331955.2). Un genoma di alta qualità di balena beluga è stato recentemente pubblicato24, ma la mappatura a una delle due specie potenziali genitori potrebbe creare distorsioni nelle analisi. Quindi abbiamo mappato le letture al genoma della balena assassina, in quanto è ben assemblato e le balene assassine sono ancora relativamente vicino a beluga e narvali (tempo di divergenza di 12 MYA)25, ma abbastanza lontano per ridurre al minimo il rischio di introgressione che avrebbe complicato le nostre analisi.
Per tutti i filtri ulteriori abbiamo usato ANGSD v0.92326, un pacchetto software che utilizza le probabilità di genotipo invece di genotipi chiamati, che è particolarmente utile quando si analizzano dati NGS a bassa copertura. Abbiamo usato il metodo SAMtools17 implementato in ANGSD per stimare le probabilità di genotipo, e inferito alleli maggiori e minori direttamente dalle probabilità di genotipo utilizzando un approccio di massima verosimiglianza come descritto in 27.
Per visualizzare i dati mappati, abbiamo tracciato la distribuzione di profondità di lettura da tutti gli individui, esclusi i siti con più di due alleli e siti con un punteggio Phred sotto 25. La visualizzazione dei dati di tutti gli individui combinati ci ha permesso di stimare la profondità di lettura media di 4,14x e di identificare i siti con una profondità di lettura elevata. Tali siti avevano maggiori probabilità di provenire da paraloghi e regioni ripetitive del genoma. Il dataset è stato ispezionato visivamente e ulteriormente filtrato, escludendo tutti i siti con profondità di lettura maggiore di nove (6,9% dei siti). In ANGSD, gli SNP sono stati chiamati in base alle loro frequenze alleliche. La frequenza allelica minore (MAF) è stata stimata dalle probabilità del genotipo e un test di likelihood ratio è stato utilizzato per determinare se la MAF era diversa da zero. Se il valore p del test del rapporto di verosimiglianza era <1e-4 il sito è stato considerato polimorfo e mantenuto nel dataset. Applicare questo filtro SNP significava che nessun sito con meno di quattro letture veniva mantenuto, poiché i siti coperti da meno letture non potevano ottenere valori di p così bassi. Questi filtri sono stati applicati a un set di dati comprendente tutti i 17 campioni e un set di dati senza MCE1356.
Per determinare se MCE1356 fosse un ibrido beluga/narwhal, abbiamo ulteriormente filtrato il set di dati di 17 individui, escludendo i siti senza letture in MCE1356. A questo punto la profondità media di lettura di MCE1356 era solo 1,15x, quindi al fine di garantire che non stavamo analizzando loci multicopia, abbiamo escluso i siti coperti da più di una lettura in MCE1356.
Il nostro obiettivo era quello di confrontare gli alleli trovati in MCE1356 agli alleli nel pannello di riferimento, quindi abbiamo stimato la probabilità di assegnare l’allele trovato in MCE1356 alla specie parentale sbagliata dato diverse profondità di lettura unica minima del pannello di riferimento e le frequenze alleliche. La profondità di lettura unica è stata definita come il numero di letture che coprono un sito specifico, dove tutte le letture provengono da un unico individuo. Questa probabilità P è stata calcolata come nell’equazione 1:
dove f(ps) è la frequenza allelica nella specie madre, e urd(psp) è la profondità di lettura unica specifica della specie nel pannello di riferimento.
La frequenza allelica della specie parentale che dà la massima probabilità f(ps-max) potrebbe essere descritta come nell’equazione 2:
Inserendo l’equazione 2 nell’equazione 1, la probabilità massima P(max) di assegnare l’allele trovato in MCE1356 alla specie parentale sbagliata è stata calcolata come nell’equazione 3:
I risultati hanno rivelato che con una profondità di lettura unica di due, tre e quattro in ogni specie parentale, la massima probabilità di assegnare l’allele trovato in MCE1356 alla specie parentale sbagliata era 0,148, 0,105 e 0,105.148, 0,105 e 0,082, rispettivamente. Questi valori massimi non dovrebbero essere interpretati come tassi di errore, in quanto si otterrebbero solo se tutti i siti fossero variabili all’interno delle specie parentali, e tutti i siti avessero un MAF di esattamente (1 – (urd/urd + 1)). Inoltre, assumendo che le distribuzioni MAF nei beluga e nei narvali fossero simili, l’assegnazione errata di alleli trovati in MCE1356 influenzerebbe entrambe le specie allo stesso modo, e quindi avrebbe un’influenza minima sulle inferenze di ibridazione. Un vantaggio dell’uso della profondità di lettura unica nelle equazioni 2 e 3 era che combinava il numero di individui e il numero di letture nella stima della probabilità di assegnare l’allele trovato in MCE1356 alla specie parentale sbagliata. Successivamente, la distribuzione della profondità di lettura, il MAF (con allele maggiore e minore fissi) e il numero di individui con dati in ogni sito sono stati calcolati separatamente per ogni specie parentale.
Sono stati costruiti tre set di dati, che oltre a MCE1356 includevano pannelli di specie parentali con profondità di lettura unica minima di due, tre e quattro. Il numero di siti conservati nei tre set di dati era 61.105, 11.764 e 360, rispettivamente. Come compromesso tra la massimizzazione del numero di siti e la minimizzazione dell’assegnazione errata degli alleli trovati in MCE1356, abbiamo scelto di utilizzare il set di dati con una lettura in MCE1356 e profondità di lettura unica minima di tre in ogni specie parentale. Quel set di dati comprendeva 11.764 siti, che sono stati utilizzati nelle analisi successive.
Sono state eseguite statistiche riassuntive sul set di dati con 17 individui, compreso il numero di siti che erano (i) fissati per alleli diversi nei pannelli di specie beluga e narvalo; (ii) polimorfo in beluga, ma non in narvalo; (iii) polimorfo in narvalo, ma non in beluga; (iv) polimorfo in entrambi beluga e narvalo. I siti che sono stimati per essere fissati tra i pannelli delle due specie parentali saranno arricchiti per i marcatori che sono altamente differenziati, cioè hanno grandi differenze di frequenza allelica, tra le due specie parentali. Questi marcatori, anche se non sono necessariamente differenze fisse tra le due specie parentali, sono ancora altamente informativi per l’ascendenza in MCE1356.
Abbiamo usato le probabilità del genotipo dal dataset senza MCE1356 per verificare che i beluga e i narvali nel pannello di riferimento non fossero essi stessi individui recentemente mescolati. Abbiamo stimato i loro coefficienti di commistione individuali specificando due popolazioni (K = 2) utilizzando NgsAdmix28. Sono state eseguite cento corse e la media e la deviazione standard dei coefficienti di commistione sono state utilizzate per la successiva interpretazione. Per confermare che i filtri non avevano rivelato profili genetici di commistione precedentemente nascosti nel pannello di riferimento, questa analisi è stata eseguita sia prima che dopo l’applicazione dei filtri di profondità di lettura unica.
Abbiamo analizzato le proporzioni di commistione di MCE1356 utilizzando fastNGSadmix29, applicando 1.000 bootstraps. fastNGSadmix utilizza le frequenze alleliche delle popolazioni/specie di riferimento e le probabilità di genotipo di un singolo individuo per stimare le proporzioni di commistione. Il software si è dimostrato robusto con i dati NGS con copertura a partire da 0,00015×29, ed era quindi ideale per il nostro studio.
Abbiamo ulteriormente stimato lo stato ibrido di MCE1356 indagando i siti fissati per diversi alleli nel pannello beluga e il pannello narvalo (9.178 siti), e confrontando la proporzione osservata di letture corrispondenti all’allele specifico beluga e l’allele specifico narvalo in MCE1356 alle proporzioni previste sotto diversi scenari di ibridazione. Per determinare quanto bene sette diversi scenari di ibridazione corrispondevano ai dati osservati, abbiamo calcolato un Chi-quadro di Pearson sulla bontà della statistica di adattamento. La statistica di test è calcolata come nell’equazione 4:
dove Oi ed Ei sono i conteggi osservati e previsti degli alleli derivati dalle specie parentali i, rispettivamente. Sotto l’ipotesi nulla in cui lo scenario scelto corrisponde bene ai dati osservati, la statistica di prova T segue una distribuzione centrale χ2. Quindi, lo scenario che corrisponde meglio ai dati osservati porterebbe alla statistica di test più bassa.
Per indagare ulteriormente i sette diversi scenari di ibridazione, abbiamo calcolato la probabilità degli alleli osservati nei siti che erano fissati per alleli diversi nelle popolazioni parentali. In particolare, abbiamo calcolato la probabilità degli alleli osservati nell’ibrido MCE1356, sotto un modello binomiale per l’eredità degli alleli dalle specie parentali, assumendo inoltre l’indipendenza tra i marcatori. Vorremmo notare, tuttavia, che la violazione dell’ipotesi di indipendenza porta ancora a stime imparziali. Supponendo che l’ibrido MCE1356 sia composto da una proporzione b di ascendenza beluga e (1-b) di ascendenza narvalo, possiamo scrivere la probabilità di b come prodotto della probabilità in k siti indipendenti, dato il numero di letture nib che corrispondono all’allele beluga nel sito i e nin che corrispondono all’ascendenza narvalo, come in Equazione 5.
Abbiamo calcolato la log-likelihood di b in tutto il suo intervallo da 0 a 1, e abbiamo confrontato le log-likelihoods attraverso i sette diversi scenari di ibridazione.
Determinazione del sesso
Per determinare il sesso di MCE1356 e degli individui nei pannelli di riferimento beluga e narvalo, abbiamo studiato i rapporti di copertura del cromosoma X rispetto agli autosomi. Questo è stato fatto confrontando la copertura tra gli scaffold originati putativamente dal cromosoma X e quella degli scaffold originati putativamente dagli autosomi. Abbiamo determinato quali scaffold provengono dal cromosoma X allineando il genoma della balena assassina al cromosoma X della mucca (Bos taurus) (CM008168.2) con SatsumaSynteny30. Inoltre, per rimuovere le distorsioni che possono verificarsi a causa di omologia tra alcune regioni dei cromosomi X e Y, abbiamo ulteriormente allineato le impalcature risultanti cromosoma X putativo al cromosoma Y umano (NC_000024.10) e rimosso queste impalcature da ulteriori analisi. Abbiamo poi calcolato la copertura attraverso le rimanenti impalcature X e le quattro più grandi impalcature non allineando al cromosoma X (cioè impalcature autosomiche), utilizzando la funzione di profondità in SAMtools17. Per compensare eventuali variazioni di copertura nel genoma, abbiamo selezionato a caso 10.000.000 di siti, calcolato la copertura media per questi siti e ripetuto questo passaggio per 100 volte. Per ogni individuo, sono stati calcolati gli intervalli di confidenza del 5% e del 95%, così come il primo e il terzo quartile, utilizzati per la successiva interpretazione.
Analisi degli isotopi stabili di carbonio e azoto
Campioni di ossa in polvere (~100 mg) da MCE1356, 18 crani di beluga e 18 crani di narvalo sono stati trattati con 10 ml 2:1 di cloroformio/metanolo (v/v) sotto sonicazione per 1 ora per rimuovere i lipidi. Dopo aver rimosso il solvente, i campioni sono stati essiccati a pressione atmosferica normale per 24 ore. I campioni sono stati poi demineralizzati in 10 ml di HCl 0,5 M per 4 ore mentre erano agitati da un agitatore orbitale. Dopo la demineralizzazione, i campioni sono stati risciacquati fino alla neutralità con acqua di tipo I, e poi riscaldati a 75 °C per 36 ore in 10-3 M HCl per solubilizzare il collagene. Il collagene solubile in acqua è stato poi liofilizzato. L’estratto di collagene è stato poi trattato con 10 ml 10:5:4 cloroformio/metanolo/acqua (v/v/v) sotto sonicazione per 1 h per rimuovere eventuali lipidi residui. Dopo la centrifugazione, lo strato di cloroformio/metanolo è stato rimosso e il metanolo è stato fatto evaporare dallo strato di acqua a 60 °C per 24 ore. I campioni sono stati nuovamente liofilizzati e pesati in capsule di stagno per le analisi elementari e isotopiche. Le composizioni isotopiche ed elementari stabili di carbonio e azoto sono state determinate utilizzando uno spettrometro di massa a flusso continuo Nu Horizon (Nu Instruments, Regno Unito) per il rapporto isotopico alla Trent University, Canada. Le composizioni isotopiche stabili di carbonio e azoto sono state calibrate rispetto alle scale VPDB e AIR utilizzando USGS40 e USGS41a. L’incertezza standard è stata determinata in ±0,17‰ per δ13C e ±0,22‰ per δ15N31; ulteriori dettagli analitici sono forniti nel documento supplementare S4. La significatività statistica delle differenze tra le composizioni isotopiche del beluga e del narvalo sono state valutate utilizzando test t non appaiati.