We assembled a comprehensive ancient DNA dataset from the Allen Ancient DNA Resource (AADR) maintained by the David Reich Laboratory at Harvard Medical School (version 54.1, released January 2024). The initial collection comprised 3,600 ancient DNA samples from Eurasia spanning 800 BCE to 1,200 CE. The geographic distribution of our 1,800 samples is illustrated in Figure 2, extending from the Baltic region to the Balkans and covering the period from 500 BCE to 1200 CE.

Our regional representation imcluded 400 samples from the Hungarian Plain and Carpathian Basin dating from 300 BCE–900 CE, 185 samples from the Balkan Peninsula spanning 700 BCE–1,200 CE, and 133 samples from Central and Eastern Europe covering 300 BCE-900 CE. Additionally, we included 131 samples from Baltic and Scandinavian regions dating to 300 BCE-1,100 CE, and 44 samples from Western Steppe and Pontic regions from 400 BCE–800 CE. Sample selection prioritized individuals from archaeologically well-contextualized sites with clear cultural attribution and radiocarbon dating, with complete supporting metadata including site locations, archaeological context, and dating information. We provide full list of samples we used in the supplement.

Quality control procedures adhered to established protocols for ancient DNA analysis (Patterson et al., 2012; Maier and Reich, 2023).

We implemented stringent filtering criteria at both SNP and individual levels. For SNP-level filtering, we retained only those variants with minor allele frequency of at least 0.01, Hardy-Weinberg equilibrium p-values exceeding 10^–6, and maximum missing genotype rates not exceeding 80%. We further removed SNPs in linkage disequilibrium using an r ² threshold of 0.25, with 200 SNP windows and 50 SNP steps.

At the individual level, we required minimum coverage of 50,000 autosomal SNPs and maximum missing genotype rates of 80% and we verified consistency between genetic and archaeological sex determination. To avoid pseudoreplication, we performed kinship filtering to remove first-degree relatives identified by PI_HAT values exceeding 0.25. Following these quality control procedures, our final analytical dataset comprised 1,800 individuals genotyped at 540,000 autosomal SNPs.

Principal Component Analysis has been extensively applied in population genetics since Cavalli-Sforza (1994), with traditional approaches emphasizing the first two or three principal components that capture the greatest genetic variance.

We performed principal component analysis using PLINK 2.0 (alpha 3.7). We first computed standard PC1-2 to confirm the documented Germanic-Slavic discrimination problem, as shown in Supplementary Figure S1. We then implemented our systematic screening protocol through a four-stage process. First, we generated 20 principal components from the filtered genotype matrix for comprehensive PC calculation. Second, we systematically evaluated each PC to identify populations showing extreme values exceeding two standard deviations from the mean. Third, we assessed each PC’s capacity to separate early Slavic from Baltic and Germanic populations through discrimination testing. Finally, we tested all possible three-way PC combinations to identify the optimal combination for maximum population discrimination. The complete screening protocol is available at github.com/StamovS/slavic-signatures-avar/scripts/pc_screening_analysis.R.

Formal admixture tests and ancestry proportion estimation employed qpAdm and qpWave implemented in AdmixTools 2.0.1 (Harney et al., 2021). Our reference population framework consisted of primary references including CEE_EarlyMedieval representing Early Western Slavic populations and Latvia_BA representing Baltic Bronze Age. Secondary references encompassed Hungary_Avar populations totaling 256 individuals, Balkans_IA with 89 individuals, and Byzantine populations comprising 112 individuals. The outgroup panel included Mbuti, Han, Karitiana, Papuan, Onge, Iran_N, WHG, and MA1 to ensure robust model fitting.

Our model selection criteria involved testing between two and five source population models, requiring p-values exceeding 0.05 for acceptable model fit. We used qpWave to determine the minimum number of sources needed for adequate modeling and validated results using block bootstrap with 1,000 replicates. All qpAdm configuration files and parameters are available at github.com/StamovS/slavic-signatures-avar/qpadm_configs/.

For Euclidean distance calculations in PC space, we computed genetic distances using the formula d (i,j) = √(Σ_k=1 ³ (PC_ki - PC_kj)²), where PC_ki represents the kth principal component score for individual i, using optimal PC combinations, typically PC9, PC4, and PC3. For cluster analysis, we employed Ward’s method on Euclidean distance matrices for hierarchical clustering, while multidimensional scaling used classical metric MDS implemented in R.

