Amostragem
Analisamos (i) dados de ADN de todo o genoma extraído de dentes do crânio (MCE1356, Fig. 2b, Suplemento da Fig. S1) e de amostras de tecido de oito belugas e oito narvais, amostrados em Disko Bay, West Greenland, e (ii) composições isotópicas estáveis de carbono e nitrogênio do colágeno ósseo do MCE1356 e um painel de referência de 18 belugas e 18 narvais da West Greenland (Fig. 1c). As amostras de tecido (pele) foram armazenadas em etanol a 96%. As amostras foram coletadas por cientistas do Instituto de Recursos Naturais da Groenlândia sob a licença geral para amostragem biológica dos inuit do governo da Groenlândia e exportadas para a Dinamarca sob a licença CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 e 116.2006. As informações das amostras estão detalhadas na Tabela Complementar S2.
Exploração de ADN, preparação da biblioteca e sequenciamento
Amostras de tecidos
ADN de amostras de tecidos foi extraído usando o Qiagen Blood and Tissue Kit seguindo o protocolo do fabricante. O DNA foi fragmentado usando o Covaris M220 Focused-ultrasonicator para criar ~350-550 pares de base (bp) comprimentos de fragmentos. As bibliotecas foram construídas a partir dos extratos de DNA fragmentados usando o Illumina NeoPrep seguindo o NeoPrep Library Prep System Guide aplicando as configurações padrão. A amplificação, quantificação e normalização de PCR foram todas realizadas pelo Sistema de Preparação da Biblioteca NeoPrep. As bibliotecas foram rastreadas para distribuição de tamanho em um Bioanalisador Agilent 2100 e agrupadas em proporções equimolares antes do sequenciamento em um Illumina HiSeq 2500 com tecnologia 80 bp SE.
Amostras de dente/osso
Amostras de tecido, amostras antigas de osso e dentes têm concentrações relativamente baixas de DNA, o que requer diferentes protocolos de extração e construção de bibliotecas. Aproximadamente 0,5 g de pó ósseo de cinco dentes e um fragmento de osso da amostra MCE1356 foi perfurado usando um Dremel manual. O DNA foi extraído do pó de osso/dente em um laboratório dedicado à limpeza do DNA antigo no Museu de História Natural da Dinamarca, Universidade de Copenhagen, usando o tampão de extração descrito em11 com a adição de um estágio pré-digestino de 30 minutos12. Em vez de usar colunas Zymo-Spin V (Zymo Research), o tampão de extracção foi concentrado usando Amicon Ultra 30 kDa Centrifugal Filter Units e mais concentrado e limpo usando tubos de Qiagen Minelute. Os extratos purificados foram então incorporados nas bibliotecas do Illumina seguindo o protocolo descrito por13. Utilizámos o qPCR para verificar se a construção da biblioteca foi bem sucedida, para seleccionar quais as bibliotecas a sequenciar e para calcular o número apropriado de ciclos de PCR necessários para amplificar suficientemente cada biblioteca sem causar sobreamplificação. No total, quatro bibliotecas foram amplificadas com índices únicos de 6 bp, e rastreadas para conteúdo endógeno na plataforma Illumina MiSeq usando seqüenciamento de 250 bp SE. As melhores bibliotecas foram então re-sequenciadas na plataforma Illumina HiSeq 2500 usando tecnologia 80 bp SE.
Análise bioinformática
Todos os mapeamentos e análises de danos de DNA foram realizados dentro do pipeline Paleomix 1.2.1214. As leituras foram aparadas com o AdapterRemoval 2.2.015 usando as configurações padrão, exceto o comprimento mínimo de leitura, que foi definido para 25 bp. As leituras foram inspecionadas usando FastQC e alinhadas com o BWA16 aplicando o algoritmo Backtrack e desabilitando o comprimento da semente inicial. Se as leituras fossem mapeadas para múltiplas posições ou tivessem pontuação de qualidade de mapeamento (pontuação MAPQ do BWA) inferior a 30, elas eram removidas usando o SAMtools17. Duplicações sequenciais foram removidas usando MarkDuplicates do Picard (disponível em: http://broadinstitute.github.io/picard) e o alinhamento final foi realinhado em torno de indels usando GATK18. A desaminação de citosina para uracil na amostra MCE1356 foi avaliada usando a saída do mapDamage v2.0.619. Os resultados do MapDamage não mostraram nenhum sinal claro de desaminação nas sequências (Figura Suplementar S3).
Análise mitocondrial
Para determinar a linhagem materna do MCE1356, as leituras de sequenciamento de DNA foram mapeadas para os genomas de referência mitocondrial beluga (adesão ao GenBank: KY444734) e narval (adesão ao GenBank: NC_005279) e a cobertura média foi comparada. As leituras das oito amostras de beluga e oito narvais que compunham o painel de referência foram mapeadas para seus respectivos genomas de referência mitocondrial. Foram construídas seqüências de consenso mitocondrial de MCE1356 e as 16 amostras do painel de referência com regiões cobertas por mais de cinco leituras usando BEDtools20, SAMtools17 e GATK18. Criamos dois alinhamentos sequenciais usando ClustalW21, aplicando configurações padrão, que incluíram as 16 amostras de painel de referência e tanto a versão do MCE1356 mapeada para o genoma de referência mitocondrial beluga quanto a versão do MCE1356 mapeada para o genoma de referência mitocondrial narval. Utilizamos os dois alinhamentos para construir redes de haplótipos de união mediana22 usando o Popart 1.723 (disponível em: http://popart.otago.ac.nz), excluindo quaisquer sites com indels ou dados ausentes. Posteriormente, tanto a cobertura pós mapeamento como as duas redes de haplótipos foram usadas para determinar as espécies de MCE1356′s linhagem materna.
Análises de DNA nuclear
Mapeamos as leituras sequenciais de DNA de todas as amostras para o genoma de referência da baleia assassina (Orcinus orca) (GCA_000331955.2). Um genoma de baleia beluga de alta qualidade foi recentemente publicado24 , mas o mapeamento para uma das duas espécies parentais em potencial poderia criar enviesamentos nas análises. Por isso mapeamos as leituras para o genoma da baleia assassina, uma vez que esta está bem montada e as baleias assassinas ainda estão relativamente relacionadas com belugas e narvais (tempo de divergência de 12 MYA)25, mas distantes o suficiente para minimizar o risco de introgressão que complicaria as nossas análises.
Para toda a filtragem adicional usamos o ANGSD v0.92326, um pacote de software que usa probabilidades de genótipo em vez dos chamados genótipos, o que é particularmente útil na análise de dados NGS de baixa cobertura. Usamos o método SAMtools17 implementado em ANGSD para estimar as verosimilhanças genotípicas, e inferimos alelos maiores e menores diretamente das verosimilhanças genotípicas usando uma abordagem de máxima verosimilhança como descrita em27.
Para visualizar os dados mapeados, plotamos a distribuição de profundidade de leitura de todos os indivíduos, excluindo locais com mais de dois alelos e locais com pontuação de Phred abaixo de 25. A visualização dos dados de todos os indivíduos combinados permitiu-nos estimar a profundidade de leitura média de 4,14x e identificar locais com profundidade de leitura elevada. Tais locais eram mais prováveis de vir de parálogos e regiões repetitivas do genoma. O conjunto de dados foi inspecionado visualmente e filtrado, excluindo todos os locais com profundidade de leitura superior a nove (6,9% dos locais). Em ANGSD, os SNPs foram chamados com base em suas frequências de alelos. A frequência dos alelos menores (MAF) foi estimada a partir das probabilidades do genótipo e um teste de razão de probabilidade foi usado para determinar se o MAF era diferente de zero. Se o valor de p do teste do fator de correção foi <1e-4, o local foi considerado polimórfico e retido no conjunto de dados. A aplicação deste filtro SNP significou que nenhum local com menos de quatro leituras foi retido, pois locais cobertos por menos leituras não puderam obter valores de p tão baixos. Estes filtros foram aplicados a um conjunto de dados incluindo todas as 17 amostras, e a um conjunto de dados sem MCE1356.
Para determinar se MCE1356 era um híbrido beluga/narwhal, filtramos ainda o conjunto de dados de 17 indivíduos, excluindo locais sem leituras em MCE1356. Neste ponto a profundidade média de leitura de MCE1356 foi de apenas 1,15x, então para garantir que não estávamos analisando loci multicópia, excluímos locais cobertos por mais de uma leitura em MCE1356.
Nosso objetivo era comparar os alelos encontrados em MCE1356 com os alelos do painel de referência, então estimamos a probabilidade de atribuir o alelo encontrado em MCE1356 às espécies parentais erradas, dadas as diferentes profundidades mínimas de leitura e freqüências dos alelos do painel de referência. A profundidade de leitura única foi definida como o número de leituras cobrindo um local específico, onde todas as leituras vieram de um indivíduo único. Esta probabilidade P foi calculada como na Equação 1:
Onde f(ps) é a frequência alela na espécie parental, e urd(psp) é a profundidade de leitura única específica da espécie no painel de referência.
A frequência de alelos da espécie parental que dá a maior probabilidade f(ps-max) poderia ser descrita como na Equação 2:
Inserindo a Equação 2 na equação 1, a probabilidade máxima P(max) de atribuir o alelo encontrado em MCE1356 à espécie parental errada foi calculada como na Equação 3:
Resultados revelaram que com uma profundidade de leitura única de dois, três e quatro em cada espécie parental, a probabilidade máxima de atribuir o alelo encontrado no MCE1356 à espécie parental errada era de 0.148, 0,105 e 0,082, respectivamente. Esses valores máximos não devem ser interpretados como taxas de erro, pois eles só seriam obtidos se todos os locais fossem variáveis dentro da espécie parental e todos os locais tivessem um MAF de exatamente (1 – (urd/urd + 1)). Além disso, supondo que as distribuições de MAF nas belugas e narvais fossem semelhantes, a atribuição errada dos alelos encontrados no MCE1356 afectaria igualmente ambas as espécies e, portanto, teria uma influência mínima nas inferências de hibridação. Um benefício de usar a profundidade de leitura única nas equações 2 e 3 foi que combinou o número de indivíduos e o número de leituras na estimativa da probabilidade de atribuir o alelo encontrado no MCE1356 à espécie parental errada. Posteriormente, a distribuição da profundidade de leitura, MAF (com alelo maior e menor fixo), e número de indivíduos com dados em cada local foram calculados separadamente para cada espécie parental.
Foram construídos três conjuntos de dados, que além de MCE1356 incluíam painéis de espécies parentais com profundidade de leitura mínima única de dois, três e quatro. O número de locais retidos nos três conjuntos de dados foi de 61.105, 11.764 e 360, respectivamente. Como um compromisso entre maximizar o número de locais e minimizar a atribuição errada dos alelos encontrados em MCE1356, escolhemos usar o conjunto de dados com uma leitura em MCE1356 e profundidades mínimas de leitura únicas de três em cada espécie parental. Esse conjunto de dados incluiu 11.764 locais, que foram utilizados em análises subseqüentes.
As estatísticas sumárias foram realizadas no conjunto de dados com 17 indivíduos, incluindo o número de locais que foram (i) fixados para diferentes alelos nos painéis das espécies beluga e narval; (ii) polimórficos em belugas, mas não em narvais; (iii) polimórficos em narvais, mas não em belugas; (iv) polimórficos tanto em belugas como em narvais. Os locais que se estimam fixos entre os painéis das duas espécies parentais serão enriquecidos com marcadores altamente diferenciados, ou seja, com grandes diferenças de frequência alélica, entre as duas espécies parentais. Estes marcadores, embora não necessariamente diferenças fixas entre as duas espécies parentais, ainda são altamente informativos para a ancestralidade em MCE1356.
Usamos as probabilidades genotípicas do conjunto de dados sem MCE1356 para verificar que as belugas e os narvais no painel de referência não eram, eles próprios, indivíduos recentemente misturados. Estimamos seus coeficientes de mistura individuais especificando duas populações (K = 2) usando NgsAdmix28. Foram realizadas cem corridas e a média e o desvio padrão dos coeficientes de mistura foram utilizados para interpretação posterior. Para confirmar que os filtros não tinham revelado perfis genéticos de mistura previamente ocultos no painel de referência, esta análise foi realizada antes e depois da aplicação dos filtros de profundidade de leitura únicos.
Analisamos as proporções de mistura de MCE1356 usando fastNGSadmix29, aplicando 1.000 bootstraps. fastNGSadmix usa frequências de alelos de populações/espécies de referência e as probabilidades genotípicas de um único indivíduo para estimar suas proporções de mistura. O software provou ser robusto com dados NGS com cobertura tão baixa quanto 0,00015×29, e portanto foi ideal para nosso estudo.
Estimamos ainda o status híbrido do MCE1356 investigando locais fixados para diferentes alelos no painel beluga e no painel narval (9.178 locais), e comparando a proporção observada de leituras combinando o alelo específico beluga e o alelo específico narval em MCE1356 com as proporções esperadas sob diferentes cenários de hibridização. Para determinar quão bem sete diferentes cenários de hibridização corresponderam aos dados observados, calculamos uma estatística de adequação do Qui-quadrado de Pearson. A estatística do teste é calculada como na Equação 4:
onde Oi e Ei são as contagens observadas e esperadas de alelos derivados da espécie parental i, respectivamente. Sob a hipótese nula onde o cenário escolhido corresponde bem aos dados observados, a estatística de teste T segue uma distribuição central χ2. Assim, o cenário que melhor corresponde aos dados observados levaria à estatística de teste mais baixa.
Para investigar melhor os sete diferentes cenários de hibridação, calculamos a probabilidade dos alelos observados em locais que foram fixados para diferentes alelos nas populações parentais. Especificamente, calculamos a probabilidade dos alelos observados no híbrido MCE1356, sob um modelo binomial para a herança dos alelos da espécie parental, assumindo ainda mais a independência entre os marcadores. Gostaríamos de notar, no entanto, que a violação da suposição da independência ainda leva a estimativas imparciais. Assumindo que o híbrido MCE1356 é composto por uma proporção b de ancestralidade beluga e (1-b) de ancestralidade narval, podemos escrever a probabilidade de b como produto da probabilidade em k locais independentes, dado o número de bicos de leitura que correspondem ao alelo beluga no local i e nin que correspondem à ancestralidade narval, como na Equação 5.
Calculamos a probabilidade logarítmica de b em toda a sua amplitude, de 0 a 1, e comparou os log-likelihoods através dos sete diferentes cenários de hibridação.
Determinação do sexo
Para determinar o sexo do MCE1356 e os indivíduos nos painéis de referência beluga e narval, investigamos as relações entre o cromossoma X e a cobertura autossômica. Isto foi feito comparando a cobertura entre os andaimes putativamente originários do cromossomo X, com a dos andaimes putativamente originários dos autossomos. Nós determinamos quais andaimes se originam do cromossomo X alinhando o genoma da baleia assassina com o cromossomo X da Vaca (Bos taurus) (CM008168.2) com o SatsumaSynteny30. Além disso, para remover os vieses que podem ocorrer devido à homologia entre algumas regiões dos cromossomos X e Y, nós alinhamos ainda mais os supostos andaimes X resultantes ao cromossomo Y humano (NC_000024.10) e removemos esses andaimes de análises posteriores. Em seguida calculamos a cobertura através dos andaimes X restantes e dos quatro maiores andaimes não alinhados ao cromossoma X (ou seja, andaimes autossômicos), usando a função de profundidade em SAMtools17. Para compensar qualquer variação na cobertura ao longo do genoma, selecionamos aleatoriamente 10.000.000 locais, calculamos a cobertura média para estes locais, e repetimos este passo 100 vezes. Para cada indivíduo, foram calculados intervalos de confiança de 5% e 95%, assim como o primeiro e terceiro quartis, e utilizados para interpretação posterior.
Análise de isótopos estáveis de carbono e nitrogênio
Amostras de osso em pó (~100 mg) de MCE1356, 18 crânios de beluga e 18 crânios de narval foram tratados com 10 ml de clorofórmio 2:1/metanol (v/v) sob sonicação durante 1 h para remoção de lipídios. Após a remoção do solvente, as amostras foram secas sob pressão atmosférica normal durante 24 h. As amostras foram então desmineralizadas em 10 ml de HCl 0,5 M durante 4 h enquanto eram agitadas por um agitador orbital. Após a desmineralização, as amostras foram enxaguadas para neutralidade com água tipo I, e então aquecidas a 75 °C por 36 h em HCl 10-3 M para solubilizar o colágeno. O colágeno solúvel em água foi então liofilizado. O extrato de colágeno foi então tratado com 10 ml 10:5:4 clorofórmio/metanol/água (v/v/v) sob sonicação por 1 h para remover quaisquer lipídios residuais. Após a centrifugação, a camada de clorofórmio/metanol foi removida e o metanol foi evaporado da camada de água a 60 °C durante 24 h. As amostras foram novamente liofilizadas e pesadas em cápsulas de estanho para análise elementar e isotópica. As composições isotópicas e elementares estáveis em carbono e nitrogênio foram determinadas usando um espectrômetro de massa de razão isotópica de fluxo contínuo Nu Horizon (Nu Instruments, Reino Unido) na Universidade de Trent, Canadá. As composições isotópicas estáveis de carbono e nitrogênio foram calibradas em relação às escalas VPDB e AIR usando USGS40 e USGS41a. A incerteza padrão foi determinada em ±0.17‰ para δ13C e ±0.22‰ para δ15N31; detalhes analíticos adicionais são fornecidos no Documento Suplementar S4. A significância estatística das diferenças entre as composições isotópicas beluga e narval foi avaliada utilizando testes t não pareados.