Échantillonnage
Nous avons analysé (i) les données d’ADN à l’échelle du génome extraites des dents du crâne (MCE1356, Fig. 2b, Fig. supplémentaire S1) et d’échantillons de tissus de huit bélugas et huit narvals échantillonnés à Disko Bay, Groenland occidental, et (ii) les compositions isotopiques stables du carbone et de l’azote du collagène osseux de MCE1356 et d’un panel de référence de 18 crânes de bélugas et 18 crânes de narvals du Groenland occidental (Fig. 1c). Les échantillons de tissus (peau) ont été conservés dans de l’éthanol à 96%. Les échantillons ont été collectés par des scientifiques de l’Institut des ressources naturelles du Groenland en vertu du permis général d’échantillonnage biologique des Inuits du gouvernement du Groenland et exportés au Danemark en vertu des permis CITES IM 0401-897/04, IM 0721-199/08, IM 0330-819/09 et 116.2006. Les informations sur les échantillons sont détaillées dans le tableau supplémentaire S2.
Extraction d’ADN, préparation de la bibliothèque et séquençage
Échantillons de tissus
L’ADN des échantillons de tissus a été extrait à l’aide du kit de sang et de tissus Qiagen en suivant le protocole du fabricant. L’ADN a été fragmenté à l’aide de l’appareil à ultrasons focalisés Covaris M220 pour créer des fragments d’une longueur de ~350-550 paires de bases (pb). Les librairies ont été construites à partir des extraits d’ADN fragmentés à l’aide du système NeoPrep d’Illumina en suivant le guide du système NeoPrep Library Prep et en appliquant les paramètres par défaut. L’amplification par PCR, la quantification et la normalisation ont toutes été effectuées par le système NeoPrep Library Prep. Les bibliothèques ont été criblées pour la distribution de taille sur un bioanalyseur Agilent 2100 et regroupées dans des rapports équimolaires avant le séquençage sur un Illumina HiSeq 2500 avec une technologie SE de 80 pb.
Échantillons de dents/os
Contrairement aux échantillons de tissus, les échantillons de vieux os et de dents ont des concentrations d’ADN relativement faibles, ce qui nécessite des protocoles d’extraction et de construction de bibliothèque différents. Environ 0,5 g de poudre d’os provenant de cinq dents et un tesson d’os du spécimen MCE1356 ont été percés à l’aide d’un Dremel manuel. L’ADN a été extrait de la poudre de dent/d’os dans un laboratoire dédié à l’ADN ancien au Musée d’histoire naturelle du Danemark, Université de Copenhague, en utilisant le tampon d’extraction décrit dans11 avec l’ajout d’une étape de prédigestion de 30 minutes12. Au lieu d’utiliser des colonnes Zymo-Spin V (Zymo Research), le tampon d’extraction a été concentré à l’aide d’unités de filtration centrifuge Amicon Ultra 30 kDa, puis concentré et nettoyé à l’aide de tubes Minelute de Qiagen. Les extraits purifiés ont ensuite été construits en bibliothèques Illumina en suivant le protocole décrit par13. Nous avons utilisé la qPCR pour vérifier que la construction de la bibliothèque était réussie, pour sélectionner les bibliothèques à séquencer et pour calculer le nombre approprié de cycles PCR nécessaires pour amplifier suffisamment chaque bibliothèque sans provoquer de suramplification. Au total, quatre librairies ont été amplifiées avec des index uniques de 6 pb, et criblées pour le contenu endogène sur la plateforme Illumina MiSeq en utilisant un séquençage SE de 250 pb. Les meilleures librairies ont ensuite été reséquencées sur la plateforme Illumina HiSeq 2500 à l’aide de la technologie SE 80 pb.
Analyse bioinformatique
Toutes les analyses de cartographie et de dommages à l’ADN ont été effectuées dans le pipeline Paleomix 1.2.1214. Les lectures ont été rognées avec AdapterRemoval 2.2.015 en utilisant les paramètres par défaut, sauf la longueur minimale de lecture qui a été fixée à 25 pb. Les lectures ont été inspectées avec FastQC et alignées avec BWA16 en appliquant l’algorithme Backtrack et en désactivant la longueur de la graine de départ. Si les lectures étaient mappées à des positions multiples ou avaient des scores de qualité de mappage (score MAPQ de BWA) inférieurs à 30, elles étaient supprimées à l’aide de SAMtools17. Les doublons de séquence ont été supprimés à l’aide de MarkDuplicates de Picard (disponible à partir de : http://broadinstitute.github.io/picard) et l’alignement final a été réaligné autour des indels à l’aide de GATK18. La désamination de la cytosine en uracile dans le spécimen MCE1356 a été évaluée en utilisant la sortie de mapDamage v2.0.619. Les résultats de mapDamage n’ont pas montré de signal clair de désamination dans les séquences (figure supplémentaire S3).
Analyse mitochondriale
Pour déterminer la lignée maternelle de MCE1356, les lectures de séquençage de l’ADN ont été mappées aux génomes de référence mitochondriaux du béluga (accession GenBank : KY444734) et du narval (accession GenBank : NC_005279) et la couverture moyenne a été comparée. Les lectures provenant des huit échantillons de bélugas et des huit échantillons de narvals qui constituaient le panel de référence ont été mises en correspondance avec leurs génomes mitochondriaux de référence respectifs. Nous avons construit des séquences consensus mitochondriales de MCE1356 et des 16 échantillons du panel de référence avec des régions couvertes par plus de cinq lectures en utilisant BEDtools20, SAMtools17 et GATK18. Nous avons créé deux alignements de séquences en utilisant ClustalW21, en appliquant les paramètres par défaut, qui incluaient les 16 échantillons du panel de référence et soit la version du MCE1356 cartographiée au génome mitochondrial de référence du béluga, soit la version du MCE1356 cartographiée au génome mitochondrial de référence du narval. Nous avons utilisé les deux alignements pour construire des réseaux d’haplotypes à jonction médiane22 à l’aide de Popart 1.723 (disponible auprès de : http://popart.otago.ac.nz), en excluant tous les sites comportant des indels ou des données manquantes. Par la suite, la couverture post-mapping et les deux réseaux d’haplotypes ont été utilisés pour déterminer l’espèce de la lignée maternelle de MCE1356′.
Analyses de l’ADN nucléaire
Nous avons mappé les lectures de séquençage de l’ADN de tous les échantillons au génome de référence de l’orque (Orcinus orca) (GCA_000331955.2). Un génome de béluga de haute qualité a été récemment publié24, mais la cartographie sur l’une des deux espèces parentales potentielles pourrait créer des biais dans les analyses. Nous avons donc mappé les lectures sur le génome de l’épaulard, car il est bien assemblé et les épaulards sont encore relativement proches des bélugas et des narvals (temps de divergence de 12 MYA)25, tout en étant suffisamment éloignés pour minimiser le risque d’introgression qui compliquerait nos analyses.
Pour tous les autres filtrages, nous avons utilisé ANGSD v0.92326, un logiciel qui utilise les probabilités de génotype au lieu des génotypes appelés, ce qui est particulièrement utile lors de l’analyse des données NGS à faible couverture. Nous avons utilisé la méthode SAMtools17 mise en œuvre dans ANGSD pour estimer les vraisemblances de génotype, et inféré les allèles majeurs et mineurs directement à partir des vraisemblances de génotype en utilisant une approche de vraisemblance maximale comme décrit dans27.
Pour visualiser les données cartographiées, nous avons tracé la distribution de la profondeur de lecture de tous les individus, en excluant les sites avec plus de deux allèles et les sites avec un score Phred inférieur à 25. La visualisation des données de tous les individus combinés nous a permis d’estimer la profondeur de lecture moyenne de 4,14x et d’identifier les sites avec une profondeur de lecture élevée. Ces sites étaient plus susceptibles de provenir de paralogues et de régions répétitives du génome. L’ensemble de données a été inspecté visuellement et filtré davantage, excluant tous les sites avec une profondeur de lecture supérieure à neuf (6,9% des sites). Dans ANGSD, les SNP ont été appelés sur la base de leurs fréquences alléliques. La fréquence des allèles mineurs (MAF) a été estimée à partir des vraisemblances des génotypes et un test de rapport de vraisemblance a été utilisé pour déterminer si la MAF était différente de zéro. Si la valeur p du test du rapport de vraisemblance était <1e-4, le site était considéré comme polymorphe et retenu dans l’ensemble de données. En appliquant ce filtre SNP, aucun site avec moins de quatre lectures n’a été retenu, car les sites couverts par moins de lectures ne pouvaient pas obtenir des valeurs p aussi basses. Ces filtres ont été appliqués à un ensemble de données comprenant les 17 échantillons, et un ensemble de données sans MCE1356.
Pour déterminer si MCE1356 était un hybride béluga/narval, nous avons encore filtré l’ensemble de données de 17 individus, en excluant les sites sans lectures dans MCE1356. A ce stade, la profondeur de lecture moyenne de MCE1356 n’était que de 1,15x, donc pour s’assurer que nous n’analysions pas des loci multicopies, nous avons exclu les sites couverts par plus d’une lecture dans MCE1356.
Notre objectif était de comparer les allèles trouvés dans MCE1356 aux allèles du panel de référence, donc nous avons estimé la probabilité d’attribuer l’allèle trouvé dans MCE1356 à la mauvaise espèce parentale compte tenu des profondeurs de lecture uniques minimales et des fréquences alléliques différentes du panel de référence. La profondeur de lecture unique a été définie comme le nombre de lectures couvrant un site spécifique, où toutes les lectures proviennent d’un individu unique. Cette probabilité P a été calculée comme dans l’équation 1 :
Où f(ps) est la fréquence de l’allèle dans l’espèce parentale, et urd(psp) est la profondeur de lecture unique spécifique à l’espèce dans le panel de référence.
La fréquence de l’allèle de l’espèce parentale donnant la probabilité la plus élevée f(ps-max) pourrait être décrite comme dans l’équation 2 :
En insérant l’équation 2 dans l’équation 1, la probabilité maximale P(max) d’attribuer l’allèle trouvé dans MCE1356 à la mauvaise espèce parentale a été calculée comme dans l’équation 3 :
Les résultats ont révélé qu’avec une profondeur de lecture unique de deux, trois et quatre dans chaque espèce parentale, la probabilité maximale d’attribuer l’allèle trouvé dans MCE1356 à la mauvaise espèce parentale était de 0.148, 0,105 et 0,082, respectivement. Ces valeurs maximales ne doivent pas être interprétées comme des taux d’erreur, car elles ne seraient obtenues que si tous les sites étaient variables au sein de l’espèce parentale, et que tous les sites avaient un MAF d’exactement (1 – (urd/urd + 1)). De plus, en supposant que les distributions de MAF chez les bélugas et les narvals étaient similaires, l’assignation erronée des allèles trouvés dans MCE1356 affecterait les deux espèces de manière égale, et aurait donc une influence minimale sur les inférences d’hybridation. Un avantage de l’utilisation de la profondeur de lecture unique dans les équations 2 et 3 est qu’elle combine le nombre d’individus et le nombre de lectures dans l’estimation de la probabilité d’attribuer l’allèle trouvé dans MCE1356 à la mauvaise espèce parentale. Par la suite, la distribution de la profondeur de lecture, le MAF (avec un allèle majeur et mineur fixe) et le nombre d’individus avec des données dans chaque site ont été calculés séparément pour chaque espèce parentale.
Trois ensembles de données ont été construits, qui en plus de MCE1356 comprenaient des panels d’espèces parentales avec des profondeurs de lecture uniques minimales de deux, trois et quatre. Le nombre de sites retenus dans les trois ensembles de données était de 61 105, 11 764 et 360, respectivement. Comme compromis entre la maximisation du nombre de sites et la minimisation de l’attribution erronée des allèles trouvés dans MCE1356, nous avons choisi d’utiliser le jeu de données avec une lecture dans MCE1356 et des profondeurs de lecture uniques minimales de trois dans chaque espèce parentale. Cet ensemble de données comprenait 11 764 sites, qui ont été utilisés dans les analyses ultérieures.
Des statistiques sommaires ont été effectuées sur l’ensemble de données avec 17 individus, y compris le nombre de sites qui étaient (i) fixés pour différents allèles dans les panels d’espèces de bélugas et de narvals ; (ii) polymorphes chez les bélugas, mais pas chez les narvals ; (iii) polymorphes chez les narvals, mais pas chez les bélugas ; (iv) polymorphes à la fois chez les bélugas et les narvals. Les sites qui sont estimés être fixes entre les deux panels d’espèces parentales seront enrichis de marqueurs qui sont hautement différenciés, c’est-à-dire qui présentent de grandes différences de fréquences alléliques, entre les deux espèces parentales. Ces marqueurs, même s’ils ne présentent pas nécessairement des différences fixes entre les deux espèces parentales, sont tout de même très informatifs pour l’ascendance dans MCE1356.
Nous avons utilisé les vraisemblances de génotype de l’ensemble de données sans MCE1356 pour vérifier que les bélugas et les narvals du panel de référence n’étaient pas eux-mêmes des individus récemment mélangés. Nous avons estimé leurs coefficients d’admixtion individuels en spécifiant deux populations (K = 2) à l’aide de NgsAdmix28. Cent passages ont été effectués et la moyenne et l’écart type des coefficients d’admixtion ont été utilisés pour une interprétation ultérieure. Pour confirmer que les filtres n’avaient pas révélé des profils génétiques mixtes cachés dans le panel de référence, cette analyse a été effectuée avant et après l’application des filtres de profondeur de lecture unique.
Nous avons analysé les proportions de mélange de MCE1356 en utilisant fastNGSadmix29, en appliquant 1 000 bootstraps. fastNGSadmix utilise les fréquences alléliques des populations/espèces de référence et les probabilités de génotype d’un individu unique pour estimer ses proportions de mélange. Le logiciel s’est avéré robuste avec les données NGS avec une couverture aussi faible que 0,00015×29, et était donc idéal pour notre étude.
Nous avons en outre estimé le statut hybride de MCE1356 en examinant les sites fixés pour différents allèles dans le panel de bélugas et le panel de narvals (9 178 sites), et en comparant la proportion observée de lectures correspondant à l’allèle spécifique au béluga et à l’allèle spécifique au narval dans MCE1356 aux proportions attendues dans différents scénarios d’hybridation. Pour déterminer dans quelle mesure sept scénarios d’hybridation différents correspondent aux données observées, nous avons calculé une statistique de qualité d’ajustement du chi carré de Pearson. La statistique de test est calculée comme dans l’équation 4 :
où Oi et Ei sont les comptes observés et attendus des allèles dérivés de l’espèce parentale i, respectivement. Sous l’hypothèse nulle où le scénario choisi correspond bien aux données observées, la statistique de test T suit une distribution centrale χ2. Ainsi, le scénario qui correspond le mieux aux données observées conduirait à la statistique de test la plus faible.
Pour approfondir les sept différents scénarios d’hybridation, nous avons calculé la vraisemblance des allèles observés sur des sites qui étaient fixés pour différents allèles dans les populations parentales. Plus précisément, nous avons calculé la probabilité des allèles observés dans l’hybride MCE1356, selon un modèle binomial pour l’héritage des allèles des espèces parentales, en supposant en outre l’indépendance entre les marqueurs. Nous tenons cependant à noter que la violation de l’hypothèse d’indépendance conduit toujours à des estimations non biaisées. En supposant que l’hybride MCE1356 est composé d’une proportion b d’ascendance béluga et (1-b) d’ascendance narval, nous pouvons écrire la vraisemblance de b comme un produit de la vraisemblance à k sites indépendants, étant donné le nombre de lectures nib qui correspondent à l’allèle béluga au site i et nin qui correspondent à l’ascendance narval, comme dans l’équation 5.
Nous avons calculé la log-vraisemblance de b sur toute sa plage, de 0 à 1, et avons comparé les log-vraisemblances à travers les sept différents scénarios d’hybridation.
Détermination du sexe
Pour déterminer le sexe de MCE1356 et des individus des panels de référence de bélugas et de narvals, nous avons étudié les rapports de couverture du chromosome X par rapport aux autosomes. Ceci a été fait en comparant la couverture à travers les échafaudages provenant putativement du chromosome X, à celle des échafaudages provenant putativement des autosomes. Nous avons déterminé quels échafaudages provenaient du chromosome X en alignant le génome de l’épaulard sur le chromosome X de la vache (Bos taurus) (CM008168.2) avec SatsumaSynteny30. De plus, afin d’éliminer les biais qui peuvent se produire en raison de l’homologie entre certaines régions des chromosomes X et Y, nous avons aligné les échafaudages putatifs du chromosome X sur le chromosome Y humain (NC_000024.10) et avons retiré ces échafaudages de l’analyse. Nous avons ensuite calculé la couverture des échafaudages X restants et des quatre plus grands échafaudages non alignés sur le chromosome X (c’est-à-dire les échafaudages autosomiques), en utilisant la fonction de profondeur de SAMtools17. Pour compenser toute variation de la couverture dans le génome, nous avons sélectionné au hasard 10 000 000 de sites, calculé la couverture moyenne de ces sites et répété cette étape 100 fois. Pour chaque individu, des intervalles de confiance à 5% et 95%, ainsi que les premiers et troisièmes quartiles, ont été calculés et utilisés pour l’interprétation ultérieure.
Analyse des isotopes stables du carbone et de l’azote
Des échantillons d’os en poudre (~100 mg) provenant de MCE1356, 18 crânes de béluga et 18 crânes de narval ont été traités avec 10 ml de chloroforme/méthanol 2:1 (v/v) sous sonication pendant 1 h pour éliminer les lipides. Après avoir éliminé le solvant, les échantillons ont été séchés sous pression atmosphérique normale pendant 24 h. Les échantillons ont ensuite été déminéralisés dans 10 ml de HCl 0,5 M pendant 4 h tout en étant agités par un agitateur orbital. Après déminéralisation, les échantillons ont été rincés jusqu’à neutralité avec de l’eau de type I, puis chauffés à 75 °C pendant 36 h dans du HCl 10-3 M pour solubiliser le collagène. Le collagène soluble dans l’eau a ensuite été lyophilisé. L’extrait de collagène a ensuite été traité avec 10 ml de 10:5:4 chloroforme/méthanol/eau (v/v/v) sous sonication pendant 1 h pour éliminer tout lipide résiduel. Après centrifugation, la couche de chloroforme/méthanol a été retirée et le méthanol a été évaporé de la couche d’eau à 60 °C pendant 24 h. Les échantillons ont été à nouveau lyophilisés et pesés dans des capsules d’étain pour l’analyse élémentaire et isotopique. Les compositions isotopiques et élémentaires stables du carbone et de l’azote ont été déterminées à l’aide d’un spectromètre de masse de rapport isotopique à flux continu Nu Horizon (Nu Instruments, UK) à l’Université de Trent, Canada. Les compositions isotopiques stables du carbone et de l’azote ont été calibrées par rapport aux échelles VPDB et AIR en utilisant USGS40 et USGS41a. L’incertitude standard a été déterminée comme étant de ±0,17‰ pour δ13C et de ±0,22‰ pour δ15N31 ; des détails analytiques supplémentaires sont fournis dans le document supplémentaire S4. La signification statistique des différences entre les compositions isotopiques des bélugas et des narvals a été évaluée à l’aide de tests t non appariés.