Échantillonnage
Le consentement éclairé écrit a été obtenu de l’individu sâme dont le génome a été analysé dans cette étude, Celle-ci a été approuvée par le comité d’éthique du district hospitalier d’Helsinki et d’Uusimaa (décision 329/13/03/00/2013) et par le comité d’éthique de l’université de Leipzig (numéro de référence d’approbation 398-13-16122013).
L’échantillonnage et l’extraction d’ADN ancien nécessitent une procédure stricte afin d’éviter toute contamination introduite par du matériel génétique contemporain. Pour les 13 individus de l’âge du fer originaires de Finlande dont nous disposions, le prélèvement a eu lieu dans une installation en salle blanche dédiée aux travaux sur l’ADN ancien, à l’Institut des sciences archéologiques de Tübingen. Le flux de travail préliminaire comprenait la documentation, la photographie et le stockage des échantillons dans des tubes en plastique et des sacs en plastique individuels portant un code d’identification. À la suite d’une première étude pilote, les échantillons de dents que nous avons utilisés étaient fragmentés et une partie de la dentine a été retirée. La dentine restante a été recueillie en la séparant soigneusement de l’émail à l’aide d’une perceuse de dentiste et de têtes de forage en diamant refroidies, tournées à une vitesse inférieure à 15 tr/min, afin d’éviter d’éventuels dommages causés par la chaleur à l’ADN ancien.
Pour les échantillons provenant des sites de Bolshoy et de Chalmny Varre, nous avons utilisé des restes de poudre de dent qui avaient été traités à l’origine à l’Institut d’anthropologie de l’Université de Mayence à des fins de réplication, comme décrit dans la réf. 24. En bref, les étapes de préparation des échantillons comprenaient une irradiation aux UV pendant 30 à 45 minutes, suivie d’un essuyage doux de la surface avec des agents de blanchiment commerciaux dilués. Les dents ont ensuite été sablées à l’aide d’un abrasif à l’oxyde d’aluminium (Harnisch & Rieth) et broyées en poudre fine à l’aide d’un mélangeur-mélangeur (Retsch).
Calibrage de la date radiocarbone
Nous avons calibré la date radiocarbone de Bolshoy, rapportée dans les réfs 24,55 comme 3473 ± 42 ans BP, en utilisant Intcal1356 comme courbe de calibrage, en utilisant OxCal 4.357.
Extraction d’ADN et préparation de la bibliothèque
L’ADN des six échantillons de Bolshoy et des deux échantillons de Chalmny Varre a été extrait dans les installations d’ADN ancien de l’Institut Max Planck pour la science de l’histoire humaine (MPI-SHH) à Iéna, en Allemagne. L’extraction des échantillons de Levänluhta a été réalisée de la même manière dans les salles blanches de l’Institut des sciences archéologiques de Tübingen. Pour chaque spécimen, ~ 50 mg de poudre de dentine ont été utilisés pour une procédure d’extraction spécifiquement conçue pour la récupération d’ADN ancien58. Un tampon d’extraction contenant 0,45 M EDTA, pH 8,0 (Life Technologies) et 0,25 mg/ml de protéinase K (Sigma-Aldrich) a été ajouté à la poudre d’os et incubé à 37 °C avec rotation pendant la nuit. Le surnageant a été séparé du culot de la poudre d’os par centrifugation (14 000 rpm). Un tampon de liaison composé de GuHCL 5 M (Sigma Aldrich) et d’Isopropanol 40 % (Merck), ainsi que 400 μl d’acétate de sodium 1 M (pH 5,5) ont été ajoutés au surnageant, et la solution a été purifiée en la faisant tourner à travers une colonne de purification fixée à un entonnoir High Pure Extender Assembly (8 min à 1500 rpm, avec une accélération lente). La colonne a ensuite été centrifugée dans un tube collecteur (1 min 14 000 rpm) 1 à 2 fois pour maximiser le rendement. Cette étape a été suivie de deux étapes de lavage ultérieures de 450 μl de tampon de lavage (High Pure Viral Nucleic Acid Large Volume Kit) et de deux étapes de spin sec de 1 min de centrifugation à 14 000 rpm. Le volume total final de 100 μl d’éluat a été atteint par deux tours d’élution séparés de 50 μl de TET (10 mM Tris-HCL, 1 mM EDTA pH 8,0, 0,1% Tween20), chacun étant essoré pendant 1 min à 14 000 rpm dans un tube Eppendorf de 1,5 ml frais. Les contrôles négatifs (tampon au lieu de l’échantillon) ont été traités en parallèle à raison d’un contrôle pour 7 échantillons.
Sur les 100 μl d’extrait, 20 μl ont été utilisés pour immortaliser l’ADN de l’échantillon sous forme de bibliothèque double brin. La procédure comprenait une étape de réparation de l’extrémité émoussée, de ligature de l’adaptateur et de remplissage de l’adaptateur, comme décrit par Meyer et Kircher59. Au cours de l’étape de réparation de l’extrémité émoussée, un mélange de 0,4 U/μl de PNK (polynucléotide kinase) T4 et de 0,024 U/μl d’ADN polymérase T4, de tampon 2 1× NEB (NEB), de mélange de dNTP 100 μM (Thermo Scientific), d’ATP 1 mM (NEB) et de 0.8 mg/ml de BSA (NEB) a été ajouté à l’ADN matrice, suivi d’une incubation dans un thermocycleur (15 min 15 °C, 15 min 25 °C) et d’une purification avec un kit MinElute (QIAGEN). Le produit a été élué dans 18 μl de tampon TET. L’étape de ligature des adaptateurs comprenait un mélange de 1× tampon Quick Ligase (NEB), 250 nM d’adaptateurs Illumina (Sigma-Aldrich) et 0,125 U/μl de Quick Ligase (NEB), ajouté aux 18 μl d’éluat, suivi d’une incubation de 20 min, et d’une deuxième étape de purification avec des colonnes MinElute, cette fois dans 20 μl d’éluat. Pour l’étape de remplissage, un mélange de 0,4 U/μl de Bst-polymérase et de 125 μM de mélange dNTP a été ajouté, puis le mélange a été incubé dans un thermocycleur (30 min 37 °C, 10 min 80 °C). Des bibliothèques sans traitement à l’Uracil-ADN-glycosylase (UDG) ont été produites pour les 13 extraits de Levänluhta. En outre, des librairies semi-traitées à l’UDG ont été produites pour sept des 13 extraits originaux de Levänluhta, et pour tous les extraits de Bolshoy et de Chalmny Varre. Pour introduire le traitement UDG-half, une étape initiale a été incluse dans la préparation de la bibliothèque, dans laquelle 250 U d’enzyme USER (NEB) ont été ajoutés dans les 20 μl d’extrait, suivis d’une incubation à 37 °C pendant 30 min, puis à 12 °C pendant 1 min. Cette opération a de nouveau été suivie de l’ajout de 200 UGI (inhibiteur de l’Uracil Glycosylase, par NEB) et d’une autre incubation identique pour arrêter l’excision enzymatique des sites désaminés, comme décrit dans60. Pour chaque bibliothèque, une paire unique d’index de huit pb de long a été incorporée à l’aide d’une ADN polymérase Pfu Turbo Cx Hotstart et d’un programme de thermocyclage dont le profil de température est le suivant : dénaturation initiale (98 °C pendant 30 sec), cycle de dénaturation/recuit/allongement (98 °C pendant 10 sec/ 60 °C pendant 20 sec/ 72 °C pendant 20 sec) et extension finale à 72 °C pendant 10 min61. De la poudre d’os provenant d’un ours des cavernes a été traitée en parallèle, servant de contrôle positif. Les contrôles négatifs pour les deux étapes d’extraction et de préparation des bibliothèques ont été conservés à côté des échantillons tout au long du flux de travail.
L’efficacité de l’expérience a été assurée en quantifiant la concentration des bibliothèques sur qPCR (Roche) en utilisant des aliquotes de bibliothèques avant et après l’indexation. Le nombre de copies moléculaires dans les bibliothèques pré-indexées allait de ~ 10E8 à ~ 10E9 copies/μl, indiquant une préparation réussie de la bibliothèque, tandis que les bibliothèques indexées allaient de ~ 10E10 à ~ 10E12 copies/μl, indiquant une efficacité d’indexation admissible. Les contrôles négatifs présentaient une concentration inférieure de 4 à 5 ordres de grandeur à celle des échantillons, indiquant de faibles niveaux de contamination provenant des étapes de traitement en laboratoire.
Les bibliothèques ont été amplifiées par PCR, pour la quantité de cycles correspondant aux concentrations des bibliothèques indexées, à l’aide de la polymérase AccuPrime Pfx (5 μl de matrice de bibliothèque, 2 U d’ADN polymérase AccuPrime Pfx d’Invitrogen, 1 U de mélange maître PCR 10× prêt à l’emploi et 0.3 μM d’amorces IS5 et IS6, pour chaque réaction de 100 μl) avec un profil thermique de 2 min de dénaturation à 95 °C, 3 à 9 cycles consistant en 15 sec de dénaturation à 95 °C, 30 sec de recuit à 60 °C, 2 min d’élongation à 68 °C et 5 min d’élongation à 68 °C. Les bibliothèques amplifiées ont été purifiées à l’aide de colonnes de spin MinElute avec le protocole standard fourni par le fabricant (Qiagen), et quantifiées pour le séquençage à l’aide d’une puce ADN 1000 du bioanalyseur Agilent 2100.
Pour l’individu sâme moderne, l’ADN total a été extrait au phénol-chloroforme et cisaillé physiquement en utilisant la fragmentation COVARIS. Une préparation de bibliothèque Illumina modifiée a été effectuée en utilisant la réparation de l’extrémité émoussée suivie de la queue A de l’extrémité 3′ et de la ligature d’adaptateurs fourchus. La PCR d’indexation a été suivie de l’excision de fragments allant de 500 à 600 pb à partir d’un gel d’agarose à 2%.
Capture & séquençage
Nous avons utilisé la procédure de capture en solution de la réf. 62 pour enrichir nos bibliothèques de fragments d’ADN chevauchant les 1 237 207 positions variables du génome humain4. Les séquences utilisées comme appât, attachées à des billes magnétiques, ont été mélangées à la matrice de l’échantillon d’ADN en solution, et laissées pour s’hybrider avec l’ADN cible pendant une incubation de 24 heures à 60 °C dans un four rotatif. 4 à 6 échantillons ont été regroupés dans des rapports de masse égaux pour un total de 2000 ng d’ADN. Les bibliothèques capturées ont été séquencées (75 pb à l’unité, plus des paires supplémentaires pour les trois bibliothèques non-UDG des individus de Levänluhta) sur une plateforme Illumina HiSeq 4000 à l’Institut Max Planck pour la science de l’histoire humaine à Jena. Sur les 13 échantillons de l’âge du fer originellement traités en Finlande, sept se sont révélés d’une qualité suffisante pour être utilisés dans des analyses en aval. Le génome sâme moderne a été séquencé dans sur un Genome Analyser II (8 voies, 125 bp paired-end) à l’Institut Max Planck d’anthropologie évolutive de Leipzig.
Traitement des lectures séquencées
Nous avons utilisé EAGER63 (version 1.92.50) pour traiter les lectures séquencées, en utilisant les paramètres par défaut (voir ci-dessous) pour les données de séquençage à une extrémité d’origine humaine, traitées par UDG-half, lors du traitement des bibliothèques UDG-half pour tous les individus. Plus précisément, AdapterRemoval a été utilisé pour éliminer les adaptateurs de séquençage de nos lectures, avec un chevauchement minimal de 1 pb, et en utilisant une qualité de base minimale de 20 et une longueur de séquence minimale de 30 pb. BWA aln (version 0.7.12-r1039, https://sourceforge.net/projects/bio-bwa/files)64 a été utilisé pour mapper les lectures à la séquence de référence humaine hs37d5, avec une longueur de graine (-l) de 32, un nombre maximal de différences (-n) de 0,01 tout en ne faisant aucun filtrage de qualité. La suppression des doublons a été effectuée en utilisant DeDup v0.12.1. Le calcul des dommages de désamination des bases terminales a été effectué à l’aide de mapDamage, en spécifiant une longueur (-l) de 100 pb (tableau supplémentaire 1).
Pour les analyses en aval, nous avons utilisé bamutils (version 1.0.13, https://github.com/statgen/bamUtil.git) TrimBam pour rogner deux bases au début et à la fin de toutes les lectures. Cette procédure élimine les positions qui sont affectées par la désamination, supprimant ainsi les erreurs de génotypage qui pourraient survenir en raison de dommages anciens à l’ADN.
Pour trois individus de Levänluhta dépassant le seuil de couverture de 1 % lors du dépistage préliminaire, nous avons utilisé les bibliothèques non traitées par l’UDG pour confirmer l’authenticité des données anciennes. Pour ces bibliothèques non traitées, deux tours de séquençage ont été effectués, qui ont été traités en utilisant EAGER avec les paramètres ci-dessus, mais en spécifiant un mode de traitement non-UDG et en définissant le type de séquençage correct entre les bibliothèques. Les lectures fusionnées ont été extraites des fichiers bam résultants, et fusionnées avec le fichier bam contenant les lectures du cycle de séquençage à extrémité unique à l’aide de samtools merge (version 1.3)65.
Le génome moderne Saami a été généré à l’aide d’Ibis pour l’appel de base et d’un script interne de rognage d’adaptateur. Les lectures résultantes ont ensuite été alignées sur le génome de référence humain hs37d5 en utilisant bwa 0.5.9-r16 (paramètres -e 20 -o 2 -n 0,01).
Génotypage
Nous avons utilisé un programme personnalisé (pileupCaller) pour génotyper les 15 individus anciens. Un fichier pileup a été généré à l’aide de samtools mpileup avec les paramètres -q 30 -Q 30 -B contenant uniquement les sites se chevauchant avec notre panel de capture. À partir de ce fichier, pour chaque individu et chaque SNP du panel 1240K, une lecture couvrant le SNP a été tirée au hasard, et un appel pseudo-haploïde a été effectué, c’est-à-dire que l’individu ancien a été supposé homozygote pour l’allèle sur la lecture tirée au hasard pour le SNP en question. PileupCaller est disponible à l’adresse https://github.com/stschiff/sequenceTools.git.
Pour les trois bibliothèques de Levänluhta qui n’ont pas subi de traitement UDG, nous avons uniquement génotypé les transversions, éliminant ainsi les artefacts des dommages post-mortem C- > T des analyses ultérieures.
Le génome shotgun de l’individu moderne Saami a été génotypé en utilisant GATK (version 1.3-25-g32cdef9) Unified Genotyper après réalignement des indels. Les appels de variants ont été filtrés pour les variants ayant un score de qualité supérieur à 30, et un script personnalisé a été utilisé pour convertir les variants au format EigenStrat.
Les données ont été fusionnées avec un grand ensemble de données composé de 3871 individus anciens et modernes génotypés sur les matrices Human Origins et/ou 1240K SNP, en utilisant mergeit.
Détermination du sexe
Pour déterminer le sexe génétique de chaque individu ancien, nous avons calculé la couverture sur les autosomes ainsi que sur chaque chromosome sexuel. Nous avons utilisé un script personnalisé (https://github.com/TCLamnidis/Sex.DetERRmine) pour le calcul de chaque couverture relative ainsi que leurs barres d’erreur associées (figure supplémentaire 1, note supplémentaire 3 pour plus d’informations sur le calcul des erreurs). Les femmes devraient avoir un taux x de 1 et un taux y de 0, tandis que les hommes devraient avoir un taux x et un taux y de 0,5 (réf. 49).
Authentification
Nous avons d’abord confirmé que le motif de désamination aux bases terminales des lectures d’ADN était caractéristique de l’ADN ancien (1.04-4,5 % pour les bibliothèques non-UDG, et 4,7-9,5 % pour les bibliothèques non-UDG, voir le tableau supplémentaire 1), en utilisant mapDamage (version 2.0.6)66. Nous avons effectué un certain nombre de tests différents pour garantir l’authenticité de nos données anciennes. Pour les individus mâles, nous avons étudié les polymorphismes sur le chromosome X27 en utilisant le progiciel ANGSD (version 0.910)67. Cela a révélé des estimations de contamination robustes pour 2 individus masculins de Bolshoy, et 1 individu masculin de Chalmny-Varre. Toutes ces estimations étaient inférieures à 1,6 % de contamination (tableau 1). Pour les individus femelles de ces deux sites, nous notons qu’ils sont projetés près des mâles dans l’espace PCA (Fig. 2a, Figure supplémentaire 3), ce qui suggère des effets limités de la contamination potentielle. En outre, nous avons généré un ensemble de données filtrées PMD pour tous les individus en utilisant pmdtools (version 0.60)30. Le filtrage PMD a été effectué en utilisant un génome de référence masqué pour toutes les positions du panel de capture 1240K, afin d’éviter les biais alléliques systématiques sur les positions SNP analysées. Nous avons fixé un seuil pmd de 3, ce qui, selon la publication originale30, élimine efficacement les contaminants modernes potentiels sur la base de l’absence de modifications de base compatibles avec la désamination.
Pour fournir une estimation plus quantitative de la contamination possible chez les femelles, nous avons utilisé le programme ContamMix (version 1.0-10)29 pour estimer la contamination mitochondriale. Nous avons extrait les reads mappant à la référence mitochondriale pour chacun des individus anciens en utilisant samtools (version 1.3)65. Nous avons ensuite généré une séquence consensus mitochondriale pour chacun des individus anciens en utilisant Geneious (version 10.0.9, http://www.geneious.com,68), et en appelant N tous les sites dont la couverture est inférieure à 5. Enfin, toutes les lectures mitochondriales ont été alignées sur leur séquence consensus respective, en utilisant bwa aln (version 0.7.12-r1039)64 avec un nombre maximal de différences dans la graine (-k) fixé à 5 et le nombre maximal de différences (-n) à 10, et bwa samse (version 0.7.12-r1039)64. Un alignement multiple de la séquence consensus et d’un ensemble de référence de 311 génomes mitochondriaux69 a été généré, en utilisant mafft (version v7.305b)70,71,72 avec le paramètre –auto. L’alignement de lecture, ainsi que l’alignement multiple de la séquence consensus et des 311 génomes mitochondriaux de référence ont ensuite été fournis à ContamMix. Nous rapportons ici le mode de contamination a posteriori, ainsi que les limites supérieures et inférieures de l’intervalle postérieur à 95 % (tableau 1).
Pour une authentification supplémentaire, nous avons exécuté ADMIXTURE28 supervisé (version 1.3.0) pour tous les échantillons en utilisant les six populations actuelles (Atayal, Français, Kalash, Karitiana, Mbuti et Papou) comme groupes génétiques définis, afin de localiser toute grande différence de regroupement génétique entre les individus d’un même site (figure supplémentaire 2). Nous avons testé la puissance de cette méthode pour détecter la contamination et nous avons constaté qu’elle peut détecter une contamination qui est lointainement liée aux ascendances présentes chez les individus testés déjà à des taux de 5-8%, mais elle n’a pas la puissance nécessaire pour identifier une contamination étroitement liée aux individus testés (voir la note supplémentaire 2). Nous n’avons pas observé de différences significatives (dans le cadre de notre résolution) dans les modèles d’ascendance entre les individus anciens d’un même site, à l’exception de Levänluhta, où l’échantillon individuel JK2065 semble dériver d’une ascendance différente. Nous lui avons donc attribué une étiquette de population distincte, Levänluhta_B dans cette étude.
Enfin, à l’aide de smartpca, nous avons projeté les ensembles de données filtrés par PMD et non filtrés sur le même ensemble de composantes principales construites sur les populations européennes modernes, afin de nous assurer que les individus anciens restent projetés dans les positions à peu près équivalentes, indépendamment du filtrage PMD. Cela a été possible pour tous les échantillons avec traitement UDG-demi, à l’exception des individus de Levänluhta, qui représentaient trop peu de dommages pour que le filtrage PMD puisse être appliqué. En ce qui concerne ce site, nous nous sommes donc appuyés sur les bibliothèques non-UDG (utilisant uniquement les transversions) qui ont été générées pour les trois individus utilisés dans l’analyse principale. Nous avons constaté que, dans les limites du bruit attendu en raison du faible nombre de SNP, tous les échantillons présentent une cohérence entre les ensembles de données filtrées et non filtrées, ce qui suggère un faible niveau de contamination dans tous les échantillons (figure supplémentaire 3a, b). Quatre individus supplémentaires de Levänluhta ont été exclus de l’analyse principale et de ce test d’authentification en raison de la faible couverture (< 15 000 SNP couverts) et de l’absence de bibliothèques non-UDG.
Statistiques F
Tous les programmes utilisés pour calculer les statistiques F dans cette étude se trouvent dans le cadre du paquet Admixtools (https://github.com/DReichLab/AdmixTools)2,42.
Nous avons utilisé qp3Pop (version 412) pour tous les calculs F3.
Les statistiques F4 ont été calculées à l’aide de qpDstat (version 711), et qpAdm (version 632)2 a été utilisé pour estimer les proportions de mélange en utilisant ce qui suit : Sources (Populations gauches) : Nganasan ; WHG ; EHG ; Yamnaya_Samara ; LBK_EN. Outgroups (populations de droite) : OG1 : Mbuti ; CHG ; Israel_Natufian ; Onge ; Villabruna ; Ami ; Mixe. OG2 : Mbuti ; CHG ; Onge ; Villabruna ; Ami ; Mixe. OG3 : Mixe ; CHG ; Israel_Natufian ; Villabruna ; Onge ; Ami. OG4 : Mbuti ; Israel_Natufian ; Onge ; Villabruna ; Ami ; Mixe. OG5 : Mbuti ; Samara_HG ; CHG ; Israel_Natufian ; Villabruna ; Ami.
Pour nous assurer que les ensembles d’outgroupes avaient suffisamment de puissance pour distinguer les ascendances présentes dans les sources, nous avons exécuté qpWave (version 410) en utilisant uniquement les sources comme populations de gauche et chaque ensemble d’outgroupes comme droits. Toutes ces exécutions de qpWave n’ont été cohérentes qu’avec le rang maximum, ce qui signifie que tous les ensembles d’outgroupes avaient suffisamment de puissance pour distinguer les cinq sources différentes. Tous les modèles qpWave et qpAdm ont été exécutés en utilisant l’option allsnps : YES. Lorsque les modèles de mélange à cinq voies fournis par qpAdm présentaient des valeurs p supérieures à 0,05, mais incluaient des proportions de mélange infaisables et que l’une des sources se voyait attribuer une proportion de mélange négative, nous avons exécuté à nouveau le modèle en excluant cette source. Pour chaque population test, si l’ensemble d’outgroupes OG1 n’a pas produit un modèle complet fonctionnel (p < 0,05), nous avons essayé d’autres ensembles d’outgroupes en retirant une population droite. Cela a donné les ensembles d’outgroupes OG2-4. Dans le cas de Levänluhta, plusieurs ensembles d’outgroupes ont produit des modèles fonctionnels, qui sont énumérés dans les données supplémentaires 4. Le modèle de mélange nécessitant le nombre minimum de sources tout en fournissant des proportions de mélange réalisables est toujours indiqué. Dans le cas de PWC de Suède, où aucun des ensembles de sous-groupes OG1-4 n’a produit un modèle fonctionnel, un ensemble révisé de populations correctes a été utilisé (OG5), qui inclut Samara_HG afin de fournir plus de puissance pour distinguer les ancêtres des chasseurs-cueilleurs. Nous avons préféré les modèles avec OG1-4 à OG5 en général, parce que OG5 contient plus de génomes anciens avec des biais potentiels dans les bonnes populations, ce qui provoque plus souvent l’échec des modèles pour les populations de test modernes. Les sources exclues dans les modèles minimaux ont été spécifiées comme N/A (Supplementary Data 4). Si Yamnaya ou EHG pouvaient être abandonnés (comme c’est le cas pour Levänluhta), nous montrons le modèle qui est plus cohérent avec les publications précédentes3,7,8,45 dans la figure 4, mais nous montrons les deux modèles dans les données supplémentaires 4.
Analyse en composantes principales
Nous avons utilisé smartpca (version #16000)73 (https://github.com/DReichLab/EIG) pour effectuer une analyse en composantes principales (ACP), en utilisant les paramètres lsqproject : YES et shrinkmode : YES.
Pour l’ACP eurasienne (Fig. 2a), les populations suivantes ont été utilisées pour construire les composantes principales : Abkhasian, Adygei, Albanais, Altaian, Ami, Arménien, Atayal, Avar.SG, Azeri_WGA, Balkar, Balochi, Basque, BedouinA, BedouinB, Belarusian, Borneo, Brahui, Bulgare, Buryat.SG, Cambodgien, Canariens, Tchétchène, Tchouvache, Croate, Chypriote, Tchèque, Dai, Daur, Dolgan, Druze, Anglais, Estonien, Even, Finnois, Français, Géorgien, Grec, Han, Hazara, Hezhen, Hongrois, Islandais, Iranien, Italien_Nord, Italien_Sud, Japonais, Juif_Ashkenazi, Juif_Georgien, Juif_Iranien, Juif_Irakien, Juif_Libyen, Juif_Marocain, Juif_Tunisien, Juif_Turc, Juif_Yéménite, Jordanien, Kalach, Kalmyk, Kinh, Coréen, Kumyk, Kurd_WGA, Kirghiz, Lahu, Libanais, Lezgin, Lituanien, Makrani, Mala, Maltais, Mansi, Miao, Mongola, Mordovian, Naxi, Nganasan, Nogai, North_Ossetian.DG, Norvégien, Orcadien, Oroqen, Palestinien, Pathan, Russe, Saami.DG, Saami_WGA, Sarde, Écossais, Selkup, Semende, She, Sherpa.DG, Sicilien, Espagnol, Espagnol_Nord, Syrien, Tadjik, Thaï, Tibétain.DG, Tu, Tubalar, Tujia, Turc, Turkmène, Tuvinian, Ukrainien, Ulchi, Uygur, Uzbek, Xibo, Yakut, Yi, Yukagir.
Pour l’ACP de l’Eurasie occidentale (figure supplémentaire 3a, b), les populations suivantes ont été utilisées pour construire les composantes principales : Abkhasian, Adygei, Albanais, Arménien, Balkar, Basque, BédouinA, BédouinB, Biélorusse, Bulgare, Canary_Islander, Tchétchène, Tchouvache, Croate, Chypriote, Tchèque, Druze, Anglais, Estonien, Finnois, Français, Géorgien, Grec, Hongrois, Islandais, Iranien, Italien_Nord, Italien_Sud, Juif_Ashkenazi, Juif_Georgien, Juif_Iranien, Juif_Irakien, Juif_Libyen, Juif_Marocain, Juif_Tunisien, Juif_Turc, Juif_Yéménite, Jordanien, Kumyk, Libanais, Lezgin, Lituanien, Maltais, Mordovien, Nord_Ossète, Norvégien, Orcadien, Palestinien, Polonais, Russe, Sarde, Saoudien, Écossais, Sicilien, Espagnol, Espagnol_Nord, Syrien, Turc, Ukrainien.
Analyse ADMIXTURE
ADMIXTURE28 a été exécuté avec la version 1.3.0, après exclusion des variants dont la fréquence de l’allèle mineur était de 0,01 et après élagage LD à l’aide de plink (version 1.90b3.29)74 avec une fenêtre de 200, un pas de 5 et un seuil R2 de 0,5 (https://www.genetics.ucla.edu/software/admixture/download.html). Cinq répétitions ont été effectuées pour chaque valeur K, les valeurs K étant comprises entre 2 et 15. Les populations utilisées étaient : Ami, Ami.DG, Arménien, Atayal, Atayal.DG, Balochi, Basque, BedouinB, Biélorusse, Brahmin_Tiwari, Brahui, Chuvash, Croate, Chypriote, Tchèque, Anglais, Estonien, Even, Finnois, Finnois.DG, Français, Grec, GujaratiB, Hadza, Han, Hongrois, Islandais, Kalash, Karitiana, Lituanien, Makrani, Mala, Mansi, Mansi.DG, Mari.SG, Mbuti, Mbuti.DG, Mixe, Mordovien, Nganasan, Norvégien, Onge, Orcadien, Papou, Pima, Russe, Saami.DG, ModernSaami, Sarde, Écossais, Selkup, Espagnol, Ukrainien, Ulchi, Yoruba, ALPC_Hongrie_MN, Baalberge_MN, Baltic_BA, Baltic_CCC, Baltic_CWC, Baltic_LN, BolshoyOleniOstrov, Bu_kk_Culture_Hongrie_MN, ChalmnyVarre, CHG, EHG, Esperstedt_MN, Ganj_Dareh_Iran_néolithique, Hongrie_MN, Hongrie_néolithique, Iran_Chalcolithique, JK2065, Koros_Hongrie_EN, Kunda, Lettonie_HG3, Lettonie_MN1, LBK_EN, LBK_Hongrie_EN, Levanluhta, Narva, PWC_Suède_NHG.SG, Scandinavie_LNBA, SHG, Suède_HG.SG, TRB, Ukraine_HG1, Ukraine_N1, WHG, Yamnaya_Samara.
Nous constatons que K = 11 entraîne l’erreur de validation croisée la plus faible, comme le montre la figure supplémentaire 4b.
Haplotypage du chromosome Y
Nous avons assigné les hommes anciens à des haplogroupes Y en utilisant le programme yHaplo (https://github.com/23andMe/yhaplo)75. En bref, ce programme fournit une recherche automatisée dans l’arbre des haplogroupes Y (tel que fourni dans yHaplo, tel qu’accédé à partir de l’ISOGG le 04 Jan 2016) de la racine à la branche en aval en fonction de la présence d’allèles dérivés et attribue l’haplogroupe le plus en aval avec des allèles dérivés. Pour environ 15 000 SNP du chromosome Y présents à la fois dans notre panel de capture et dans deux ensembles de données publiés76,77, nous avons échantillonné au hasard une seule base et l’avons utilisée comme génotype haploïde. Nous avons utilisé un script personnalisé pour convertir les génotypes EigenStrat au format yHaplo. Nous rapportons l’haplogroupe assigné le plus en aval fourni par le programme (Tableau 1). Nous avons également vérifié manuellement le statut dérivé et l’absence de mutations définissant l’haplogroupe désigné, car des informations manquantes pourraient entraîner un arrêt prématuré de sa recherche automatisée.
Haplotypage mitochondrial
Nous avons importé les lectures mitochondriales rognées pour chaque individu avec une qualité de cartographie >30 dans Geneious (version 10.0.9, https://www.geneious.com)68 et réassemblé ces lectures au génome de référence RSRS78, en utilisant le mappeur Geneious, avec une sensibilité moyenne et 5 itérations. Nous avons utilisé l’appelant automatique de variants intégré à Geneious pour trouver les polymorphismes mitochondriaux avec une couverture minimale de 3 et une fréquence minimale de variants de 0,67. Les variants résultants ont été exportés vers Excel et comparés manuellement aux SNP rapportés dans la phylogénie ADNmt en ligne (arbre ADNmt Build 17, 18 févr. 2016, http://www.phylotree.org/). Les positions nucléotidiques 309.1 C(C), 315.1C, les indels AC à 515-522, 16182C, 16183C, 16193.1C(C) et 16519 ont été masqués et non inclus dans nos appels d’haplotype.
SNP phénotypiques
Nous avons utilisé samtools mpileup (version 1.3)65, en filtrant pour une qualité de carte (-Q) et de base (-q) de 30, en désactivant la qualité d’alignement par base (-B), sur les fichiers bam rognés, pour générer un pileup de lectures correspondant à un ensemble de 43 SNP phénotypiques4,40,41,79 qui font partie de notre panel de capture du génome. Un script python personnalisé a été utilisé pour analyser le pileup dans un tableau contenant le nombre de lectures soutenant chaque allèle (Données supplémentaires 2).
Disponibilité du code
Tous les logiciels décrits pour la première fois dans cette étude sont librement disponibles dans des dépôts en ligne. Sex.DetERRmine : https://github.com/TCLamnidis/Sex.DetERRmine
Les génotypes contaminés : https://github.com/TCLamnidis/ContaminateGenotypes
.