All analyses were implemented in R version 4.3.1, utilizing PLINK 2.0 for genotype processing and PCA, AdmixTools 2.0.1 for formal admixture modeling, and tidyverse 2.0.0 with ggplot2 3.4.2 for data manipulation and visualization. We developed custom R scripts for systematic PC screening and population analysis. The complete analytical pipeline comprising five integrated modules is available at github. com/StamovS/slavic-signatures-avar with comprehensive documentation, error handling, and example datasets enabling full reproducibility.

Standard principal component analysis using PC1-2 (Supplementary Figure S1) confirmed the well-documented challenge of distinguishing Germanic from Slavic populations due to shared Baltic Bronze Age ancestry. These primary components, which capture the greatest genetic variance, place early Germanic and early Slavic individuals in overlapping clusters without clear separation, forcing previous studies to employ unexpected reference populations when modeling Balkan Slavic ancestry (Olalde et al., 2023).

Our systematic screening of all twenty principal components identified PC9 as capturing a distinctive ancestry signal that effectively discriminates early Slavic populations within European ancient DNA samples (Supplementary Table S1). While PC9 reaches maximum values globally in northern Eurasian populations including Nganasan (0.0567), Ket, and European Saami, within European samples this component consistently maximizes in early Western Slavic populations from Slovakia and Slovenia dating to the 6th-9th centuries CE.

Analysis of PC9 extremes (Supplementary Table S1) revealed CEE_EarlyMedieval samples show the highest European values. The systematic screening approach identified PC9 combined with PC3 and PC4 as optimal for Slavic ancestry detection (Figures 3A,B).

Analysis of PC9 extremes (Supplementary Table S1) revealed that CEE_EarlyMedieval samples show the highest European values (0.0186), followed by various Croatian medieval samples (Croatia_Brekinjova: 0.0162, Croatia_Jagodnjak: 0.0158, Croatia_Sibeník: 0.0153). The presence of Latvia_BA samples with elevated PC9 values (0.0143–0.0128) confirms the Baltic Bronze Age connection to early Slavic populations. The systematic screening approach identified PC9 combined with PC3 and PC4 as optimal for Slavic ancestry detection (Figures 3A,B).

The heat map analysis (Supplementary Figure S2) of genetic distances in PC9,4,3 space demonstrates clear population clustering that is absent in conventional PC1-2 analysis. Hungarian Avar-period Slavic samples (Hungary_EarlyAvarslav, Hungary_LateAvarslav) cluster tightly with Croatian medieval samples, while maintaining distinct separation from Germanic and other non-Slavic populations.

Our methodology successfully identified samples previously classified as Slavic by independent research teams, validating the systematic PC screening approach. These validation samples included Byzantine_oEuropean outlier samples) identified by Reich and Lazaridis as “Slavic admixed,” Czech_Medieval samples, displaying expected temporal and geographic patterns, Hungary_Avar context samples AV1 and AV5 independently identified as early Slavic individuals, and Croatia_Medieval_o samples considered by the Reich team as representing early Slavic presence.

The discriminating power of PC9 appears to capture an ancient Paleolithic Siberian ancestry component that was absorbed differentially by various Bronze Age populations. This component is absent in Mesolithic European hunter-gatherers but present in Bronze Age populations from the Sintashta-Alakul complex (2,100–1,800 BCE) and Kazakhstan MLBA cultures (1800–1500 BCE), suggesting its introduction through steppe-mediated gene flow.

Analysis of thirty individuals from the Steppe Barbaricum, the 150-km flatland strip between the Danube and Tisza rivers, revealed extraordinary genetic complexity that predates traditional accounts of Slavic presence by several centuries.

The population composition shows remarkable balance among five ancestral components (Supplementary Figure S3): Germanic (21.5%), Balkan Iron Age (20.8%), Proto-Slavic (20.2%), Byzantine (19.1%), and Steppe Sarmatian-Hunnic (16.9%). This balanced distribution, confirmed through qpAdm modeling with high statistical confidence (p = 0.73), suggests extensive intermarriage and cultural exchange rather than segregated ethnic communities.

Significantly, Proto-Slavic signatures absent in second-century samples emerge clearly in 5th-century individuals, suggesting Slavic presence in the Hungarian Plain by at least 400 CE—over a century before historical documentation. Temporal analysis reveals progressive increase in Proto-Slavic components from 0% (2nd century) to 20.2% (5th century), while Sarmatian components decrease from 35% to 16.9% over the same period.

