Sampling
Die schriftliche, informierte Zustimmung wurde von dem Saami Individuum eingeholt, dessen Genom in dieser Studie analysiert wurde, Die Studie wurde von der Ethikkommission des Krankenhausbezirks Helsinki und Uusimaa (Entscheidung 329/13/03/00/2013) und der Ethikkommission der Universität Leipzig (Genehmigungsreferenznummer 398-13-16122013) genehmigt.
Die Entnahme und Extraktion antiker DNA erfordert ein strenges Verfahren, um Kontaminationen durch zeitgenössisches genetisches Material zu vermeiden. Für die 13 eisenzeitlichen Individuen aus Finnland, die uns zur Verfügung standen, fand die Probenahme in einer Reinraumanlage für antike DNA-Arbeiten am Institut für Archäologische Wissenschaften in Tübingen statt. Der vorläufige Arbeitsablauf beinhaltete die Dokumentation, das Fotografieren und die Lagerung der Proben in individuellen, ID-codierten Plastikröhrchen und Plastikbeuteln. Als Ergebnis einer frühen Pilotstudie waren die von uns verwendeten Zahnproben fragmentiert, und ein Teil des Dentins wurde entfernt. Das verbleibende Dentin wurde durch vorsichtiges Trennen vom Zahnschmelz mit einem Zahnarztbohrer und gekühlten Diamantbohrköpfen, die mit einer Geschwindigkeit von weniger als 15 Umdrehungen pro Minute rotierten, entnommen, um eine mögliche hitzebedingte Schädigung der antiken DNA zu vermeiden.
Für die Proben von den Fundorten Bolshoy und Chalmny Varre verwendeten wir übrig gebliebenes Zahnpulver, das ursprünglich am Institut für Anthropologie der Universität Mainz zu Replikationszwecken aufbereitet wurde, wie in Ref. 24 beschrieben. Kurz gesagt beinhalteten die Probenvorbereitungsschritte eine UV-Bestrahlung für 30-45 min, gefolgt von einem sanften Abwischen der Oberfläche mit verdünnten kommerziellen Bleichmitteln. Die Zähne wurden dann mit Aluminiumoxid-Schleifmittel (Harnisch & Rieth) sandgestrahlt und mit einer Rührwerksmühle (Retsch) zu feinem Pulver gemahlen.
Kalibrierung des Radiokarbondatums
Wir kalibrierten das Radiokarbondatum von Bolshoy, das in Ref. 24,55 als 3473 ± 42 Jahre BP angegeben wird, unter Verwendung von Intcal1356 als Kalibrierungskurve, mit OxCal 4.357.
DNA-Extraktion und Bibliotheksvorbereitung
DNA aus den sechs Bolshoy- und den zwei Chalmny Varre-Proben wurde in den Antiken-DNA-Einrichtungen des Max-Planck-Instituts für die Erforschung der Menschheitsgeschichte (MPI-SHH) in Jena, Deutschland, extrahiert. Die Extraktion für die Levänluhta-Proben wurde in ähnlicher Weise in den Reinraumanlagen des Instituts für Archäologische Wissenschaften in Tübingen durchgeführt. Für jede Probe wurden ~ 50 mg Dentinpulver für ein Extraktionsverfahren verwendet, das speziell für die Gewinnung antiker DNA entwickelt wurde58. Extraktionspuffer mit 0,45 M EDTA, pH 8,0 (Life Technologies) und 0,25 mg/ml Proteinase K (Sigma-Aldrich) wurde zum Knochenpulver gegeben und über Nacht bei 37 °C unter Rotation inkubiert. Der Überstand wurde durch Zentrifugation (14.000 rpm) vom Pellet des Knochenpulvers getrennt. Ein Bindungspuffer, bestehend aus 5 M GuHCL (Sigma Aldrich) und 40 % Isopropanol (Merck), wurde zusammen mit 400 μl 1 M Natriumacetat (pH 5,5) zum Überstand gegeben und die Lösung durch Schleudern durch eine Reinigungssäule, die an einem High Pure Extender Assembly-Trichter befestigt war, gereinigt (8 Min. bei 1500 U/min, mit langsamer Beschleunigung). Die Säule wurde dann zur Maximierung der Ausbeute 1-2 Mal in ein Sammelröhrchen geschleudert (1 min bei 14.000 U/min). Es folgten zwei weitere Waschschritte mit 450 μl Waschpuffer (High Pure Viral Nucleic Acid Large Volume Kit) und zwei Trockenschleuderschritte mit 1 min Zentrifugation bei 14.000 U/min. Das endgültige Gesamtvolumen von 100 μl Eluat wurde durch zwei separate Elutionsrunden von 50 μl TET (10 mM Tris-HCL, 1 mM EDTA pH 8,0, 0,1% Tween20) erreicht, die jeweils für 1 min bei 14.000 rpm in ein frisches Eppendorf 1,5 ml Röhrchen geschleudert wurden. Negativkontrollen (Puffer statt Probe) wurden parallel in einem Verhältnis von 1 Kontrolle pro 7 Proben prozessiert.
Von den 100 μl Extrakt wurden 20 μl verwendet, um die Proben-DNA als doppelsträngige Bibliothek zu immortalisieren. Die Prozedur beinhaltete einen Blunt-End-Reparatur-, Adapter-Ligations- und Adapter-Fill-In-Schritt, wie von Meyer und Kircher59 beschrieben. Während des Blunt-End-Reparaturschritts wurde eine Mischung aus 0,4 U/μl T4 PNK (Polynukleotidkinase) und 0,024 U/μl T4 DNA-Polymerase, 1× NEB-Puffer 2 (NEB), 100 μM dNTP-Mix (Thermo Scientific), 1 mM ATP (NEB) und 0.8 mg/ml BSA (NEB) wurden der Template-DNA zugesetzt, gefolgt von der Inkubation im Thermocycler (15 min 15 °C, 15 min 25 °C) und der Aufreinigung mit einem MinElute-Kit (QIAGEN). Das Produkt wurde in 18 μl TET-Puffer eluiert. Der Adapter-Ligationsschritt umfasste eine Mischung aus 1× Quick Ligase Buffer (NEB), 250 nM Illumina Adaptern (Sigma-Aldrich) und 0,125 U/μl Quick Ligase (NEB), die dem 18 μl Eluat zugegeben wurde, gefolgt von einer 20-minütigen Inkubation und einem zweiten Reinigungsschritt mit MinElute-Säulen, diesmal in 20 μl Eluat. Für den Auffüllschritt wurde eine Mischung aus 0,4 U/μl Bst-Polymerase und 125 μM dNTP-Mix zugegeben und die Mischung anschließend im Thermocycler inkubiert (30 min 37 °C, 10 min 80 °C). Bibliotheken ohne Uracil-DNA-Glykosylase (UDG)-Behandlung wurden für alle 13 Extrakte aus Levänluhta hergestellt. Zusätzlich wurden UDG-half behandelte Bibliotheken für sieben der ursprünglichen 13 Extrakte aus Levänluhta und für alle Bolshoy- und Chalmny Varre-Extrakte hergestellt. Um die UDG-half-Behandlung einzuführen, wurde ein erster Schritt in die Bibliotheksvorbereitung aufgenommen, bei dem 250 U USER-Enzym (NEB) in die 20 μl Extrakt gegeben wurden, gefolgt von einer Inkubation bei 37 °C für 30 min und dann 12 °C für 1 min. Darauf folgte wiederum die Zugabe von 200 U UGI (Uracil-Glykosylase-Inhibitor, von NEB) und eine weitere identische Inkubation, um die enzymatische Exzision der deaminierten Stellen zu stoppen, wie in60 beschrieben. Für jede Bibliothek wurde ein einzigartiges Paar von acht bp langen Indizes unter Verwendung einer Pfu Turbo Cx Hotstart DNA Polymerase und eines Thermocycling-Programms mit folgendem Temperaturprofil eingearbeitet: anfängliche Denaturierung (98 °C für 30 Sek.), Zyklus von Denaturierung/Verhärtung/Verlängerung (98 °C für 10 Sek./ 60 °C für 20 Sek./ 72 °C für 20 Sek.) und abschließende Verlängerung bei 72 °C für 10 Min.61. Knochenpulver von einem Höhlenbären wurde parallel als Positivkontrolle bearbeitet. Negativkontrollen sowohl für die Extraktions- als auch für die Bibliotheksvorbereitung wurden während des gesamten Arbeitsablaufs neben den Proben aufbewahrt.
Die Effizienz des Experiments wurde durch die Quantifizierung der Konzentration der Bibliotheken mittels qPCR (Roche) unter Verwendung von Aliquots der Bibliotheken vor und nach der Indizierung sichergestellt. Die molekulare Kopienzahl in den vorindizierten Bibliotheken reichte von ~ 10E8 bis ~ 10E9 Kopien/μl, was auf eine erfolgreiche Bibliothekspräparation hinweist, während die indizierten Bibliotheken von ~ 10E10 bis ~ 10E12 Kopien/μl reichten, was auf eine zulässige Indizierungseffizienz hindeutet. Die Negativkontrollen zeigten 4-5 Größenordnungen niedrigere Konzentrationen als die Proben, was auf eine geringe Kontamination durch die Verarbeitungsschritte im Labor hindeutet.
Die Bibliotheken wurden mit PCR amplifiziert, und zwar für die Anzahl der Zyklen, die den Konzentrationen der indizierten Bibliotheken entsprachen, unter Verwendung von AccuPrime Pfx Polymerase (5 μl Bibliotheks-Template, 2 U AccuPrime Pfx DNA Polymerase von Invitrogen, 1 U vorgefertigter 10× PCR-Mastermix und 0.3 μM der Primer IS5 und IS6, für jeweils 100 μl Reaktion) mit einem thermischen Profil von 2 min Denaturierung bei 95 °C, 3-9 Zyklen bestehend aus 15 s Denaturierung bei 95 °C, 30 s Annealing bei 60 °C, 2 min Elongation bei 68 °C und 5 min Elongation bei 68 °C. Die amplifizierten Bibliotheken wurden mit MinElute-Spinsäulen nach dem Standardprotokoll des Herstellers (Qiagen) aufgereinigt und für die Sequenzierung mit einem Agilent 2100 Bioanalyzer DNA 1000 Chip quantifiziert.
Für das moderne Saami-Individuum wurde die Gesamt-DNA mit Phenol-Chloroform extrahiert und mit COVARIS-Fragmentierung physikalisch geschert. Eine modifizierte Illumina-Bibliothekspräparation wurde mit Blunt-End-Repair, gefolgt von A-Tailing des 3′-Endes und Ligation von Forked-Adaptern durchgeführt. Nach der Indizierungs-PCR erfolgte die Exzision von Fragmenten im Bereich von 500 bis 600 bp aus einem 2%igen Agarosegel.
Capture & Sequenzierung
Wir verwendeten die In-Solution-Capture-Prozedur aus Ref. 62, um unsere Bibliotheken mit DNA-Fragmenten anzureichern, die sich mit 1.237.207 variablen Positionen im menschlichen Genom4 überlappen. Die als Köder verwendeten Sequenzen, die an Magnetbeads gebunden waren, wurden mit der DNA-Probenvorlage in Lösung gemischt und während einer 24-stündigen Inkubation bei 60 °C in einem Rotationsofen mit der Ziel-DNA hybridisiert. 4-6 Proben wurden in gleichen Massenverhältnissen zu insgesamt 2000 ng DNA gepoolt. Die erfassten Bibliotheken wurden auf einer Illumina HiSeq 4000-Plattform am Max-Planck-Institut für die Erforschung der Menschheitsgeschichte in Jena sequenziert (75 bp single-end, plus zusätzlich paired-end für die drei Nicht-UDG-Bibliotheken der Levänluhta-Individuen). Von den 13 ursprünglich prozessierten eisenzeitlichen Proben aus Finnland erwiesen sich sieben als von ausreichender Qualität, um in nachgeschalteten Analysen verwendet zu werden. Das moderne Saami-Genom wurde auf einem Genome Analyser II (8 Lanes, 125 bp paired-end) am Max-Planck-Institut für evolutionäre Anthropologie in Leipzig sequenziert.
Processing der sequenzierten Reads
Wir verwendeten EAGER63 (Version 1.92.50), um die sequenzierten Reads zu verarbeiten, wobei wir die Standardparameter (siehe unten) für von Menschen stammende, mit UDG-half behandelte Single-End-Sequenzierdaten verwendeten, als wir die UDG-half-Bibliotheken für alle Individuen verarbeiteten. Insbesondere wurde AdapterRemoval verwendet, um die Sequenzierungsadapter aus unseren Reads zu entfernen, mit einer minimalen Überlappung von 1 bp und unter Verwendung einer minimalen Basenqualität von 20 und einer minimalen Sequenzlänge von 30 bp. BWA aln (Version 0.7.12-r1039, https://sourceforge.net/projects/bio-bwa/files)64 wurde verwendet, um die Reads auf die menschliche Referenzsequenz hs37d5 zu mappen, mit einer Seed-Länge (-l) von 32, einer maximalen Anzahl von Unterschieden (-n) von 0,01 und ohne Qualitätsfilterung. Die Entfernung von Duplikaten wurde mit DeDup v0.12.1 durchgeführt. Die Berechnung des terminalen Basendeaminationsschadens wurde mit mapDamage durchgeführt, wobei eine Länge (-l) von 100 bp angegeben wurde (ergänzende Tabelle 1).
Für die Downstream-Analysen verwendeten wir bamutils (Version 1.0.13, https://github.com/statgen/bamUtil.git) TrimBam, um zwei Basen am Anfang und Ende aller Reads zu trimmen. Diese Prozedur eliminiert die Positionen, die von der Desaminierung betroffen sind, und beseitigt so Genotypisierungsfehler, die durch alte DNA-Schäden entstehen könnten.
Für drei Levänluhta-Individuen, die im vorläufigen Screening die Schwellenabdeckung von 1 % überschritten, verwendeten wir die nicht-UDG-behandelten Bibliotheken, um die Authentizität der alten Daten zu bestätigen. Für diese unbehandelten Bibliotheken wurden zwei Sequenzierungsrunden durchgeführt, die mit EAGER mit den oben genannten Parametern prozessiert wurden, wobei jedoch ein nicht-UDG-Behandlungsmodus angegeben und der richtige Sequenzierungstyp zwischen den Bibliotheken eingestellt wurde. Die fusionierten Reads wurden aus den resultierenden bam-Dateien extrahiert und mit der bam-Datei, die die Reads aus dem Single-End-Sequenzierungslauf enthielt, unter Verwendung von samtools merge (Version 1.3)65 zusammengeführt.
Das moderne Saami-Genom wurde mit Ibis für das Base-Calling und einem hauseigenen Adapter-Trimming-Skript erzeugt. Die resultierenden Reads wurden dann mit bwa 0.5.9-r16 (Parameter -e 20 -o 2 -n 0.01) an das menschliche Referenzgenom hs37d5 angeglichen.
Genotypisierung
Wir verwendeten ein eigenes Programm (pileupCaller), um die 15 alten Individuen zu genotypisieren. Eine Pileup-Datei wurde mit samtools mpileup mit den Parametern -q 30 -Q 30 -B generiert, die nur Sites enthielt, die sich mit unserem Capture-Panel überschnitten. Aus dieser Datei wurde für jedes Individuum und jeden SNP auf dem 1240K-Panel ein Read, der den SNP abdeckt, zufällig gezogen und ein pseudohaploider Call durchgeführt, d. h. es wurde angenommen, dass das alte Individuum homozygot für das Allel auf dem zufällig gezogenen Read für den betreffenden SNP ist. PileupCaller ist verfügbar unter https://github.com/stschiff/sequenceTools.git.
Für die drei Levänluhta-Bibliotheken, die keiner UDG-Behandlung unterzogen wurden, haben wir nur Transversionen genotypisiert, um Artefakte von post-mortem C- > T-Schäden aus den weiteren Analysen zu eliminieren.
Das Shotgun-Genom des modernen Saami-Individuums wurde mit GATK (Version 1.3-25-g32cdef9) Unified Genotyper nach Indel-Realignment genotypisiert. Die Varianten-Calls wurden nach Varianten mit einem Qualitätsscore über 30 gefiltert, und ein benutzerdefiniertes Skript wurde verwendet, um die Varianten in das EigenStrat-Format zu konvertieren.
Die Daten wurden mit einem großen Datensatz fusioniert, der aus 3871 antiken und modernen Individuen besteht, die auf den Human Origins und/oder 1240K SNP Arrays genotypisiert wurden, unter Verwendung von mergeit.
Geschlechtsbestimmung
Um das genetische Geschlecht eines jeden antiken Individuums zu bestimmen, berechneten wir die Abdeckung auf den Autosomen sowie auf jedem Geschlechtschromosom. Wir verwendeten ein eigenes Skript (https://github.com/TCLamnidis/Sex.DetERRmine) für die Berechnung jeder relativen Abdeckung sowie der zugehörigen Fehlerbalken (Ergänzende Abbildung 1, Ergänzende Anmerkung 3 für weitere Informationen zur Fehlerberechnung). Bei Frauen wird eine x-Rate von 1 und eine y-Rate von 0 erwartet, während bei Männern sowohl eine x- als auch eine y-Rate von 0,5 erwartet wird (Ref. 49).
Authentifizierung
Wir bestätigten zunächst, dass das Deaminierungsmuster an den terminalen Basen der DNA-Reads charakteristisch für alte DNA ist (1.04-4,5 % für Nicht-UDG-Bibliotheken und 4,7-9,5 % für Nicht-UDG-Bibliotheken, siehe ergänzende Tabelle 1), indem wir mapDamage (Version 2.0.6)66 verwendeten. Wir führten eine Reihe verschiedener Tests durch, um die Authentizität unserer alten Daten sicherzustellen. Für männliche Individuen untersuchten wir Polymorphismen auf dem X-Chromosom27 mit dem ANGSD-Softwarepaket (Version 0.910)67. Dies ergab robuste Kontaminationsschätzungen für 2 männliche Bolshoy-Individuen und 1 männliches Chalmny-Varre-Individuum. Diese lagen alle unter 1,6 % Kontamination (Tabelle 1). Für die weiblichen Individuen von diesen beiden Standorten stellen wir fest, dass sie im PCA-Raum nahe an die männlichen Individuen projiziert werden (Abb. 2a, ergänzende Abbildung 3), was auf begrenzte Auswirkungen einer möglichen Kontamination hindeutet. Zusätzlich generierten wir einen PMD-gefilterten Datensatz für alle Individuen mit pmdtools (Version 0.60)30. Die PMD-Filterung erfolgte unter Verwendung eines Referenzgenoms, das für alle Positionen des 1240K-Capture-Panels maskiert war, um systematische allelische Verzerrungen an den analysierten SNP-Positionen zu vermeiden. Wir setzten einen pmd-Schwellenwert von 3, der gemäß der Originalveröffentlichung30 potenzielle moderne Kontaminationen aufgrund des Fehlens von Basenmodifikationen, die mit Deaminierung konsistent sind, effektiv eliminiert.
Um eine quantitativere Abschätzung möglicher Kontaminationen bei Frauen zu erhalten, verwendeten wir das Programm ContamMix (Version 1.0-10)29 zur Abschätzung der mitochondrialen Kontamination. Wir extrahierten Reads, die auf die mitochondriale Referenz abgebildet wurden, für jedes der alten Individuen mit samtools (Version 1.3)65. Dann generierten wir eine mitochondriale Konsenssequenz für jedes der alten Individuen mit Geneious (Version 10.0.9, http://www.geneious.com,68) und riefen N für alle Stellen mit einer Abdeckung von weniger als 5 auf. Abschließend wurden alle mitochondrialen Reads mit Hilfe von bwa aln (Version 0.7.12-r1039)64 mit einer maximalen Anzahl von Unterschieden im Seed (-k) auf 5 und einer maximalen Anzahl von Unterschieden (-n) auf 10 an ihrer jeweiligen Konsensussequenz ausgerichtet. Ein multiples Alignment der Konsensussequenz und eines Referenzsatzes von 311 mitochondrialen Genomen69 wurde mit mafft (Version v7.305b)70,71,72 mit dem Parameter –auto erstellt. Das Read-Alignment sowie das multiple Alignment der Konsensussequenz und der 311 mitochondrialen Referenzgenome wurden dann ContamMix zur Verfügung gestellt. Wir berichten hier über den a posteriori Modus der Kontamination, zusammen mit den oberen und unteren Grenzen des 95%-Posterior-Intervalls (Tabelle 1).
Für eine zusätzliche Authentifizierung führten wir ADMIXTURE28 (Version 1.3.0) für alle Proben unter Aufsicht durch, wobei wir die sechs heutigen Populationen (Atayal, French, Kalash, Karitiana, Mbuti und Papuan) als definierte genetische Cluster verwendeten, um große Unterschiede im genetischen Clustering zwischen Individuen vom gleichen Standort zu lokalisieren (ergänzende Abbildung 2). Wir testeten die Leistungsfähigkeit dieser Methode zur Erkennung von Kontaminationen und stellten fest, dass sie Kontaminationen, die mit den in den Testindividuen vorhandenen Abstammungen entfernt verwandt sind, bereits mit Raten von 5-8% erkennen kann, aber nicht die Leistungsfähigkeit besitzt, um Kontaminationen zu identifizieren, die eng mit den Testindividuen verwandt sind (siehe Anhang 2). Wir beobachteten keine signifikanten Unterschiede (innerhalb unserer Auflösung) in den Abstammungsmustern zwischen den alten Individuen vom gleichen Standort, mit Ausnahme von Levänluhta, wo die Einzelprobe JK2065 aus einer anderen Abstammung zu stammen scheint. Daher haben wir ihr eine separate Populationsbezeichnung zugewiesen, in dieser Studie Levänluhta_B.
Schließlich haben wir mit smartpca die PMD-gefilterten und nicht-gefilterten Datensätze auf denselben Satz von Hauptkomponenten projiziert, die für moderne europäische Populationen konstruiert wurden, um sicherzustellen, dass die alten Individuen unabhängig von der PMD-Filterung in den ungefähr gleichen Positionen projiziert werden. Dies war für alle Proben mit UDG-Halbwertbehandlung möglich, mit Ausnahme der Individuen aus Levänluhta, die für die Anwendung der PMD-Filterung zu geringe Schäden aufwiesen. Für diesen Standort stützten wir uns daher auf die Nicht-UDG-Bibliotheken (nur mit Transversionen), die für die drei in der Hauptanalyse verwendeten Individuen erzeugt wurden. Wir fanden heraus, dass innerhalb des erwarteten Rauschens aufgrund einer geringen Anzahl von SNPs alle Proben Konsistenz zwischen den gefilterten und nicht gefilterten Datensätzen zeigen, was auf eine geringe Kontamination in allen Proben hindeutet (ergänzende Abbildung 3a, b). Vier zusätzliche Individuen aus Levänluhta wurden wegen geringer Abdeckung (< 15.000 abgedeckte SNPs) und fehlender Non-UDG-Bibliotheken von der Hauptanalyse und von diesem Authentifizierungstest ausgeschlossen.
F-Statistik
Alle Programme, die zur Berechnung der F-Statistiken in dieser Studie verwendet wurden, sind als Teil des Admixtools-Pakets (https://github.com/DReichLab/AdmixTools)2,42 zu finden.
Wir verwendeten qp3Pop (Version 412) für alle F3-Berechnungen.
F4-Statistiken wurden mit qpDstat (Version 711) berechnet, und qpAdm (Version 632)2 wurde zur Schätzung von Mischungsverhältnissen verwendet: Quellen (linke Populationen): Nganasan; WHG; EHG; Yamnaya_Samara; LBK_EN. Outgroups (Rechte Populationen): 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.
Um sicherzugehen, dass die Outgroup-Sets genug Power haben, um die in den Quellen vorhandenen Abstammungen zu unterscheiden, haben wir qpWave (Version 410) laufen lassen, wobei wir nur die Quellen als linke Populationen und jedes Outgroup-Set als rechte verwendet haben. Alle diese qpWave-Läufe waren nur mit maximalem Rang konsistent, d. h. alle Outgroup-Sets hatten genug Power, um zwischen den fünf verschiedenen Quellen zu unterscheiden. Alle qpWave- und qpAdm-Modelle wurden mit der Option allsnps ausgeführt: YES. Wenn die von qpAdm bereitgestellten fünfseitigen Vermischungsmodelle p-Werte über 0,05 hatten, aber nicht durchführbare Mischungsanteile enthielten und einer der Quellen ein negativer Mischungsanteil zugewiesen wurde, führten wir das Modell erneut aus, wobei diese Quelle ausgeschlossen wurde. Wenn das Outgroup-Set OG1 für jede Testpopulation kein funktionierendes vollständiges Modell ergab (p < 0,05), haben wir alternative Outgroup-Sets ausprobiert, bei denen eine rechte Population entfernt wurde. Dies führte zu den Outgroup-Sets OG2-4. Im Fall von Levänluhta führten mehrere Outgroup-Sets zu funktionierenden Modellen, die in den ergänzenden Daten 4 aufgeführt sind. Es wird immer das Vermischungsmodell gezeigt, das die geringste Anzahl von Quellen benötigt und dennoch praktikable Vermischungsanteile liefert. Im Fall von PWC aus Schweden, wo keines der Outgroup-Sets OG1-4 ein funktionierendes Modell lieferte, wurde ein überarbeitetes Set von richtigen Populationen verwendet (OG5), das Samara_HG einschließt, um mehr Macht zur Unterscheidung von Jäger- und Sammler-Vorfahren zu haben. Wir haben Modelle mit OG1-4 gegenüber OG5 im Allgemeinen bevorzugt, weil OG5 mehr alte Genome mit potenziellen Verzerrungen in den richtigen Populationen enthält, was häufiger zum Scheitern von Modellen für moderne Testpopulationen führt. Die ausgeschlossenen Quellen in den Minimalmodellen wurden als N/A angegeben (Supplementary Data 4). Wenn entweder Yamnaya oder EHG ausgeschlossen werden konnten (wie es bei Levänluhta der Fall ist), zeigen wir in Abb. 4 das Modell, das mehr mit früheren Publikationen3,7,8,45 übereinstimmt, zeigen aber beide Modelle in Supplementary Data 4.
Hauptkomponentenanalyse
Wir benutzten smartpca (Version #16000)73 (https://github.com/DReichLab/EIG), um die Hauptkomponentenanalyse (PCA) durchzuführen, wobei wir die Einstellungen lsqproject: YES und shrinkmode: YES.
Für die eurasische PCA (Abb. 2a) wurden die folgenden Populationen zur Konstruktion der Hauptkomponenten verwendet: Abchasisch, Adygei, Albanisch, Altaisch, Ami, Armenisch, Atayal, Avar.SG, Azeri_WGA, Balkar, Belutschi, Baskisch, BeduinA, BeduinB, Weißrussisch, Borneo, Brahui, Bulgarisch, Burjat.SG, Kambodschanisch, Kanarische_Insulaner, Tschetschenisch, Tschuwaschisch, Kroatisch, Zypriotisch, Tschechisch, Dai, Daur, Dolgan, Drusen, Englisch, Estnisch, Even, Finnisch, Französisch, Georgisch, Griechisch, Han, Hazara, Hezhen, Ungarisch, Isländisch, Iranisch, Italienisch_Nord, Italienisch_Süd, Japanisch, Jew_Ashkenazi, Jew_Georgian, Jude_Iraner, Jude_Iraker, Jude_Libyer, Jude_Marokkaner, Jude_Tunesier, Jude_Türke, Jude_Jemenit, Jordanier, Kalash, Kalmyk, Kinh, Koreaner, Kumyk, Kurd_WGA, Kirgisisch, Lahu, Libanesisch, Lezgin, Litauisch, Makrani, Mala, Maltesisch, Mansi, Miao, Mongola, Mordovian, Naxi, Nganasan, Nogai, North_Ossetian.DG, Norwegisch, Orcadian, Oroqen, Palästinensisch, Pathan, Russisch, Saami.DG, Saami_WGA, Sardisch, Saudi, Schottisch, Selkup, Semende, She, Sherpa.DG, Sizilianisch, Spanisch, Spanisch_Nord, Syrisch, Tadschikisch, Thai, Tibetisch.DG, Tu, Tubalar, Tujia, Türkisch, Turkmenisch, Tuvinisch, Ukrainisch, Ulchi, Uygur, Usbekisch, Xibo, Yakut, Yi, Yukagir.
Für die westeurasische PCA (ergänzende Abbildung 3a, b) wurden die folgenden Populationen zur Konstruktion von Hauptkomponenten verwendet: Abchasisch, Adygei, Albanisch, Armenisch, Balkar, Baskisch, BeduineA, BeduineB, Weißrussisch, Bulgarisch, Kanarischer_Islander, Tschetschenisch, Tschuwaschisch, Kroatisch, Zypriotisch, Tschechisch, Drusen, Englisch, Estnisch, Finnisch, Französisch, Georgisch, Griechisch, Ungarisch, Isländisch, Iranisch, Italienisch_Nord, Italienisch_Süd, Jude_Ashkenazi, Jude_Georgisch, Jew_Iranian, Jew_Iraqi, Jew_Libyan, Jew_Moroccan, Jew_Tunisian, Jew_Turkish, Jew_Yemenite, Jordanian, Kumyk, Lezgin, Lebanese, Lithuanian, Maltese, Mordovian, North_Ossetian, Norwegian, Orcadian, Palestinian, Polish, Russian, Sardinian, Saudi, Scottish, Sicilian, Spanish, Spanish_North, Syrian, Turkish, Ukrainian.
ADMIXTURE-Analyse
ADMIXTURE28 wurde mit Version 1.3.0 durchgeführt, nach Ausschluss von Varianten mit einer Minor-Allel-Häufigkeit von 0,01 und nach LD-Pruning mit plink (Version 1.90b3.29)74 mit einer Fenstergröße von 200, einer Schrittgröße von 5 und einem R2-Schwellenwert von 0,5 (https://www.genetics.ucla.edu/software/admixture/download.html). Für jeden K-Wert wurden fünf Replikate durchgeführt, wobei die K-Werte zwischen 2 und 15 lagen. Die verwendeten Populationen waren: Ami, Ami.DG, Armenisch, Atayal, Atayal.DG, Belutschi, Baskisch, BeduinB, Weißrussisch, Brahmin_Tiwari, Brahui, Tschuwaschisch, Kroatisch, Zypriotisch, Tschechisch, Englisch, Estnisch, Even, Finnisch, Finnisch.DG, Französisch, Griechisch, GujaratiB, Hadza, Han, Ungarisch, Isländisch, Kalash, Karitiana, Litauisch, Makrani, Mala, Mansi, Mansi.DG, Mari.SG, Mbuti, Mbuti.DG, Mixe, Mordovian, Nganasan, Norwegisch, Onge, Orcadian, Papuan, Pima, Russisch, Saami.DG, ModernSaami, Sardisch, Schottisch, Selkup, Spanisch, Ukrainisch, Ulchi, Yoruba, ALPC_Ungarn_MN, Baalberge_MN, Baltic_BA, Baltic_CCC, Baltic_CWC, Baltic_LN, BolshoyOleniOstrov, Bu_kk_Culture_Ungarn_MN, ChalmnyVarre, CHG, EHG, Esperstedt_MN, Ganj_Dareh_Iran_Neolithikum, Hungary_MN, Hungary_Neolithic, Iran_Chalcolithic, JK2065, Koros_Hungary_EN, Kunda, Latvia_HG3, Latvia_MN1, LBK_EN, LBK_Hungary_EN, Levanluhta, Narva, PWC_Sweden_NHG.SG, Scandinavia_LNBA, SHG, Sweden_HG.SG, TRB, Ukraine_HG1, Ukraine_N1, WHG, Yamnaya_Samara.
Wir stellen fest, dass K = 11 den geringsten Kreuzvalidierungsfehler ergibt, wie in der ergänzenden Abbildung 4b gezeigt.
Y-chromosomale Haplotypisierung
Wir ordneten antike Männer Y-Haplogruppen zu, indem wir das Programm yHaplo (https://github.com/23andMe/yhaplo)75 verwendeten. Kurz gesagt bietet dieses Programm eine automatisierte Suche durch den Y-Haplogruppen-Baum (wie in yHaplo bereitgestellt, abgerufen von ISOGG am 04.01.2016) von der Wurzel bis zum stromabwärts gelegenen Zweig, basierend auf dem Vorhandensein von abgeleiteten Allelen, und ordnet die am weitesten stromabwärts gelegene Haplogruppe mit abgeleiteten Allelen zu. Für etwa 15.000 Y-chromosomale SNPs, die sowohl in unserem Capture-Panel als auch in zwei veröffentlichten Datensätzen76,77 vorhanden sind, haben wir eine einzelne Base zufällig entnommen und als haploiden Genotyp verwendet. Wir verwendeten ein benutzerdefiniertes Skript, um EigenStrat-Genotypen in das yHaplo-Format zu konvertieren. Wir berichten über die am weitesten stromabwärts zugeordnete Haplogruppe, die vom Programm bereitgestellt wurde (Tabelle 1). Wir überprüften auch manuell den Ableitungsstatus und das Fehlen von Mutationen, die die zugewiesene Haplogruppe definieren, da fehlende Informationen zu einem vorzeitigen Stopp der automatischen Suche führen könnten.
Mitochondriale Haplotypisierung
Wir importierten die getrimmten mitochondrialen Reads für jedes Individuum mit Mapping-Qualität >30 in Geneious (Version 10.0.9, https://www.geneious.com)68 und reassemblierten diese Reads zum Referenzgenom RSRS78, indem wir den Geneious-Mapper mit mittlerer Sensitivität und 5 Iterationen verwendeten. Wir benutzten den eingebauten automatischen Varianten-Caller in Geneious, um mitochondriale Polymorphismen mit einer minimalen Abdeckung von 3 und einer minimalen Variantenhäufigkeit von 0,67 zu finden. Die resultierenden Varianten wurden nach Excel exportiert und manuell mit den SNPs verglichen, die in der online mtDNA-Phylogenie (mtDNA tree Build 17, 18 Feb 2016, http://www.phylotree.org/) berichtet wurden. Die Nukleotidpositionen 309.1 C(C), 315.1C, AC-Indels bei 515-522, 16182C, 16183C, 16193.1C(C) und 16519 wurden maskiert und nicht in unsere Haplotyp-Calls einbezogen.
Phänotypische SNPs
Wir verwendeten samtools mpileup (Version 1.3)65 , Filterung für Map- (-Q) und Base- (-q) Qualität von 30, Deaktivierung der per-Base Alignment Quality (-B), auf den getrimmten bam-Dateien, um ein Pileup von Reads zu generieren, die einem Satz von 43 phänotypischen SNPs4,40,41,79 zugeordnet sind, die Teil unseres Genome Capture Panels sind. Ein benutzerdefiniertes Python-Skript wurde verwendet, um das Pileup in eine Tabelle zu parsen, die die Anzahl der Reads enthält, die jedes Allel unterstützen (Supplementary Data 2).
Verfügbarkeit des Codes
Die gesamte Software, die in dieser Studie zuerst beschrieben wurde, ist in Online-Repositories frei verfügbar. Sex.DetERRmine: https://github.com/TCLamnidis/Sex.DetERRmine
ContaminateGenotypes: https://github.com/TCLamnidis/ContaminateGenotypes