The Yamnaya Steppe Herders from Russia

new

Admin
Administrator

Posts: 80,377

The Yamnaya Steppe Herders from Russia Feb 22, 2021 0:49:18 GMT

Quote

Post by Admin on Feb 22, 2021 0:49:18 GMT

DISCUSSION

After the last glacial period, at the end of the 10th and the beginning of the 9th millennium BC, vast areas in the Eastern Baltic, Finland, and northern parts of Russia were populated relatively quickly by HG groups (42–44). Flint originating from several places in the Eastern Baltic and the European part of Russia and the similarities in lithic and bone technologies and artifacts suggest the existence of extensive social networks in the forest belt of Eastern and Northern Europe after it was populated (43, 45, 46). This has led to the hypothesis that the descendants of Final Paleolithic HG from eastern Poland to the central areas of European Russia took part in the process (8, 47) and remained connected to their origin, creating a somewhat stretched social network (4). However, these connections ceased after a few centuries, as can be witnessed from the production of stone tools from mostly local materials, and new geographically smaller social units appeared in the middle of the 9th millennium BCE (4). aDNA studies of this area have included Mesolithic individuals that are from the 8th millennium BC or younger and reveal two genetic groups: WHG in the Eastern Baltic and EHG in northwestern Russia (14–17, 27, 34). However, no human genomes from the settlement period have been published so far, leaving the genetic ancestry/ancestries of the settlers up for discussion. The individual PES001 from around 10,700 cal BCE presented here with 5× coverage provides evidence for EHG ancestry in northwestern Russia close to the time it was populated. This, in turn, raises the question—hopefully answered by future studies—of the ancestry of the first people of the Eastern Baltic: Did they have EHG ancestry, as suggested by the shared social network of the two areas at the time of settlement, and an influx of WHG ancestry people later, not accompanied by a change in archaeological material, or did WHG ancestry people live in the Eastern Baltic from the beginning, representing a case of groups with different ancestry sharing a similar culture?

The formation of Fatyanovo Culture is one of the main factors that affected the population, culture, and lifestyle of the previously hunter/fisher-gatherer culture of the Eastern European forest belt. The Fatyanovo Culture people were the first farmers in the area, and the arrival of the culture has been associated with migration (21, 24). This is supported by our results as the Stone Age HG and the Bronze Age Fatyanovo individuals are genetically clearly distinguishable. The sample size of our HG is low, but the three individuals form a genetically homogenous group with previously reported EHG individuals (14, 16), and the newly reported PES001 is the highest-coverage whole-genome sequenced EHG individual so far, providing a valuable resource for future studies. What is more, the Fatyanovo Culture individuals (similarly to other CWC people) have not only mostly Steppe ancestry but also some EF ancestry that was not present in the area before and thus excludes the northward migration of Yamnaya Culture people with only Steppe ancestry as the source of Fatyanovo Culture population. The strongest connections for Fatyanovo Culture in archaeological material can be seen with the Middle Dnieper Culture (23, 48) spread in present-day Belarus and Ukraine (49, 50). The territory of what is now Ukraine is where the most eastern individuals with European EF ancestry and the most western Yamnaya Culture individuals are from based on published genomic data (13, 51) (Fig. 1 and data S1). Furthermore, archaeological finds show that LBK reached western Ukraine around 5300 BCE (52), and the Yamnaya complex (burial mounds) arrived in southeastern Europe around 3000 BCE and spread further as far as Romania, Bulgaria, Serbia, and Hungary (53). This is in accordance with our genetic results as the two populations that proved to be plausible mixture sources for Fatyanovo, with the other source being either of the two Russian Yamnaya groups (Kalmykia or Samara), were Globular Amphora that includes individuals from Ukraine and Poland, and Trypillia that is composed of individuals from Ukraine. These findings suggest present-day Ukraine as the possible origin of the migration leading to the formation of the Fatyanovo Culture and of the Corded Ware cultures in general.

The exact timing of and processes involved in the emergence of the Fatyanovo Culture in European Russia and the local processes following it have also remained unclear. Until recently, the Fatyanovo Culture was thought to have developed later than other CWC groups and over a longer period of time (21, 23). However, radiocarbon dates published last year (24), in combination with the 25 new dates presented here and the estimate for Yamnaya and EF admixture ~300 to 500 years before the Fatyanovo individuals of this study lived, point toward a fast process, similar in time to CWC people reaching the Eastern Baltic and southern Fennoscandia (54–56). The archaeological cultures are clearly differentiated between the areas. What is more, it has been suggested that the Fatyanovo Culture people admixed with the local Volosovo Culture HG after their arrival in European Russia (21, 57, 58). Our results do not support this as they do not reveal more HG ancestry in the Fatyanovo people compared to two other CWC groups; the three groups are shown to be similar by nonrejected one-way qpAdm models, and correlating radiocarbon dates with PC values or qpAdm ancestry proportions reveals no change in ancestry proportions of the Fatyanovo people during the period covered by our samples (2900 to 2050 BCE).

Last, allele frequency changes in western Russia and the Eastern Baltic revealed similar patterns in both areas: The frequencies of alleles in the MCM6 and FADS1-2 genes, which have been hypothesized to have changed due to dietary shifts from the Neolithic onward (16, 29, 59), change significantly during the Bronze Age, although the first signals of change can be seen already from the Neolithic. In accordance with a recent publication on lactase persistence (60), we find a low frequency of the rs4988235A allele in the initial Steppe ancestry samples [90% CI, 0 to 2.7% in (60), 0 to 33.8% in this study]. This suggests that the factors affecting these allele frequency shifts over time were complex and may have involved several environmental factors and genetic forces, as has already been suggested previously (14, 16, 18, 28, 29, 60, 61).

Last Edit: Feb 22, 2021 0:49:37 GMT by Admin

Admin
Administrator

Posts: 80,377

The Yamnaya Steppe Herders from Russia Feb 22, 2021 3:22:18 GMT

Quote

Post by Admin on Feb 22, 2021 3:22:18 GMT

MATERIALS AND METHODS
Experimental design