This admixed Barbaricum population contributed substantially to later groups: 20%–40% to Early Western Slavic samples and approximately 20% to medieval Croatian populations. The linguistic evidence, including the possibly Slavic term “strava” recorded at Attila’s court (453 CE) and the ethnonym “Limigantes” (meaning “speakers”), gains new significance in light of these genetic findings.

Comprehensive analysis of nearly 300 individuals from Avar-period cemeteries (Supplementary Figure S5) in the Carpathian Basin revealed patterns challenging traditional narratives of strict ethnic hierarchies within the Khaganate.

Twelve individuals with predominantly Slavic genetic profiles (>50% Early Western Slavic ancestry) were identified in prestigious burial contexts, including graves with weapons, horse sacrifices, and high-status ornaments. Notable examples include SZOD1-829 from the 6th-7th century with 67% Slavic ancestry, buried with sword and belt set; OBT-56 from the 8th century with 72% Slavic ancestry in an elite warrior burial with horse; and Sample Av5 from 560 CE, representing the earliest confirmed Slavic individual in an Avar context.

The multidimensional scaling plot (Supplementary Figure S5) shows clear separation between different Avar-period populations. Hungary_EarlyAvarBalt clusters near Latvia_BA, confirming the presence of individuals with predominantly Baltic ancestry (59%–75%, Supplementary Figure S7) within the Khaganate. These likely represent either direct Baltic migrants or preserved unadmixed early Slavic populations.

Avar-Slavic genetic profiles (Supplementary Figure S8) from the late period (700–800 CE) show distinct composition patterns, incuding substantial Byzantine components (15%–20%), indicating ongoing interaction with imperial territories. The composition—Early Western Slavic (59.2%), Roman-Byzantine (20.9%), and Steppe Post-Hunnic (19.9%)—demonstrates that Slavic populations within the Khaganate maintained distinct genetic identity while incorporating regional elements.

Euclidean distance analysis in PC9, 4, 3 space (Supplementary Figure S6) reveals clear genetic relationships between Baltic and Central/Eastern European populations. Plotting the distance to the nearest Latvia_BA sample against the distance to the nearest CEE_EarlyMedieval sample showed that Croatian individuals form a continuous genetic gradient linking Baltic references to the Balkan populations from antiquity. Early Slavic individuals cluster tightly near the origin, reflecting their dual genetic affinity to both Baltic Bronze Age and Early Western Slavic populations. In contrast, present-day Slavic groups appear more dispersed, indicating additional post-migration admixture. Croatian samples (green) form a continuous gradient between Baltic reference and other European populations, supporting the multi-wave migration model.

Formal qpAdm modeling revealed a clear gradient of Baltic Bronze Age ancestry (Figure 4) that traces Slavic migration routes from the Baltic homeland to the Balkans.

Modern populations show decreasing Baltic ancestry with geographic distance from the Baltic region. Lithuanian. HO retains 71% Baltic Bronze Age ancestry, while Polish. DG shows 54.2%. Early Western Slavic populations from Slovakia and Slovenia maintain 57%, and the Hungary_AvarBalt individual (Supplementary Figure S7) shows 59%–75%. Croatian populations demonstrate intermediate levels, with Croatia_Medieval_o showing 39%–51% and Hungary_AvarSlav showing 24.6%.

Critically, non-Slavic Balkan populations show no significant Baltic component, with Albanian. DG, Greek. DG, and Balkans_IA all showing 0% Baltic ancestry. This gradient provides quantitative support for the Baltic origins hypothesis and demonstrates that Croatian populations retain 39%–51% of the original Baltic genetic signature, consistent with substantial but incomplete population replacement during medieval formation.

Statistical modeling demonstrates contemporary Croatian populations formed through at least three distinct waves (Figures 5A,B). Analysis of Croatia_Medieval_o samples, representing the earliest medieval Croatian individuals (7th-ninth century), reveals complex admixture patterns (Figure 5A). The optimal four-way model identifies Latvia_BA (Baltic) contributing 50.8%, Kazakhstan Sarmatian contributing 17.3%, Balkans_IA contributing 19.1%, and Steppe components contributing 12.8%. These high Baltic ancestry levels in early samples contrast with later medieval populations, supporting a multi-wave model with progressive admixture. The temporal flow diagram illustrates component transformations (Figure 5B).

