A mosaic genetic structure of the human population living in the South Baltic region during the Iron Age

Despite the increase in our knowledge about the factors that shaped the genetic structure of the human population in Europe, the demographic processes that occurred during and after the Early Bronze Age (EBA) in Central-East Europe remain unclear. To fill the gap, we isolated and sequenced DNAs of 60 individuals from Kowalewko, a bi-ritual cemetery of the Iron Age (IA) Wielbark culture, located between the Oder and Vistula rivers (Kow-OVIA population). The collected data revealed high genetic diversity of Kow-OVIA, suggesting that it was not a small isolated population. Analyses of mtDNA haplogroup frequencies and genetic distances performed for Kow-OVIA and other ancient European populations showed that Kow-OVIA was most closely linked to the Jutland Iron Age (JIA) population. However, the relationship of both populations to the preceding Late Neolithic (LN) and EBA populations were different. We found that this phenomenon is most likely the consequence of the distinct genetic history observed for Kow-OVIA women and men. Females were related to the Early-Middle Neolithic farmers, whereas males were related to JIA and LN Bell Beakers. In general, our findings disclose the mechanisms that could underlie the formation of the local genetic substructures in the South Baltic region during the IA.


Results
aDNA extraction, sequencing and preliminary characteristics. The studied group included people whose remains were excavated from the IA Wielbark culture cemetery, located in Kowalewko ( Supplementary  Fig. S1). There were two major reasons for choosing this particular cemetery. First, it is located in the central part of the region that spreads between the Oder and Vistula rivers. The genetic structure of the population living in this region during the IA and before has been poorly characterized. Second, the cemetery was used by the local population for 200 years. Therefore, it may be assumed that the obtained genetic profile would characterize reasonably well the population that lived there for a considerable period of time. The study included 60 individuals, whose burial sites have been thoroughly characterized archaeologically and anthropologically (detailed  Supplementary Table S1). Because they represented the population living between the Oder and Vistula rivers during the IA, we called this group Kow-OVIA.
Radiocarbon dating (Supplementary Fig. S2 and Table S1) confirmed that the collected bone materials were from the first and second centuries (between AD 36 (+−30) and AD 181 (+−30)). aDNA was isolated from teeth, according to the procedure described by 24 and 25 . In addition, for three individuals, aDNA was isolated twice, each time from a different tooth (internal control). In total, 63 aDNA samples were obtained. The NGS libraries were successfully obtained for 60 samples (for 57 individuals, including two libraries for three individuals from whom two samples were extracted). To estimate the amount of human DNA, all of the obtained libraries were subjected to shallow sequencing (the average coverage of mtDNA genome was 3.62x and ranged between 1 and 11x). The average amount of human DNA in the examined samples was 12% and ranged from 0.1 to 91.9% (Supplementary Table S2). In all of the samples, human DNA showed a typical aDNA damage pattern (C > T and G > A alteration most frequently occurring at the 5′ and 3′ ends of the reads) (Supplementary Table S3). The data obtained from the shallow sequencing enabled mtDNA haplogroup definition for 23 individuals. Importantly, the haplogroups obtained for the samples isolated from the same individuals were identical. Furthermore, for 27 individuals, their genetic sex was identified. The obtained results are summarized in Supplementary Table S4.
The shallow sequencing allowed for determination of a full sequence of mtDNA genomes for 13 individuals. In addition, 23 samples (22 individuals) with high amounts of endogenous DNA (over 10%) were subjected to deep sequencing (the average coverage was 163x and ranged between 6x and 1092x). As a result, a full sequence of mtDNA genomes was determined for 22 additional individuals (Supplementary Table S4). In all cases, mtDNA haplogroups were also established. For 20 individuals, mtDNA haplogroups were determined twice, i.e., based on shallow and deep sequencing data. Importantly, in all cases, particular individuals were classified into the same haplogroups. Next we assessed the contamination levels in all samples. To this end, we used the software contam-Mix. Mean value for authenticity estimates was 96%. Two individuals (PCA0018 and PCA0063) had lower than average authenticity probability, therefore were excluded from subsequent analyzes (Supplementary Table 1). In total, the full sequence of mtDNA was determined for 33 individuals and the haplogroup was established for 40 individuals (Supplementary Table S4). All identified haplogroups had been found earlier in the populations living in Central Europe from the EBA to the IA and belonged to 9 haplogroups (H, N, I, J, K, T, U, W, and X) ( Table 1).
Intrapopulation genetic diversity of Kow-OVIA. The intrapopulation diversity of Kow-OVIA was analyzed using the full length mtDNA sequences obtained for 33 individuals (Supplementary Table S4). In a first step, we used a Minimum Spanning Network (MSN) method 26 to identify and visualize differences between particular mitochondrial genomes (Fig. 1). The obtained graph was typical for the population with a high level of genetic diversity 27 . We identified 3 pairs of individuals with identical mtDNA sequences (PCA0059, PCA0031; PCA0049, PCA0027; PCA0060, PCA0036) and 4 pairs of individuals with only one different nucleotide in the mtDNA sequences (PCA0047, PCA0003; PCA0004, PCA0028; PCA0034, PCA0059; PCA0034, PCA0031). To exclude the possibility that the presence of two potentially maternally related individuals in relatively small dataset will affect the results of our consecutive analyses, we removed from the tested group three individuals -one from each pair with identical mtDNA sequences. For the comparison, the results obtained with the dataset containing potentially related individuals were shown in the Supplementary Figures S4-S6.
Next, we applied Arlequin ver. 3.5.1 to determine the levels of haplotype diversity (HD) based on two fragments of the HVS (hypervariable sequence) region (HVS-I between 16033 and 16365 np (nucleotide positions) and HVS-II between 73 and 340 np) (Supplementary Table S5), and nucleotide diversity (π) based on the fragment of the HVS-I region (between 16000 and 16410 np) (Supplementary Table S6). The HD level calculated for Kow-OVIA (1.0000 +/−0.0082) was similar to those currently observed for the contemporary, open European populations and higher (even when considering 95% confidence interval) than the HD obtained for the population considered isolated (Supplementary Table S5) 28 . The level of π determined for Kow-OVIA (0.007922 +/−0.004666) was only slightly lower than the average value obtained for the present-day European populations (0.009 29 ) (Supplementary Table S6). Thus, the results of all three analyses, MSN, HD, and π, indicated that Kow-OVIA was not a genetically homogenous, isolated population.
Kow-OVIA relationships with other populations from European space-time. mtDNA haplogroup frequency. In the first step, we attempted to place Kow-OVIA in Central European space-time only. To this end, we compared mtDNA haplogroup frequencies in Kow-OVIA with other populations living in that region from the Mesolithic until the IA. The studied group, including Kow-OVIA and 13 other fossil populations (for details see Table 2 and Supplementary Table S7), was extended with the present-day Central Europe Metapopulation (CEM) (Supplementary Table S8). The obtained set of 15 populations called CEPT (Central European Population Transect) was subjected to unsupervised hierarchical clustering (Ward's method with Euclidean distance) (Fig. 2a). The obtained results generally complied with known chronology of the formation of the genetic structure of the Central European population. The first to be separated from other populations were the HG. Next, we observed a divergence of the Early Neolithic (EN  EPT) as earlier was subjected to hierarchical clustering and additionally to PCA ( Fig. 3a  , the most characteristic were both the mtDNA haplogroups typical for Asia (mainly C, D) and a high prevalence of haplogroups U5a and a lack of U5b (except the population of Bronze Age Kazakhstan (BAK)) for almost the entire analyzed period of time. Interestingly, the closely located in time and space Kow-OVIA and JIA occupied the opposite side of the PCA diagram ( Fig. 3b) than the preceding LN/EBA populations (the BEC, CWC, and UC). Such a placement resulted from the higher prevalence of the haplogroup H in Kow-OVIA and JIA compared to the BEC, CWC and UC. Alternatively, Kow-OVIA and JIA showed considerable differences in the prevalence of haplogroups typical for LN/EBA (I, R, T1, U2). They were less frequently observed in Kow-OVIA (2.7%) than in JIA (16.6%). As expected, the HG did not form a single group and differed from the other populations by a considerable portion of haplogroup U and a general lack of haplogroups associated with the "Neolithic package". Our analyses also revealed very interesting and unreported earlier changes in the frequency of haplogroups U5a/U5b in the Central European populations (Supplementary Table S8). Initially, during the LN/EBA, the haplogroup U5a, characteristic for the populations from eastern glacial refugia, dominated. Then, during the IA, haplogroup U5b, typical for the populations from the western glacial refugia, was more frequent. At present, haplogroup U5a is again more often found in the CEM.
We also placed Kow-OVIA on a map showing the matrilineal genetic structure of the present-day human population. To this end, we performed a PCA of haplogroup frequencies of Kow-OVIA and 73 extant worldwide populations (Supplementary Fig. S3a and Table S10, Supplementary Material Text). Kow-OVIA was located within a range of European genetic diversity, in close proximity to the present Scandinavian populations (Norwegian and Swedish).
Genetic distances. The relationships between Kow-OVIA and the earlier and later European populations were also determined by measuring genetic distances among them. In the first step, we focused on populations from Central Europe only (from the CEPT group). Our analysis involved a fragment of mtDNA HVS-I sequence (the region between 16064 and 16400 np) ( Fig. 2b and Supplementary Table S11 20 . **Anthropological assignment given in (), genetic sex assignment given after/. In the second step, we determined the genetic distances between Kow-OVIA and populations living in the whole of Europe from the Mesolithic age to the IA. To this end, we selected 26 populations from the EPT group for which sequences of the mtDNA HVS-I region between 16064 and 16400 np (Supplementary Table S12) were known for a considerable number of individuals. The calculated genetic distances (Supplementary Table S12) again showed that Kow-OVIA was closest to the BEC (Fst = −0.01097, p = 0.72712) and then to the JIA (Fst = −0.01054, p = 0.80627). The results were visualized on a multidimensional scaling (MDS) plot (Fig. 4). As in the analysis of the haplogroup frequencies, we observed that although separated by a short distance, Kow-OVIA and JIA were both differently located regarding particular LN/EBA populations (the BBC, CWC, and UC). The distance between the Kow-OVIA and BBC (Fst = 0.00342; p = 0.36707) is five times smaller than for the JIA -BBC pair (Fst = 0.02426; p = 0.07196). Simultaneously, the relationship is opposite when one compares the Kow-OVIA and JIA with the CWC and UC ( Fig. 4 and Supplementary Table S12). In addition, the performed analysis showed a previously unnoticed small genetic distance between the YAM and BBC (Fst = −0.00463, p = 0.58485), and the YAM and CWC (Fst = −0.00157, p = 0.50459).
Turning points in the formation of the Central European genetic structure. To verify the existing hypothesis describing the formation of the genetic structure of the Central European population in the period between the EN and IA, we used the analysis of molecular variance (AMOVA). In the first step, we determined an optimal division of the populations from CEPT excluding the HG and CEM. The division was considered to be optimal if intragroup variability (Fsc) was minimal and intergroup variability (Fct) maximal. In our analysis, we took advantage of earlier observations made by 2    results suggest that the JIA was closely related to the North-Central Europe populations, whereas Kow-OVIA was related to the populations of both the Eastern and Western regions of Europe. Next, we performed an analogous analysis but the studied group was extended by CEM (CEPT with excluded HG populations). Once again, we assumed that the populations of the EN/MN (STA, LBKT, LBK, RSC, SCG, BAC, and SMC) are a homogenous group, and we analyzed 33 combinations of the other fossil populations and CEM (Supplementary Table S14 Table S14). The results of the above analyses suggested that the contribution of the particular populations to the genetic structure of present-day Europe is not fully consistent with the current chronology of which populations inhabited Central Europe. The JIA, Kow-OVIA, BBC and CWC contributed more significantly to the genetic structure of present-day Europe than the BEC and UC.
Analysis of mtDNA lineages. To obtain more evidence to verify the above observations, we determined which populations gave rise to haplotypes found in Kow-OVIA (Supplementary Table S15). The same analysis of shared ancestral haplotypes 10 was performed for JIA (Supplementary Table S15). We observed a higher incidence of the HG mtDNA lineages in Kow-OVIA (9.68%) than in JIA (0%) and CEM (2%). Simultaneously, Kow-OVIA showed a lower incidence of the lineages originating from the LN/EBA (6.45%) compared to the JIA (29.16%) and CEM (13.62%).
Moreover, we performed a traditional analysis of shared haplotypes 30 (Supplementary Table S16), which demonstrated that approximately 29.4% of the haplotypes inherited by the Kow-OVIA were common for LN/ EBA populations (BBC, CWC, and UC). These haplotypes came from haplogroups H and U5a. The haplotype that belonged to U5a, according to current knowledge, first occurred in the LN/EBA, whereas the H haplotypes showed continuation from the EEF (Starčevo). The JIA inherited only approximately 23% of the haplotypes shared by LN/EBA populations (BBC, CWC, UC). Approximately 41% of Kow-OVIA and 19% of JIA haplotypes were not found in any LN/EBA population although the latter preceded them. Sex-linked factors affecting the genetic structure of Kow-OVIA. Finally, we attempted to determine whether the processes shaping the genetic structure of the Kow-OVIA population were biased by sex. To this end, based on anthropological and genetic sex assignment (Table 1), we divided the studied group into female and male subgroups -Kow-OVIA-F (17 individuals) and Kow-OVIA-M (14 individuals), respectively. Both subgroups were subjected to the analogical analyses as previously performed for the entire Kow-OVIA.
We found that the Kow-OVIA-F had a higher π value than the Kow-OVIA-M: 0.009185+/−0.005470 and 0.006987 0.004413, respectively (Supplementary Table S6). However, the difference was not statistically significant. Unsupervised hierarchical clustering involving CEPT showed that Kow-OVIA-M, similar to Kow-OVIA, grouped with the JIA and CEM whereas Kow-OVIA-F was placed in the diagram among the EN and MN populations ( Fig. 5a and Supplementary Table S8).
PCA analysis of Kow-OVIA-F, Kow-OVIA-M and EPT populations ( Fig. 5b and Supplementary Table S9) once again showed that Kow-OVIA-M was located similarly to Kow-OVIA, i.e., near the BBC and JIA due to considerable sharing of haplogroups H and U5b. However, the Kow-OVIA-F was grouped with the EN and MN populations, which was mainly a result of a high prevalence of the haplogroup K. The existence of genetic differences between females and males was also supported by the Fisher's test on the haplogroup frequencies (p = 0.013).  Table S18). The MDS analysis (Fig. 5c) demonstrated that the male individuals had smaller genetic distances to HG populations (HGCN, HGSW, HGE, and PWC) than did women and were closer to YAM than women (Fst = −0.00171, p = 0.48936; Fst = 0.01696, p = 0.21600, respectively). The analysis of shared haplotypes supported the described above tendencies -women from Kowalewko grouped with the EN and MN populations and were less connected with the BBC than men (Supplementary Table S19).

Discussion
Despite tremendous progress that has recently been made in archaeogenomic studies of European populations, our knowledge of the changes that occurred in the genetic structure of the peoples inhabiting Central Europe between the LN/EBA and Middle Ages is still very limited. One major problem is the scarcity of genetic data describing individuals who occupied the region east of the Oder river. To fill this gap, we performed an analysis of the mitochondrial genomes of a large group of people who lived in the area between the Oder and Vistula rivers (present-day western Poland) during the IA (the 1 st and 2 nd centuries AD). The studied group initially comprised 60 individuals; for 40 individuals, we assigned the mtDNA haplogroup, and for 33 of them, we determined the sequence of the entire mitochondrial genome. This dataset is particularly important to understand demographic processes that occurred in Central Europe during the IA and in the context of migrations between the 3 rd and 6 th centuries AD that are believed to have shaped the genetic landscape of contemporary Europe 31 .
According to common knowledge, the region between the Oder and Vistula rivers 2 tya was densely wooded with single isolated human settlements scattered in the forest. It is thought that inhabitants of these settlements had very limited contact not only with the so-called "civilized world" but also with neighboring populations. Interestingly, the picture emerging from the latest archaeological studies of burials and artifacts seems to be different. Recently, it has been suggested that the settlement in Greater Poland was compact and regular, with the appearance and disappearance of the sites occupied by the arriving groups interacting with the surrounding populations 32 and references therein. Our results strongly support the hypothesis presented by 32 . The analyses of intrapopulation variability revealed that the studied group was not an isolated population. The calculated HD and π levels fell within the range of values characteristic for contemporary open European populations.
Based on our results of mtDNA haplogroup frequency analysis, we were able to update a phylogenetic tree describing the history of H. sapiens in Central Europe (see Fig. 2a). Compared with the previous phylogenetic tree 2 , the new tree was expanded with the IA populations inhabiting the region between Oder and Vistula rivers and the region of Jutland. In the earlier constructed tree, CEM was closest to BBC despite the latter coexisting in Central Europe with the CWC (LN) and followed by the UC (EBA). In the new tree, Kow-OVIA was clustered with the LN/EBA populations. The present-day European population (CEM) was clustered with the JIA that was contemporary to Kow-OVIA and inhabited Jutland.
Based on the results presented above one can assume that there are some discrepancies between the hierarchical clustering (Fig. 3a) and PCA (Fig. 3b) of haplogroup frequencies in the EPT (for example, the placement of the Kow-OVIA in relation to the Treilles Culture population (TRE)). These differences indicate that in contrast to dendrogram, the PCA plot (capturing only first two principal components) did not show the full spectrum of EPT variability. We explored this issue by plotting additional principal components (Supplementary Figure S7a  and b). It demonstrated that their contribution to the explanation of the genetic variability of EPT populations was non-trivial.
Analysis of genetic distances (see Fig. 2b) showed that both JIA and Kow-OVIA, are the closest to the CEM. However, it should be mentioned that many of the resulting genetic proximities did not reach statistical significance at the alpha level 0.05 (mainly due to the multiple comparisons), thus they should be interpreted with  Table S11). Higher prevalence of the mtDNA haplogroup H in Kow-OVIA and JIA (its high level is also characteristic for the BBC) than in the preceding CWC and UC supports the hypothesis assuming significant demographic changes in Central Europe after the LN/EBA period 2 . This hypothesis is additionally strengthened by the results of AMOVA analysis indicating that there is some inconsistency between genetic distances and the chronology of the appearance of the studied populations in Central Europe, i.e., the older populations (BBC, CWC) contributed more to the genetic structure of CEM than the younger ones (UC).
Changes in the occurrence of mtDNA haplogroups U5a/U5b in Central Europe are also worth noting. At LN and EBA, the prevailing haplogroup was U5a for BBC/CWC/UC. Next, there was a dominance of U5b for the Kow-OVIA/JIA during IA and now U5a is again more popular (CEM). The first alteration in the U5a/U5b prevalence between the LN/EBA and the IA supports the hypothesis of demographic changes right after the LN, proposed by 2 . The second conversion indicated by our results suggests another crucial demographic event that should occur between the IA and present.
On the basis of the above observations, one may assume that in the IA, specific genetic substructures were formed in Central Europe. Because the demographic history of fossil populations often has a local character 33,34 , it is worth considering the range of the observed changes. These considerations should also take into account the hypothesis on the migrations that most likely occurred between the 3 rd and 6 th century AD. In this context, it seems necessary to compare Kow-OVIA and JIA with other populations from the IA, in particular those located east of Vistula, and with the populations that inhabited this region during the Middle Ages.
The process of demographic change in Central and Northern Europe after the LN/EBA appears to have a complex nature. Genetic proximity of the Kow-OVIA and JIA is consistent with the Kow-OVIA's affiliation to the Wielbark archaeological culture, strictly connected with Baltic regions 32 . However, despite close genetic distance between Kow-OVIA and the JIA, PCA placed them asymmetrically in relation to other ancient populations (see Fig. 3b). According to MDS (Fig. 4), JIA was close to the North-East European populations (CWC, BEC, UC) whereas the relationships of Kow-OVIA with other populations (earlier and contemporary to Kow-OVIA) were not so obvious. Interestingly, a small genetic distance between the JIA and UC was correlated with a high prevalence of the mtDNA haplogroup I in both populations. This result is consistent with earlier hypotheses suggesting that the genetic structure of a contemporary Danish population was formed not later than in the IA 35 . Unfortunately, our knowledge of haplogroup I prevalence in the Nordic Bronze Age is still scarce because of the small number of analyzed individuals; thus, it cannot be unambiguously stated that the observed proximity of the JIA and UC was, indeed, a result of the demographic changes after the LN/EBA. However, the above conclusion postulating a close connection between JIA and UC is also supported by the result of shared haplotype analysis. We found that 25% of the ancestral haplotypes found in the JIA were first reported in the UC and were not common in any earlier population. To better understand this phenomenon, the analysis of Y chromosome haplogroups is required, as the Nordic Bronze Age is characterized by the occurrence of the Y-haplogroups I and I1 18 , whereas the UC is mainly characterized by a higher prevalence of I2 13 . More detailed Y chromosome analyses involving a larger number of individuals would also shed more light on the process that resulted in a high prevalence of the mtDNA H haplogroup in the JIA, which is another signal of post-LN demographic changes 2 . In the case of Kow-OVIA, its genetic root in multiple European populations is evidenced by the fact that 55% of the ancestral haplotypes that were identified in this group were common to populations of the LN/EBA period (the CWC, BBC, and UC). For the JIA, the same origin is observed for only 19% of ancestral haplotypes.
Finally, we found that the genetic structures of female and male subpopulations of Kow-OVIA were significantly different. This fact cannot be explicitly determined based on the results of individual analyses; however, it is quite evident if one considers the whole set of data presented here including the Fisher test on haplogroup frequencies. The analyses of both mtDNA haplogroups and genetic distances indicated that women from Kowalewko were related closer to the EN/MN populations, and the men were closer to the CWC and UC. This observation may explain why the genetic relationships of Kow-OVIA with other ancient European populations were more complex and more difficult to define as it was in the case of JIA. In analyzing Kow-OVIA, we observed multiple overlapping effects of two subpopulations with different genetic affinities. One would speculate that the genetic profile of Kow-OVIA-F resulted from exogamy that was described for the CWC population 36 . This is, however, not the case. We found that the genetic differences between women and men were maintained for the entire observation period, i.e., for 200 years (approximately 8 generations). Such a composition of the genetic structure of Kow-OVIA could exist only if at least one subgroup (Kow-OVIA-F or -M) was periodically exchanged. It would further mean that Kowalewko played some specific roles in that region. According to the recent archaeological studies, the colonization pattern in IA Greater Poland could be linked with the existence of a centralized organization system 32 . Kowalewko could have been one of the important elements of this system. For example, it could have functioned as a garrison for the population closely associated with the JIA, such that warriors stayed in the garrison for only a few years and were then replaced by others. Other scenarios are also possible; however, verification of any hypothesis requires more detailed studies.

Materials and Methods
aDNA extraction and library preparation. For the purposes of this study, 63 samples (teeth) were obtained from 60 individuals from the Kowalewko archaeological site (Supplementary Table S1). After being transported to a clean aDNA lab, the teeth were cleaned with 5% NaOCl and rinsed with sterile water, followed by UV irradiation (254 nm) for 2 hours per site. The roots of the teeth were drilled using Dremel ® drill bits. Bone powder (approximately 250 mg) was digested with proteinase K and DNA-containing extract was purified using a silica-based method following 24 and 37 . Genomic libraries were prepared following a protocol from 38 (Günther et al. 2015). The PCR profile was as follows: initial denaturation (94 °C, 12 min), 12-16 cycles of 94 °C (30 s), 60 °C (30 s), 72 °C (45 s) and final extension (72 °C, 10 min). PCR reactions for the same library were pooled and purified with AMPure ® XP beads (Agencourt-Beckman Coulter) following 39 . Quality and size distribution of the libraries were verified with a High Sensitivity DNA kit and 2100 Bioanalyzer system (Agilent) while DNA concentration was determined with a Qubit fluorimeter and Qubit dsDNA HS Assay Kit (ThermoFisher Scientific), according to the manufacturers' protocols.
Next generation sequencing of aDNA libraries. For  Human DNA damage patterns. To examine data authenticity, we used mapDamage 2.0 46 to estimate human DNA damage parameters typical for aDNA for each sample: (i) λ, the fraction of nucleotides positioned in the single-stranded DNA overhang context; (ii) δs, C→T substitution rate in the single-stranded overhang context; and (iii) δd, C→T substitution rate in the double-stranded DNA context (Supplementary Table S3).

Contamination assessment.
To estimate the level of contamination with contemporary human DNA or with other human aDNA we used the software contamMix to compare the alignment rates between estimated sample's consensus mtDNA sequence and 311 potential contaminant mtDNA sequences.  Table S5). π was calculated for the fragment of mtDNA HVS-I region (16000-16410 np) (Supplementary Table S6). mtDNA haplogroup frequency analyses. Haplogroups were predicted with HaploFind 48 with respect to Phylotree build 17 (http://www.phylotree.org/) 49 . Only samples with a haplogroup score ≥ 0.8 and average coverage ≥ 3 were used in downstream analyses. To assess temporal mtDNA haplogroup frequency changes from Mesolithic to the present day, we performed Ward clustering with Euclidean distance on 23 haplogroups (H, H5, HV, HV0, V, I, J, K, N, N1a, R, T1, T2, U, U2, U3, U4, U5a, U5b, U8, W, X and others). We used a set of 40 samples from this study and 13 populations from Central/North Europe (Supplementary Table S7), and a generated Central European Metapopulation (CEM) composed of 500 random individuals sampled from the extant populations of Poland, Czech Republic, Germany and Austria, as in 2 (Supplementary Table S8).

Analysis of intrapopulation genetic diversity. To assess intrapopulation diversity of
To compare mtDNA variability of Kow-OVIA in a broader geographical context, we also applied unsupervised hierarchical clustering with the Ward method and Euclidean distance, and PCA on haplogroup frequencies of Kow-OVIA and published ancient mtDNA data from Mesolithic, Neolithic, Bronze Age and Iron Age populations from across Europe and West Asia (Supplementary Table S7 Table S10).
Cluster significance was tested by performing 10,000 permutations with the pvclust package in R ver. 3.3.0. PCA was conducted with prcomp of the vegan package in R ver. 3.3.0 (http://R-project.org).

Genetic distance analyses.
For sequence based analyses, the longest mtDNA HVS-I fragment present in the biggest fraction of published samples was selected (16064-16400 np). Additionally, for newly reported samples from this study, no missing nucleotides were allowed in the selected HVS-I range, and at least 3x coverage was expected for 95% of the nucleotides. To examine genetic affinities on the sequence level, we calculated genetic distances (Fst) 50 between two sample sets: CEPT (Supplementary Table S11 and S17) and EPT (Supplementary  Table S12 and S18), including only those individuals for which the 16064-16400 fragment of mtDNA HVS-I was present.
Pairwise and Slatkin's Fst 51 values were calculated in Arlequin ver. 3.5.1 for both datasets separately, with associated substitution model and gamma values selected with jModel test 0.1 AIC and BIC 52 . P values were calculated by performing 10,000 permutations. Genetic distances were visualized on an MDS plot with metaMDS function from the vegan package in R ver. 3.3.0.
Analysis of genetic structure. To examine if genetic affinities between particular populations from CEPT were a result of shared genetic structure, we conducted AMOVA analysis in 18 combinations of CEPT populations (excluding HG and CEM) (Supplementary Table S13). Next, we tested the genetic contributions of ancient populations to the extant mtDNA variability by analyzing AMOVA results from 33 combinations of CEPT populations (excluding HG) (Supplementary Table S14). Statistical significance was obtained by performing 10,000 permutation tests. mtDNA lineage analyses. To further test genetic affinities and account for temporal succession of archaeological cultures in Central Europe, we conducted an analysis of shared ancestral haplotypes as described in 10 in CEPT (Supplementary Table S15) and classical shared haplotype analysis 30 in LN/EBA and IA populations (Supplementary Table S16 and S19).
For details, see the Supplementary Materials.