The teeth used for DNA extraction were obtained with relevant institutional permissions from the Institute of Ethnology and Anthropology of Russian Academy of Sciences (Russia), Cherepovets Museum Association (Russia), and Archaeological Research Collection of Tallinn University (Estonia). DNA was extracted from the teeth of 48 individuals: 3 from Stone Age HGs from western Russia (WeRuHG; 10,800 to 4250 cal BCE), 44 from Bronze Age Fatyanovo Culture individuals from western Russia (Fatyanovo; 2900 to 2050 cal BCE), and 1 from a Corded Ware Culture individual from Estonia (EstCWC; 2850 to 2500 cal BCE) (Fig. 1, data S1, table S1, and text S1). Petrous bones of 13 of the Fatyanovo Culture individuals have been sampled for another project. More detailed information about the archaeological periods and the specific sites and burials of this study is given below.

All of the laboratory work was performed in dedicated aDNA laboratories of the Institute of Genomics, University of Tartu. The library quantification and sequencing were performed at the Institute of Genomics Core Facility, University of Tartu. The main steps of the laboratory work are detailed below.
Laboratory methods

DNA extraction. The teeth of 48 individuals were used to extract DNA. One individual was sampled twice from different teeth. Apical tooth roots were cut off with a drill and used for extraction because root cementum has been shown to contain more endogenous DNA than crown dentine (62). The root pieces were used whole to avoid heat damage during powdering with a drill and to reduce the risk of cross-contamination between samples. Contaminants were removed from the surface of tooth roots by soaking in 6% bleach for 5 min, then rinsing three times with Milli-Q water (Millipore), and lastly soaking in 70% ethanol for 2 min, shaking the tubes during each round to dislodge particles. Last, the samples were left to dry under an ultraviolet light for 2 hours.

Next, the samples were weighed, [20 * sample mass (mg)] μl of EDTA and [sample mass (mg) / 2] μl of proteinase K were added, and the samples were left to digest for 72 hours on a rotating mixer at 20°C to compensate for the smaller surface area of the whole root compared to powder. Undigested material was stored for a second DNA extraction if need be.

The DNA solution was concentrated to 250 μl (Vivaspin Turbo 15, 30,000 MWCO PES, Sartorius) and purified in large-volume columns (High Pure Viral Nucleic Acid Large Volume Kit, Roche) using 2.5 ml of PB buffer, 1 ml of PE buffer, and 100 μl of EB buffer (MinElute PCR Purification Kit, QIAGEN).

Library preparation. Sequencing libraries were built using NEBNext DNA Library Prep Master Mix Set for 454 (E6070, New England Biolabs) and Illumina-specific adaptors (63) following established protocols (63–65). The end repair module was implemented using 30 μl of DNA extract, 12.5 μl of water, 5 μl of buffer, and 2.5 μl of enzyme mix, incubating at 20°C for 30 min. The samples were purified using 500 μl of PB and 650 μl of PE buffer and eluted in 30 μl of EB buffer (MinElute PCR Purification Kit, QIAGEN). The adaptor ligation module was implemented using 10 μl of buffer, 5 μl of T4 ligase, and 5 μl of adaptor mix (63), incubating at 20°C for 15 min. The samples were purified as in the previous step and eluted in 30 μl of EB buffer (MinElute PCR Purification Kit, QIAGEN). The adaptor fill-in module was implemented using 13 μl of water, 5 μl of buffer, and 2 μl of Bst DNA polymerase, incubating at 37°C for 30 min and at 80°C for 20 min. The libraries were amplified, and both the indexed and universal primers (NEBNext Multiplex Oligos for Illumina, New England Biolabs) were added by polymerase chain reaction (PCR) using HGS Diamond Taq DNA polymerase (Eurogentec). The samples were purified and eluted in 35 μl of EB buffer (MinElute PCR Purification Kit, QIAGEN). Three verification steps were implemented to make sure library preparation was successful and to measure the concentration of double-stranded DNA/sequencing libraries—fluorometric quantitation (Qubit, Thermo Fisher Scientific), parallel capillary electrophoresis (Fragment Analyzer, Agilent Technologies), and quantitative PCR. One sample (TIM004) had a DNA concentration lower than our threshold for sequencing and was hence excluded, leaving 48 samples from 47 individuals to be sequenced.

DNA sequencing. DNA was sequenced using the Illumina NextSeq 500 platform with the 75–base pair (bp) single-end method. First, 15 samples were sequenced together on one flow cell. Later, additional data were generated for some samples to increase coverage.
Statistical analysis

Mapping. Before mapping, the sequences of adaptors and indexes and poly-G tails occurring due to the specifics of the NextSeq 500 technology were cut from the ends of DNA sequences using cutadapt 1.11 (66). Sequences shorter than 30 bp were also removed with the same program to avoid random mapping of sequences from other species. The sequences were mapped to reference sequence GRCh37 (hs37d5) using Burrows-Wheeler Aligner (BWA 0.7.12) (67) and command mem with reseeding disabled.