Five different admixture models for contemporary Croatian populations converge on consistent patterns (Figure 6). Model 1 (AvarSlav focus) shows 22.9% Baltic, 37.1% Hungary_AvarSlav, and 40% Balkans/Roman ancestry. Model 2 (Medieval focus) reveals 46% Croatia_Medieval_o, 29% Roman, 20% Proto-Slavic, and 5% Baltic components. Model 3 (Steppe focus) identifies 51.6% Steppe Barbaricum, 27% Early Slavic, and 21.4% Balkans ancestry. Model 5 (Full model, best fit) demonstrates 53.7% combined Avar-period ancestry and 46.3% other components.

The consistency across models, showing 20%–25% pre-Slavic Balkan ancestry and 50%–60% various Slavic components, supports multiple migration waves with varying degrees of admixture. The temporal flow diagram illustrates how these components transformed across periods.

The evidence supports at least three distinct migration phases.

1. Early wave (6th century): High Baltic ancestry (>50%), minimal local admixture, documented by Gepid sources mentioning Slavs crossing at Iron Gates

Our systematic analysis of 1,800 ancient genomes provides unprecedented resolution of the complex population dynamics underlying Croatian ethnogenesis. By developing a methodological framework that overcomes the Germanic-Slavic discrimination problem through systematic principal component screening, we demonstrate that South Slavic genetic structure emerged through multiple waves of migration and extensive population interaction within the Avar Khaganate, rather than through simple replacement or single migration events. The identification of PC9 as a discriminating component represents a methodological advance applicable to other cases where populations share dominant ancestry components but differ in marginal genetic signatures.

The evidence for 50%–60% Slavic ancestry in modern Croatian populations, retained through at least three distinct migration waves between the 6th and 8th centuries CE, fundamentally revises our understanding of South Slavic ethnogenesis. Rather than the single migration event suggested by Constantine VII Porphyrogenitus’s account, the genetic data reveals a complex process of sequential population movements, each contributing distinct genetic and cultural elements to the emerging Croatian identity.

Our findings align with recent large-scale ancient DNA studies (Sarkissian et al., 2022), while providing novel regional insights. The 50%–60% Slavic ancestry identified in modern Croatians exceeds the 30%–60% range reported by Olalde et al. for general Balkan populations, likely due to our improved detection method using marginal principal components. This level of genetic turnover parallels other major demographic transformations (Török, 2022) in European prehistory.

Similar population replacement levels have been documented in other contexts: the Anglo-Saxon migration to Britain showed 40%–75% continental ancestry (Gretzinger et al., 2022), while the Hungarian Conquerors contributed 30%–40% ancestry to the Carpathian Basin (Maár et al., 2021). Unlike the strongly male-biased migrations documented in Bronze Age Iberia (Olalde et al., 2023) and Anglo-Saxon Britain, our evidence suggests more balanced Slavic migrations with both male and female contributors, consistent with family-group movements rather than military campaigns alone.

The Baltic ancestry gradient we document (71% in Lithuanians decreasing to 39%–51% in early Croatians) provides quantitative support for the linguistic homeland hypothesis, comparable to the steppe ancestry gradients that track Indo-European expansions (Haak et al., 2015; Allentoft et al., 2015). This pattern resembles the dilution of Yamnaya ancestry from 75% in Corded Ware populations to 30%–50% in modern Europeans (Lazaridis et al., 2022), suggesting similar processes of sequential admixture during population expansions.

The extraordinary genetic diversity (Supplementary Figure S3) documented in the Steppe Barbaricum between the 1st and fifth centuries CE reveals a previously unrecognized zone of intensive population interaction. The balanced representation of five distinct ancestral components—Germanic, Proto-Slavic, Sarmatian, Byzantine, and Balkan Iron Age—each contributing approximately 20% ancestry, has few parallels in ancient DNA studies.

This level of population integration exceeds even the cosmopolitan profiles documented in Roman frontier cities (Veeramah et al., 2018) or trading centers like Viking Age Birka (Krzewińska et al., 2018). The presence of Proto-Slavic components in fifth-century Barbaricum samples predates traditional historical accounts by over a century, suggesting either earlier migrations than documented or in situ development of Slavic identity through ethnogenesis processes similar to those proposed for other barbarian groups (Geary, 1999; Pohl, 2018).

The Barbaricum’s role as a demographic source contributing 20%–40% ancestry to multiple later populations parallels the function of other frontier zones in facilitating ethnic transformations. Similar processes have been documented in the Middle Danube region during the Bronze Age (Gerber et al., 2023) and in the Pontic-Caspian steppe during Scythian formations (Järve et al., 2019).

The identification of individuals with predominantly Slavic genetic profiles in elite Avar burials fundamentally challenges traditional interpretations of rigid ethnic hierarchies within the Khaganate. This finding parallels recent discoveries in other nomadic empires: the Xiongnu confederation included genetically diverse elites (Jeong et al., 2020), while the Mongol Empire incorporated local elites across its territories (Damgaard et al., 2018).

The substantial Byzantine genetic component (15%–20%) in Avar-period Slavic profiles indicates ongoing interaction with imperial territories despite political tensions. This pattern resembles the genetic continuity documented across the Roman-Medieval transition in other frontier regions (Amorim et al., 2018) and suggests that political boundaries did not prevent population movement and intermarriage.

The presence of individuals with predominantly Baltic ancestry (59%–75%) in Avar contexts, exemplified by sample CSB-9. SG, provides crucial evidence for understanding pre-expansion Slavic genetic structure. These individuals may represent specialized military recruits, similar to the Varangians in Byzantine service or the diverse warrior bands documented in Viking Age Scandinavia (Margaryan et al., 2020).

The multi-wave migration model emerging from genetic data aligns with fragmented historical accounts while revealing previously unknown complexity. The three distinct waves we identify - early (6th century, >50% Baltic ancestry), Avar-associated (6th-7th century, 25%–30% Baltic), and late (7th-8th century, 35%–40% Baltic) suggest that Croatian ethnogenesis involved sequential integration of related but distinct Slavic groups.

This pattern contrasts with single-event models of ethnic formation and instead supports processual models emphasizing gradual transformation. However, unlike purely constructivist approaches that minimize demographic change, our data demonstrates substantial population replacement accompanied by cultural transformation. The persistence of 20%–25% pre-Slavic Balkan ancestry indicates that Croatian formation involved integration rather than complete replacement, similar to patterns documented in Anglo-Saxon England (Gretzinger et al., 2022) and Viking Age Ireland (Margaryan et al., 2020).

The emergence of key South Slavic genetic features within the Avar Khaganate, rather than after Balkan settlement, suggests that interactions in the Carpathian Basin fundamentally shaped South Slavic identity. This finding supports recent archaeological arguments for the Carpathian Basin as a crucible of early medieval ethnic formations (Csaki et al., 2020) and highlights the importance of the Avar period in European demographic history.

Our systematic screening approach for marginal principal components addresses a fundamental challenge in population genetics: distinguishing closely related populations that share dominant ancestry sources. The success of PC9 in discriminating Slavic from Germanic populations, despite their shared Baltic Bronze Age heritage, demonstrates that informative genetic variation often resides in components capturing minimal overall variance.

This methodology has broad applicability for studying other cases of populations with shared major ancestry but distinct demographic histories.

The genetic evidence provides new context for interpreting historical sources and archaeological patterns. The presence of Slavic genetic signatures in the Barbaricum by the fifth century CE correlates with Jordanes’s mentions of “Sclaveni” north of the Danube (551 CE) and supports earlier presence than traditionally assumed. The identification of Slavs in elite Avar burials aligns with Byzantine sources mentioning Slavic leaders operating within the Avar confederation (Menander Protector, 6th century).

Archaeological patterns of material culture change in the 6th-7th centuries, particularly the spread of Prague-type pottery and sunken-floored dwellings, can now be understood as reflecting substantial demographic change rather than merely cultural diffusion. The genetic evidence for multiple migration waves correlates with the archaeological identification of distinct pottery traditions and burial customs appearing sequentially rather than simultaneously (Djino, 2010).

While our analysis provides unprecedented resolution of Croatian ethnogenesis, several limitations should be acknowledged. The systematic PC screening approach, while effective for this case, requires validation in other population contexts. The specific principal component capturing discriminating ancestry may vary depending on the populations included in analysis, necessitating systematic screening rather than assuming PC9 will universally capture Slavic ancestry.

Temporal resolution remains limited by the availability of well-dated samples from crucial transitional periods, particularly the fifth-6th centuries CE. Future studies incorporating more samples from this period could refine our understanding of the timing and tempo of population changes. Additionally, integrating uniparental markers (Y-chromosome and mitochondrial DNA) could provide insights into sex-biased migration patterns and social organization.

The complex demographic history revealed here suggests that other European populations traditionally viewed as resulting from simple migration events may similarly show evidence of multiple waves and extensive admixture when analyzed with appropriate methods. Future studies should apply similar systematic approaches to investigate ethnogenesis processes in other regions where historical and archaeological evidence suggests complex population formations.