After mapping, the sequences were converted to BAM format, and only sequences that mapped to the human genome were kept with samtools 1.3 (68). Next, data from different flow cell lanes were merged and duplicates were removed with picard 2.12 (http://broadinstitute.github.io/picard/index.html). Indels were realigned with GATK 3.5 (69), and lastly, reads with mapping quality under 10 were filtered out with samtools 1.3 (68).

The average endogenous DNA content (proportion of reads mapping to the human genome) for the 48 samples is 29% (table S1). The endogenous DNA content is variable as is common in aDNA studies, ranging from under 1 to around 78% (table S1).

aDNA authentication. As a result of degrading over time, aDNA can be distinguished from modern DNA by certain characteristics: short fragments and a high frequency of C→T substitutions at the 5′ ends of sequences due to cytosine deamination. The program mapDamage2.0 (70) was used to estimate the frequency of 5′ C→T transitions.

mtDNA contamination was estimated using the method from (71).This included calling an mtDNA consensus sequence based on reads with mapping quality of at least 30 and positions with at least 5× coverage, aligning the consensus with 311 other human mtDNA sequences from (71), mapping the original mtDNA reads to the consensus sequence, and running contamMix 1.0-10 with the reads mapping to the consensus and the 312 aligned mtDNA sequences while trimming seven bases from the ends of reads with the option trimBases. For the male individuals, contamination was also estimated on the basis of chrX using the two contamination estimation methods first described in (72) and incorporated in the ANGSD software (73) in the script contamination.R.

The samples show 10% C→T substitutions at the 5′ ends on average, ranging from 6 to 17% (table S1). The mtDNA contamination point estimate for samples with >5× mtDNA coverage ranges from 0.03 to 2.02% with an average of 0.4% (table S1). The average of the two chrX contamination methods of male individuals with average chrX coverage of >0.1× is between 0.4 and 0.87% with an average of 0.7% (table S1).

Kinship analysis. A total of 4,375,438 biallelic single-nucleotide variant sites, with minor allele frequency (MAF) > 0.1 in a set of more than 2000 high-coverage genomes of Estonian Genome Center (EGC) (74), were identified and called with ANGSD (73) command --doHaploCall from the 25 BAM files of 24 Fatyanovo individuals with coverage of >0.03×. The ANGSD output files were converted to .tped format as an input for the analyses with READ script to infer pairs with first- and second-degree relatedness (41).

The results are reported for the 100 most similar pairs of individuals of the 300 tested, and the analysis confirmed that the two samples from one individual (NIK008A and NIK008B) were indeed genetically identical (fig. S6). The data from the two samples from one individual were merged (NIK008AB) with samtools 1.3 option merge (68).

Calculating general statistics and determining genetic sex. Samtools 1.3 (68) option stats was used to determine the number of final reads, average read length, average coverage, etc. Genetic sex was calculated using the script sexing.py from (75), estimating the fraction of reads mapping to chrY out of all reads mapping to either X or Y chromosome.

The average coverage of the whole genome for the samples is between 0.00004× and 5.03× (table S1). Of these, 2 samples have an average coverage of >0.01×, 18 samples have >0.1×, 9 samples have >1×, 1 sample have around 5×, and the rest are lower than 0.01× (table S1). Genetic sexing confirms morphological sex estimates or provides additional information about the sex of the individuals involved in the study. Genetic sex was estimated for samples with an average genomic coverage of >0.005×. The study involves 16 females and 20 males (Table 1 and table S1).

Determining mtDNA hgs. The program bcftools (76) was used to produce VCF files for mitochondrial positions; genotype likelihoods were calculated using the option mpileup, and genotype calls were made using the option call. mtDNA hgs were determined by submitting the mtDNA VCF files to HaploGrep2 (77, 78). Subsequently, the results were checked by looking at all the identified polymorphisms and confirming the hg assignments in PhyloTree (78). Hgs for 41 of the 47 individuals were successfully determined (Table 1, fig. S1, and table S1).

No female samples have reads on the chrY consistent with a hg, indicating that levels of male contamination are negligible. Hgs for 17 (with coverage of >0.005×) of the 20 males were successfully determined (Table 1 and tables S1 and S2).

chrY variant calling and hg determination. In total, 113,217 haplogroup informative chrY variants from regions that uniquely map to chrY (36, 79–82) were called as haploid from the BAM files of the samples using the --doHaploCall function in ANGSD (73). Derived and ancestral allele and hg annotations for each of the called variants were added using BEDTools 2.19.0 intersect option (83). Hg assignments of each individual sample were made manually by determining the hg with the highest proportion of informative positions called in the derived state in the given sample. chrY haplogrouping was blindly performed on all samples regardless of their sex assignment.

Genome-wide variant calling. Genome-wide variants were called with the ANGSD software (73) command --doHaploCall, sampling a random base for the positions that are present in the 1240K dataset (https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers).

Preparing the datasets for autosomal analyses. The HO array dataset (https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers) was used as the modern DNA background. Individuals from the 1240K dataset (https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers) were used as the aDNA background.

The data of the comparison datasets and of the individuals of this study were converted to BED format using PLINK 1.90 (http://pngu.mgh.harvard.edu/purcell/plink/) (84), and the datasets were merged. Two datasets were prepared for analyses: one with HO and 1240K individuals and the individuals of this study, where 584,901 autosomal SNPs of the HO dataset were kept; the other with 1240K individuals and the individuals of this study, where 1,136,395 autosomal and 48,284 chrX SNPs of the 1240K dataset were kept.

Individuals with <10,000 SNPs overlapping with the HO autosomal dataset were removed from further autosomal analyses, leaving 30 individuals of this study to be used in autosomal analyses. These included 3 from WeRuHG, 26 from Fatyanovo, and 1 from EstCWC (table S1).

Principal components analysis. To prepare for PCA, a reduced comparison sample set composed of 813 modern individuals from 53 populations of Europe, Caucasus, and Near East and 737 ancient individuals from 107 populations was assembled (tables S3 and S4). The data were converted to EIGENSTRAT format using the program convertf from the EIGENSOFT 7.2.0 package (85). PCA was performed with the program smartpca from the same package, projecting ancient individuals onto the components constructed based on the modern genotypes using the option lsqproject and trying to account for the shrinkage problem introduced by projecting by using the option autoshrink.

Admixture analysis. For Admixture analysis (86), the same ancient sample set was used as for PCA, and the modern sample set was increased to 1861 individuals from 144 populations from all over the world (tables S3 and S4). The analysis was carried out using ADMIXTURE 1.3 (86) with the P option, projecting ancient individuals into the genetic structure calculated on the modern dataset due to missing data in the ancient samples. The HO dataset of modern individuals was pruned to decrease linkage disequilibrium using the option indep-pairwise with parameters 1000 250 0.4 in PLINK 1.90 (http://pngu.mgh.harvard.edu/purcell/plink/) (84). This resulted in a set of 269,966 SNPs. Admixture was run on this set using K = 3 to K = 18 in 100 replicates. This enabled us to assess convergence of the different models. K = 10 and K = 9 were the models with the largest number of inferred genetic clusters for which >10% of the runs that reached the highest log likelihood values yielded very similar results. This was used as a proxy to assume that the global likelihood maximum for this particular model was indeed reached. Then, the inferred genetic cluster proportions and allele frequencies of the best run at K = 9 were used to run Admixture to project the aDNA individuals, for which the intersection with the LD pruned modern dataset yielded data for more than 10,000 SNPs, on the inferred clusters. The same projecting approach was taken for all models for which there is good indication that the global likelihood maximum was reached (K3 to 18). We present all ancient individuals in fig. S2 but only population averages in Fig. 2B. The resulting membership proportions to K genetic clusters are sometimes called “ancestry components,” which can lead to overinterpretation of the results. The clustering itself is, however, an objective description of genetic structure and hence a valuable tool in population comparisons.

Outgroup f3 statistics. For calculating autosomal outgroup f3 statistics, the same ancient sample set as for previous analyses was used, and the modern sample set included 1177 individuals from 80 populations from Europe, Caucasus, Near East, Siberia and Central Asia, and Yoruba as outgroup (tables S3 and S4). The data were converted to EIGENSTRAT format using the program convertf from the EIGENSOFT 5.0.2 package (85). Outgroup f3 statistics of the form f3(Yoruba; West_Siberia_N/EHG/CentralRussiaHG/Fatyanovo/ Yamnaya_Samara/Poland_CWC/Baltic_CWC/Central_CWC, modern/ancient) were computed using the ADMIXTOOLS 6.0 program qp3Pop (87).

To allow chrX versus autosome comparison for ancient populations, outgroup f3 statistics using chrX SNPs were computed. To allow the use of the bigger number of positions in the 1240K over the HO dataset, Mbuti from the Simons Genome Diversity Project (88) was used as the outgroup. The outgroup f3 analyses of the form f3(Mbuti; West_Siberia_N/EHG/CentralRussiaHG/Fatyanovo/ Yamnaya_Samara/Poland_CWC/Baltic_CWC/Central_CWC, ancient) were run both using not only 1,136,395 autosomal SNPs but also 48,284 chrX positions available in the 1240K dataset. Because all children inherit half of their autosomal material from their father but only female children inherit their chrX from their father, then in this comparison chrX data give more information about the female and autosomal data about the male ancestors of a population.

The autosomal outgroup f3 results of the two different SNP sets were compared to each other and to the results based on the chrX positions of the 1240K dataset to see whether the SNPs used affect the trends seen. Outgroup f3 analyses were also run with the form f3(Mbuti; PES001/I0061/Sidelkino, Paleolithic/Mesolithic HG) and admixture f3 analyses with the form f3(Fatyanovo; Yamnaya, EF) using the autosomal positions of the 1240K dataset.

D statistics. D statistics of the form D(Yoruba, West_Siberia_N/EHG/CentralRussiaHG/Fatyanovo/ Yamnaya_Samara/Poland_CWC/Baltic_CWC/Central_CWC; Russian, modern/ancient) were calculated on the same dataset as outgroup f3 statistics (tables S3 and S4) using the autosomal positions of the HO dataset. The ADMIXTOOLS 6.0 package program qpDstat was used (87).

In addition, D statistics of the form D(Mbuti, ancient; Yamnaya_Samara, Fatyanovo/Baltic_CWC/ Central_CWC) and D(Mbuti, ancient; Poland_CWC/Baltic_CWC/ Central_CWC, Fatyanovo) were calculated using the autosomal positions of the 1240K dataset. However, comparing very similar populations directly using D statistics seems to be affected by batch biases—Central_CWC comes out as significantly closer to almost all populations than Fatyanovo, while this is not the case when comparing less similar Fatyanovo and Yamnaya_Samara. Because of this, the results of D(Mbuti, ancient; Poland_CWC/Baltic_CWC/Central_CWC, Fatyanovo) are not discussed in the main text, but the data are included in table S19.

FST. Weir and Cockerham pairwise average FST (89) was calculated for the dataset used for outgroup f3 and D statistics using the autosomal positions of the HO dataset using a custom script.

qpAdm. The ADMIXTOOLS 6.0 (87) package programs qpWave and qpAdm were used to estimate which populations and in which proportions are suitable proxies of admixture to form the populations or individuals of this study. The autosomal positions of the 1240K dataset were used. Only samples with more than 100,000 SNPs were used in the analyses. Mota, Ust-Ishim, Kostenki14, GoyetQ116, Vestonice16, MA1, AfontovaGora3, ElMiron, Villabruna, WHG, EHG, CHG, Iran_N, Natufian, Levant_N, and Anatolia_N (and Volosovo in some cases indicated in table S15) were used as right populations. Yamnaya_Samara or Yamnaya_Kalmykia was used as the left population representing Steppe ancestry. Levant_N, Anatolia_N, LBK_EN, Central_MN, Globular_Amphora, Trypillia, Ukraine_Eneolithic, or Ukraine_Neolithic was used as the left population representing EF ancestry. In some cases, WHG, EHG, WesternRussiaHG, or Volosovo was used as the left population representing HG ancestry. Alternatively, one-way models between Fatyanovo, Baltic_CWC, and Central_CWC were tested. Also, PES001 was modeled as a mixture of WHG and AfontovaGora3, MA1, or CHG.

To look at sex bias, four models that were not rejected using autosomal data were also tested using the 48,284 chrX positions of the 1240K dataset. The same samples were used as in the autosomal modeling.

ChromoPainter/NNLS. To infer the admixture proportions of ancient individuals, the ChromoPainter/NNLS pipeline was applied (28). Because of the low coverage of the ancient data, it is not possible to infer haplotypes, and the analysis was performed in unlinked mode (option -u). The autosomal positions of the HO dataset were used. Only samples with more than 20,000 SNPs were used in the analyses. Because ChromoPainter (90) does not tolerate missing data, every ancient target individual was iteratively painted together with one representative individual from potential source populations as recipients. All the remaining modern individuals from the sample set used for Admixture analysis were used as donors (tables S3 and S4). Subsequently, we reconstructed the profile of each target individual as a combination of two or more ancient individuals, using the non-negative least square approach. Let Xg and Yp be vectors summarizing the proportion of DNA that source and target individuals copy from each of the modern donor groups as inferred by ChromoPainter. Yp = β1X1 + β2X2 + … + βzXz was reconstructed using a slight modification of the nnls function in R (91) and implemented in GlobeTrotter (92) under the conditions βg ≥ 0 and ∑βg = 1. To evaluate the fitness of the NNLS estimation, we inferred the sum of the squared residual for every tested model (93). Models identified as plausible with qpAdm with Yamnaya_Samara and Globular_Amphora/Trypillia as sources were used. The resulting painting profiles, which summarize the fraction of the individual’s DNA inherited by each donor individual, were summed over individuals from the same population.

DATES. The time of admixture between Yamnaya and EF populations forming the Fatyanovo Culture population was estimated using the program DATES (37). The autosomal positions of the 1240K dataset were used.

Phenotyping. To predict eye, hair, and skin color in the ancient individuals (tables S20 to S22), the HIrisPlex-S variants from 19 genes in nine autosomes were selected (94–96), and the region to be analyzed was selected adding 2 Mb around each SNP, collapsing in the same region the variants separated by less than 5 Mb. A total of 10 regions (2 for chromosome 15 and 1 for each of the remaining autosomes) were obtained, ranging from about 6 to about 1.5 Mb. Similarly, to analyze the other phenotype-informative markers (diet, immunity, and diseases), 2 Mb around each variant was selected, and the overlapping regions were merged, for a total of 47 regions (45 regions in 17 autosomes and 2 regions on chrX). For the local imputation, we used a two-step pipeline (97) as follows: (i) variant calling, (ii) first imputation step using a reference panel as much similar as possible to the target samples, (iii) variant filtering, (iv) second imputation step using a larger worldwide reference panel, and (v) final variant filtering. This pipeline has been validated by randomly downsampling a high-coverage Neolithic sample (NE1) (98) to 0.05× and comparing the imputed variants in the low-coverage version with the called variants from the original genome. For a local imputation approach on 2 Mb, we obtained a concordance rate higher than 90% for all the variants, a figure that increased to 99% for frequent variants (MAF ≥ 0.3). The variants were called using ATLAS v0.9.0 (99) (task = call and method = MLE commands) (step 1) at biallelic SNPs with a MAF ≥ 0.1% in a reference panel composed of more than 2000 high-coverage Estonian genomes (EGC) (74). The variants were called separately for each sample and merged in one VCF file per chromosomal region. The merged VCFs were used as input for the first step of our two-step imputation pipeline [genotype likelihood update; -gl command on Beagle 4.1 (100)], using the EGC panel as reference (step 2). Then, the variants with a genotype probability (GP) less than 0.99 were discarded (step 3), and the missing genotype was imputed with the -gt command of Beagle 5.0 (101) using the large HRC as reference panel (102), with the exception of variants rs333 and rs2430561 [not present in the HRC (Haplotype Reference Consortium)], imputed using the 1000 Genomes as reference panel (step 4) (103). Last, a second GP filter was applied to keep variants with GP ≥ 0.85 (step 5). Then, the 113 phenotype-informative SNPs were extracted, recoded, and organized in tables, using VCFtools (104), PLINK 1.9 (http://pngu.mgh.harvard.edu/purcell/plink/) (84), and R (91) (tables S21 and S22). The HIrisPlex-S variants were uploaded on the HIrisPlex webtool (https://hirisplex.erasmusmc.nl/) to perform the pigmentation prediction, after tabulating them according to the manual of the tool. Out of 41 variants of the HIrisPlex-s system, two markers were not analyzed, namely, the rs312262906 indel and the rare (MAF = 0 in the HRC) rs201326893 SNP, because of the difficulties in the imputation of such variants.

The 28 samples analyzed here were compared with 34 ancient samples from surrounding geographical regions from literature, gathering them in seven groups according to their region and/or culture: (i) 3 Western Russian Stone Age HGs (present study); (ii) 5 Latvian Mesolithic HGs (34); (iii) 7 Estonian and Latvian Corded Ware Culture farmers [present study and (27, 34)]; (iv) 24 Fatyanovo Culture individuals (present study); (v) 10 Estonian Bronze Age individuals (28); (vi) 9 Estonian and Ingrian Iron Age individuals (28); (vii) 4 Estonian Middle Age individuals (28). For each variant, an analysis of variance (ANOVA) test was performed between the seven groups, applying Bonferroni’s correction by the number of tested variants to set the significance threshold (table S20). For the significant variant, a Tukey test was performed to identify the significant pairs of groups.

Admin
Administrator

Posts: 80,377

The Yamnaya Steppe Herders from Russia Mar 6, 2021 3:59:54 GMT

Quote

Post by Admin on Mar 6, 2021 3:59:54 GMT

Ancient DNA reveals that people of European ancestry have lost a gene linked to tuberculosis (TB) susceptibility over centuries.

TB is one of the world’s deadliest diseases and is caused by Mycobacterium tuberculosis bacteria. People with two copies of a genetic variant called P1104A are more likely to develop symptoms of TB after being infected with the bacteria.

To trace the frequency of P1104A over time, Gaspard Kerner at the Pasteur Institute in France and his team analysed modern human DNA from around the world and compared it to more than 1000 samples of ancient DNA from Europeans from the past 10,000 years.

They found that the variant first appeared in ancient DNA in low numbers around 8500 years ago in Western Eurasia. Using simulations and demographic models to date the origins and movements of this variant, the team predicted it may have originated in the same region around 30,000 years ago, long before the existence of TB in Europe. “It may have appeared randomly, like when animals have mutations in their genome,” says Kerner.

It then spread across central Europe 5000 years ago, and reached its highest frequency 3000 years ago, with around 10 per cent of the population carrying P1104A.

Kerner says it was able to spread without affecting an individual’s susceptibility to TB during that time as many people would only have one copy of the variant.

Read more: DNA reveals Romans helped spread TB across three continents: www.biorxiv.org/content/10.1101/210161v1

The frequency of the variant drastically decreased 2000 years ago, around the time modern TB bacteria became common. This may be because it was under strong negative selection from TB, Kerner says, as increasing migration made people more likely to inherit two copies of the variant and therefore become more susceptible to TB.

“Individuals carrying this mutation may have died faster than other individuals,” he says. The spread of TB during this time may have been aided by human migrations increasing populations and bringing new bacteria and diseases to Europe.

In modern Europeans and Americans, the variant appears in low frequencies, but it is absent in African and Eastern Asians populations. Kerner says this is consistent with the findings P1104A emerged in Eurasia, and that other genes may be behind the prevalence of TB in Africa and Asia today.

“People still get sick from TB, both in Europe and elsewhere,” says Vegard Eldholm at the Norwegian Institute of Public Health. Around 10 per cent of those infected with the bacteria develop TB. “This might reflect a long history of co-evolution, and humans having adapted to contain the infection. But it takes time for evolution to purge the gene,” Eldholm says.

Human ancient DNA analyses reveal the high burden of tuberculosis in Europeans over the last 2,000 years

Summary
Tuberculosis (TB), usually caused by Mycobacterium tuberculosis bacteria, is the first cause of death from an infectious disease at the worldwide scale, yet the mode and tempo of TB pressure on humans remain unknown. The recent discovery that homozygotes for the P1104A polymorphism of TYK2 are at higher risk to develop clinical forms of TB provided the first evidence of a common, monogenic predisposition to TB, offering a unique opportunity to inform on human co-evolution with a deadly pathogen. Here, we investigate the history of human exposure to TB by determining the evolutionary trajectory of the TYK2 P1104A variant in Europe, where TB is considered to be the deadliest documented infectious disease. Leveraging a large dataset of 1,013 ancient human genomes and using an approximate Bayesian computation approach, we find that the P1104A variant originated in the common ancestors ofWest Eurasians 30,000 years ago. Furthermore, we show that, following large-scale population movements of Anatolian Neolithic farmers and Eurasian steppe herders into Europe, P1104A has markedly fluctuated in frequency over the last 10,000 years of European history, with a dramatic decrease in frequency after the Bronze Age. Our analyses indicate that such a frequency drop is attributable to strong negative selection starting 2,000 years ago, with a relative fitness reduction on homozygotes of 20%, among the highest in the human genome. Together, our results provide genetic evidence that TB has imposed a heavy burden on European health over the last two millennia.

Journal reference: The American Journal of Human Genetics, 10.1016/j.ajhg.2021.02.009., DOI: 10.1016/j.ajhg.2021.02.009.

Last Edit: Mar 6, 2021 6:00:26 GMT by Admin

Admin
Administrator

Posts: 80,377

The Yamnaya Steppe Herders from Russia Mar 6, 2021 21:06:47 GMT

Quote

Post by Admin on Mar 6, 2021 21:06:47 GMT

Infectious diseases have been the leading cause of mortality
since the origin of modern humans in Africa and
throughout their subsequent dispersals around the
world.1–5 Tuberculosis (TB [MIM: 607948]) is considered
to be the deadliest infection of the common era, with
more than one billion deaths over the last 2,000 years,6–8
and still responsible for more than 1.5 million deaths
annually according to the WHO. The human genetic basis
of TB susceptibility has remained elusive until the turn of
the 21st century, when two rare inborn errors of immunity,
autosomal-recessive interleukin-12 receptor b1 (IL-12Rb1)
and tyrosine kinase 2 (TYK2) deficiencies, were identified
in children with severe TB.9,10 It was only in 2018 that
the first common, monogenic predisposition to TB was
identified. Homozygotes for the TYK2 (MIM: 611521)
P1104A polymorphism (rs34536443) were found to be at
higher risk of developing clinical forms of TB, due to the selective
disruption of IL-23-dependent antimycobacterial
IFN-g immunity, underlying a recessive trait.11 A subsequent
study revealed an enrichment in P1104A homozygotes
among TB cases of a case-control cohort from the
United Kingdom, where the allele is most prevalent today
(4%).7 The frequency of P1104A, together with its high
penetrance for TB in the homozygous state (>0.8),11 suggests
that about 1/600 British individuals would develop
TB during their lifetime because of the mutation, if TB
were still highly endemic in Europe.

Pathogen-imposed selective pressures have been paramount
during human evolution.2,4,5 Over the last decade,
population genetic studies have documented strong,
distinct selection signatures among host defense genes,
helping to delineate immunological mechanisms of major
importance,12 and supporting the notion that microbes
have had an overwhelming impact on human genome diversity.
4,5 While several studies have provided insight into
the periods when malaria has exerted pressure on humans,
13–17 little is known about the historical burden of
other infectious diseases associated with past epidemics.
Yet, TB appears to have been more lethal than malaria in
the common era,6 making it a stronger selective pressure
in endemic regions. Recent evidence based on mycobacterial
ancient DNA (aDNA) suggests a Holocene dispersal of
M. tuberculosis <6,000 years ago (ya),18,19 a time frame
that coincides with the growth of agricultural communities
and anthropogenic environmental changes, which
may have favored infectious disease transmission.20
To investigate the historical burden of TB in humans, we
sought to reconstruct the evolutionary history of the TYK2
P1104A variant. Indeed, this mutation, in the homozygous
state, underlies the only known common, monogenic predisposition
to TB.7,11 Moreover, TYK2 P1104A does not
affect the risk for other infectious diseases except, to a
milder degree, rare cases of infection by environmental
mycobacteria in otherwise healthy individuals.11 Whereas
disease-protective variants may rapidly increase in frequency
owing to positive Darwinian selection, diseaserisk
alleles are expected to evolve under strong negative
selection and be gradually purged from the population.
Because negatively selected variants have become rare,
very rare, or even extinct, they are harder to study using genetic
data from modern human populations. However,
with the increasing availability of genomes from ancient
individuals, direct measurements of the intensity of selection
are now possible, as significant increases or decreases
of allele frequencies can be captured with aDNA from
time transects.21 Thus, the study of the P1104A variant offers
an unprecedented opportunity to shed light on the
evolutionary history of a deadly human disease such as
TB. Of note, P1104A homozygotes have also been shown
to enjoy from a protective effect against various autoimmune
and inflammatory diseases.22,23 While this effect
could have provided a fitness advantage opposed to that
attributable to TB infection, the general late onset manifestation
of autoimmune and inflammatory disorders makes
unlikely the occurrence of a large counteractive effect.

We therefore examined the frequency trajectory of
P1104A over the last 10,000 years of European history,
by screening a collection of 1,013 genomes that cover a
time transect from the Mesolithic period to the Middle
Ages (Figure 1A; Table S1). We partitioned the aDNA data
into seven epochs and incorporated data from presentday
populations (supplemental material and methods).
The P1104A variant, which we found to be the result of a
single mutational event (Figure S1), appeared for the first
time in our dataset during the early Neolithic 8,500 ya
in the Anatolian peninsula, and then spread to Central Europe
where it remained at frequencies lower than 3% until
5,000 ya (Figures 1A–1C). During the Bronze Age,
P1104A increased in frequency, reaching its maximum frequency
3,000 ya at nearly 10%. After the Iron Age, we
observed a strong and consistent decrease in frequency of
P1104A, resulting in an average frequency of 2.9% among
contemporary Europeans.24

We estimated the age of the TYK2 P1104A mutation
(Tage), tested whether the mutation has been the substrate
of natural selection, and inferred the onset (Tonset) and
strength (s) of negative selection acting on homozygotes,
using an approximate Bayesian computation (ABC)
approach25 that considers large prior assumptions (Tage
U[8.5–100,000] ya, Tonset U[500–10,000] ya and s
U[0–1]; supplemental material and methods). We first
determined the extent to which our approach could determine
the evolutionary model of P1104A that best explains
the observed aDNA data, by comparing the fit of the simulated
to the observed data (supplemental material and
methods). We assumed a validated demographic model
for Europeans,26 to which we added gene flow from both
Near Easterners and Central Asians (Table S2), to account
for the large-scale migrations of early farmer populations
of the Anatolian plateau and Eurasian steppe populations
associated with the Yamnaya culture inferred from
aDNA.27 In doing so, considering the aforementioned
large prior assumptions, we obtained simulated frequency
trajectories that closely reproduce that of P1104A, similarly
to other genome-wide variants (Figure S2).We also noted a
similar, or higher, increase in frequency as that observed
for P1104A until the Bronze Age for more than 20% of
other aDNA variants within the uncertainty frequency interval
of P1104A in the Mesolithic ([0.00–0.10]; Table S3),
highlighting the marked impact of the aforementioned
migratory events on the frequency of a large fraction of
genomic variants, including P1104A. Furthermore, simulated
neutral variants closely matched observed frequency
distributions of non-coding variants for all epochs
(Figure S3), indicating that the demographic model
used—present-day Europeans are a mixture of Mesolithic
hunter-gatherers, Anatolian Neolithic farmers, and Eastern
steppe-related groups28,29—well reproduces the neutral
patterns of European diversity.

We then estimated the origin of the TYK2 P1104A mutation,
based on its frequency in K ¼ 12 populations sampled
at different epochs, including European aDNA data (Paleolithic,
Mesolithic, Early Neolithic, Late Neolithic, Bronze
Age, Iron Age, and Middle Ages; supplemental material
and methods) and present-day Europeans, Middle Easterners,
Central Asians (from 1% to 4%), Sub-Saharan Africans
(0%), and East Asians (0%) (Figure 2A; Table S3). We
found the age of P1104A to be 30,000 years old (mode
¼ 29,182; 95% CI [20,636–57,285]) (Figure 2B; supplemental
material and methods), which is consistent with
a previous estimate.30 Using cross-validation, we found
that parameter estimation was accurate across all ages,
with 96% of 1,000 estimated 95% CIs including the true
simulated value, and also robust to the choice of the summary
statistics used (Figures S4A–S4D). While the 95% CI
for the age of P1104A overlaps with the divergence time
between West and East Eurasians (35–45 kya), the proportion
of best-fitting simulated variants originating in the
common ancestors of West Eurasians was significantly
higher than that of the rest of simulated variants (OR ¼
7.00, 95% CI [5.70–8.53], p < 10
10; Figure 2B; supplemental
material and methods). This suggests that P1104A
originated in the common ancestors of West Eurasians after
the split with East Eurasians, but before the divergence
of Europeans, Middle Easterners, and Central Asians.
Together, our results provide robust evidence that TYK2
P1104A appeared during the Upper Paleolithic inWest Eurasia,
largely predating the estimated emergence of TB in
Europe.18,19,3

We then estimated the origin of the TYK2 P1104A mutation,
based on its frequency in K ¼ 12 populations sampled
at different epochs, including European aDNA data (Paleolithic,
Mesolithic, Early Neolithic, Late Neolithic, Bronze
Age, Iron Age, and Middle Ages; supplemental material
and methods) and present-day Europeans, Middle Easterners,
Central Asians (from 1% to 4%), Sub-Saharan Africans
(0%), and East Asians (0%) (Figure 2A; Table S3). We
found the age of P1104A to be 30,000 years old (mode
¼ 29,182; 95% CI [20,636–57,285]) (Figure 2B; supplemental
material and methods), which is consistent with
a previous estimate.30 Using cross-validation, we found
that parameter estimation was accurate across all ages,
with 96% of 1,000 estimated 95% CIs including the true
simulated value, and also robust to the choice of the summary
statistics used (Figures S4A–S4D). While the 95% CI
for the age of P1104A overlaps with the divergence time
between West and East Eurasians (35–45 kya), the proportion
of best-fitting simulated variants originating in the
common ancestors of West Eurasians was significantly
higher than that of the rest of simulated variants (OR ¼
7.00, 95% CI [5.70–8.53], p < 1010; Figure 2B; supplemental
material and methods). This suggests that P1104A
originated in the common ancestors of West Eurasians after
the split with East Eurasians, but before the divergence
of Europeans, Middle Easterners, and Central Asians.
Together, our results provide robust evidence that TYK2
P1104A appeared during the Upper Paleolithic inWest Eurasia,
largely predating the estimated emergence of TB in
Europe.18,19,3 the trajectories of 25% of the best fitting simulated
deleterious variants (s U[0–1] and Tonset U[500–10,000]; supplemental
material and methods), relative to only 1% of
the best fitting simulated neutral variants (OR ¼ 33, 95%
CI ¼ [5–240], p < 1010; Figure S5B; Table S4). These observations
collectively support a history of negative selection
driving the evolution of the TB-risk P1104A variant after
the Bronze Age.

Last Edit: Mar 6, 2021 21:08:45 GMT by Admin

Admin
Administrator

Posts: 80,377

The Yamnaya Steppe Herders from Russia Mar 7, 2021 5:14:10 GMT

Quote

Post by Admin on Mar 7, 2021 5:14:10 GMT

To quantify the degree of deleteriousness of TYK2 P1104A
during Europeanhistory,we verified that allele frequency trajectories
were informative to assess negative selection, and,
encouragingly, we observed a strong positive correlation between
drops in allele frequencies and s values (Figure S6A).
We first hypothesized that negative selection started with
the arrival of agriculture in Europe,20 a period that includes
the upper bound estimation for the most recent common
ancestor of the M. tuberculosis complex 6,000 ya.18,19 However,
such an early onset of selection (Tonset ¼ 10,000) was
clearly rejected by our simulations (Hotelling’s T-squared
test p ¼ 5.4 3 104; Figure S6B; supplemental material and
methods; Table S4), as no simulated variants were able to
reproduce the frequency increase of P1104A until the Bronze
Age.Conversely,whenallowing theonset of selectionto vary
across the last 10,000 years, using the former large priors (Tonset
U[500–10,000] ya and s U[0–1]), our best simulations
did not significantly differ fromP1104A (i.e., the simulation
set was not rejected; Hotelling’s T-squared test p ¼ 0.09) and
revealed that scenarios with recent onsets of negative selection
were those best fitting the data (Figure S6B).

To explain the strongest frequency increase and decrease
for P1104A, we modeled allele frequencies of K ¼ 5 ancient
populations (Late Neolithic, Bronze Age, Iron Age, and
Middle Ages) and present-day Europeans, and assumed
large priors for model parameters (supplemental material
and methods, Table S3). We found that negative selection
on P1104A homozygous carriers started 1,937 ya (95% CI
[500–7,912]), with a selection coefficient of 0.21 (95% CI
[0.06–0.82]) (Figures 3A–3C). This onset of selection is
consistent with a neutral evolution for the allele until
the Bronze Age, suggesting that drift and admixture are
sufficient to explain the increase of P1104A frequency until
this epoch. These estimations should not be biased
owing to read mapping bias of the reference allele in the
ancient genome dataset,32 given that 1104A is the alternative
allele (supplemental material and methods). Furthermore,
parameter estimation was found to be robust to
the choice of the summary statistics used, with the 95%
CIs of the estimates including the true simulated value
93% of the time (Figures S6C and S6D). Although our analysis
showed that the more recent the onset of selection was
the closer the frequency trajectory estimation was to the
empirical data (Figure S6A), the fit was found to be similar
within the last 2,000 years (Figure 3B), consistent with
our estimation. With respect to the selection coefficient,
the posterior distributions of s were shifted to 1 as Tonset
became closer to 0, and the general posterior distribution
for the strength of negative selection was similar to that
of onsets of selection occurring between 1,000 and 3,000
ya (Figures 3A and 3C). Importantly, consistent ABC estimates
of the strength and the onset of selection were found
when either excluding the Iron Age, i.e., the epoch with
smallest sample size (s ¼ 0.19; 95% ¼ [0.03–0.83]; Tonset
¼ 1,670 ya; 95% CI ¼ [500–8,388] ya) or when using the
whole European frequency trajectory, i.e., from the Paleolithic
to the present (s ¼ 0.21; 95% ¼ [0.04–0.84]; Tonset ¼
1,567 ya; 95% CI ¼ [500–8,367]).
Using the same approach, we estimated the selection coefficient
of another mutation, TYK2 I684S, a missense
variant that is neither in linkage disequilibrium with
P1104A nor associated with TB risk,11 and found values
that were compatible with neutrality (s ¼ 0.02; 95% CI
[0–0.19]; Figures S7A and S7B). Thus, our analyses support
the notion that, despite the reported protective effects of
P1104A against some immune-related disorders,22,23 TB
has exerted pressure on the TYK2 P1104A variant over
the last 2,000 years, with a 20% relative fitness reduction
for homozygotes at each generation since.
Finally, we sought to apply the same approach to reported
pathogenic variants, by cross-matching the ClinVar database
33 with aDNA variants present in our cohort that fall
into the uncertainty range of P1104A in the Bronze Age
([0.04–0.10], Figure 1B). Among the resulting three variants
with a ‘‘pathogenic’’ clinical significance annotation, only
one (HFE C282Y [MIM: 613609]) presents a frequency
decrease across the last four epochs. HFE C282Y is a known
disease-causing variant underlying hemochromatosis, an
autosomal-recessive autoimmune disease (HFE1 [MIM:
235200]) that impairsmineral metabolism,which can affect
the growth and clearance of intra- and extra-cellular pathogens.
34 HFE C282Y reached its maximum frequency, of
nearly 10%,during theMiddleAges andthen decreased to its
present-day frequency of 4%. Consistent with our expectations,
we found a similarly strong selection coefficient of
0.20 (mode ¼ 0.22; 95% CI [0.03–0.76]; Figure S7B), and an
onset of negative selection during the Middle Ages (mode
¼ 724 ya; 95% CI [500–7,508]).