The rise and transformation of Bronze Age pastoralists in the Caucasus

Abstract

The Caucasus and surrounding areas, with their rich metal resources, became a crucible of the Bronze Age¹ and the birthplace of the earliest steppe pastoralist societies². Yet, despite this region having a large influence on the subsequent development of Europe and Asia, questions remain regarding its hunter-gatherer past and its formation of expansionist mobile steppe societies^3-5. Here we present new genome-wide data for 131 individuals from 38 archaeological sites spanning 6,000 years. We find a strong genetic differentiation between populations north and south of the Caucasus mountains during the Mesolithic, with Eastern hunter-gatherer ancestry^4,6 in the north, and a distinct Caucasus hunter-gatherer ancestry⁷ with increasing East Anatolian farmer admixture in the south. During the subsequent Eneolithic period, we observe the formation of the characteristic West Eurasian steppe ancestry and heightened interaction between the mountain and steppe regions, facilitated by technological developments of the Maykop cultural complex⁸. By contrast, the peak of pastoralist activities and territorial expansions during the Early and Middle Bronze Age is characterized by long-term genetic stability. The Late Bronze Age marks another period of gene flow from multiple distinct sources that coincides with a decline of steppe cultures, followed by a transformation and absorption of the steppe ancestry into highland populations.

Fig. 1. Geographical and chronological overview of individuals analysed in this study.

a, Map of the wider Caucasus region showing the locations of the 38 sampled sites and 131 individuals from this study and published data (Supplementary Tables 1 and 5). The number in parentheses indicates the studied individuals per site. The inset shows a magnified view of the area marked by the rectangle in the main image. Scale bars, 100 km (main image) and 10 km (inset). b, Timeline of the median ¹⁴C ages of newly reported (black outline) and published individuals (no outline), separated on the x axis into main genetic and geographic clusters (Steppe, Intermediate and Caucasus), as well as other published individuals from south of the Greater Caucasus for comparison. Colours in the piecharts (a), and colours and symbols in the timeline (b) represent sampled individuals associated with different archaeological cultural complexes, and these are used consistently across all figures (see Figs. 2 and 3, Extended Data Fig. 1, Supplementary Table 1 and Supplementary Fig. 1 for full site names and individual symbols). Grey bars correspond to the chronological chapters in the results section. The map was generated using Base Relief: Mapzen, OpenStreetMap, and rivers, lakes and borders were added using free vector and raster map data from Natural Earth (https://www.naturalearthdata.com). OpenStreetMap is open data, licensed under the Open Data Commons Open Database Licence by the OpenStreetMap Foundation.

Fig. 2. Genetic overview of the seventh to fourth millennium bc.

a,b, PCA of newly produced ancient individuals (with outline) and individuals from previous publications (no outline) from the seventh and fifth millennium bc (a) and from the fourth millennium bc (b), projected onto 102 modern-day populations (grey dots). The dashed arrows in a represent the observed admixture clines between central Anatolian Neolithic and CHG (brown), between CHG and EHG (light green), and between Eneolithic_Caucasus and EHG ancestry-carrying Steppe groups (dark green), respectively. The dashed pink arrow in b represents an observed cline of mixture between Maykop-associated Caucasus groups and those carrying WSHG ancestry in the steppe. The corresponding labels and groupings are listed in Supplementary Table 5. c,d, Sankey diagram of genetic ancestry modelling for the seventh to fifth millennium bc (c) and fourth millennium bc (d) individuals from the Caucasus region with temporally and geographically proximal sources. The admixture proportions (as percentages) are indicated on each ancestry flow, with sources on the left and target populations on the right, and P values for each model in brackets under the population names (Supplementary Tables 8, 10 and 13). The suffixes in the group labels present archaeological time periods and geographical regions: WHG, Western hunter-gatherer; SHG, Scandinavian hunter-gatherer; MN, Middle Neolithic; C, Chalcolithic; EC, Early Chalcolithic; LC, Late Chalcolithic. Çatalhöyük_N is shown as Anatolia Neolithic, Tell Kurdu is shown as Anatolia EC and Jordan_PPNB and Jordan_PPNC are shown as Levant PPN.

Fig. 3. Genetic overview of the third and second millennium bc.

a,b, PCA of newly produced ancient individuals (black outline) and individuals from previous publications (no outline) from the third millennium bc (a), and the second millennium bc (b), projected onto 102 modern-day populations (grey dots). The dashed arrows represent observed mixture clines between the Caucasus and Steppe groups and re-emerging gene flow from the northeast. The correspondent labels and groupings are listed in Supplementary Table 5. c,d, Sankey diagram of genetic ancestry modelling for third millennium bc (c) and second millennium bc (d) individuals from this study based on temporally and geographically proximal sources. The admixture proportions (as percentages) are indicated on each ancestry flow, with sources on the left and target populations on the right, and P values for each model in brackets under the population names (Supplementary Tables 15 and 17). The suffixes in the group labels present archaeological time periods and geographical regions: MLBA, Middle–Late BA; BIA, Bronze–Iron Age; IA, Iron Age.

Extended Data Fig. 1. Principal component analysis (PCA) of prehistoric individuals from the Caucasus region.

a, extended West Eurasian PCA with 1522 individuals from 102 populations. b, West Eurasian PCA with 1243 individuals from 82 populations. The new individuals are shown with black outline.

Extended Data Fig. 2. Results of ADMIXTURE analysis (k = 10) of the individuals from the Caucasus.

a, Individuals belonging to the Steppe cluster, b, Individuals belonging to the Caucasus cluster, and c, a representative selection of published reference individuals. Newly genotyped individuals (black outline) and published individuals are sorted by archaeological or genetic groups in chronological order from left to right.

Extended Data Fig. 3. Outgroup f₃-statistics comparing genetic affinities of Mesolithic and Neolithic groups.

a, Heatmap showing shared genetic drift estimated by outgroup f₃-statistics of the form f₃(Mbuti; X, Y) between all EHG-related groups published to date. Lighter colours indicate higher f₃ values which corresponds to higher shared genetic drift. b, Scatterplot of outgroup f₃-statistics measuring shared drift between Neolithic and Chalcolithic groups to Anatolia and Levant PPN. Groups from the Caucasus cluster are shown in red and brown symbols. The trend line is given by the red dashed line. Error bars indicate ± 1 s.e. and were calculated using a weighted block jackknife across all autosomes on the 1,240,000 panel (nSNPs = 1,150,639) and a block size of 5 Mb.

Extended Data Fig. 4. Formal test for temporal genetic changes in the Caucasus region.

a, f₄-statistic tests where X denotes Anatolia N, EHG, or Iran_N and b, where X denotes CHG, EHG, and Iran_N, and Test denotes various genetic groups from the Caucasus and Steppe clusters on the y axis. Significant Z-scores (|Z| > 3) are highlighted in orange and error bars indicate ± 3 s.e. and were calculated using a weighted block jackknife across all autosomes on the 1,240,000 panel (nSNPs = 1,150,639) and a block size of 5 Mb.

Extended Data Fig. 5. Admixture modelling of Neolithic South Caucasus groups and Admixture dating.

Ternary plot of distal three-way qpAdm models for Neolithic and Chalcolithic groups from South Caucasus, Anatolia and Mesopotamia. a, Iran_N and b, CHG together with Anatolia_PPN, Levant_PPN are used as sources. c, Admixture date estimates calculated using the method DATES and two different pools of sources (Supplementary Table 11). Radiocarbon date intervals for each genetic group are shown as grey bars, circles and diamonds are the estimated admixture dates converted to calendar dates (assuming a generation time of 28 years) based on the mean and the oldest ¹⁴C dates, respectively. Error bars indicate ± 1 s.e.

Extended Data Fig. 6. Key f₄-statistics and qpAdm results of Eneolithic and Yamnaya groups.

a, Variation in genetic affinities of Eneolithic groups to EHG or CHG. Significant Z-scores (|Z| > 3) are highlighted in orange and error bars indicate ±1 (black) and ±3 (light gray) s.e. and were calculated using a weighted block jackknife across all autosomes on the 1,240,000 panel (nSNPs = 1,150,639) and a block size of 5 Mb. b, Ternary plot of qpAdm modelling for Steppe Eneolithic and Yamnaya groups (Supplementary Tables 10 and 15).

Extended Data Fig. 7. Results of genetic sex analysis and uniparentally inherited markers.

a, Total count of XX and XY individuals per genetic grouping from this study and Wang et al.. The chrono-cultural and genetic groups are chronologically ordered within the Steppe and Caucasus clusters. b, Y-chromosome haplogroups and c, mitochondrial haplogroups compared across the three main genetic clusters of Caucasus, Intermediate and Steppe. The x axis shows the total number of individuals per genetic cluster.

Extended Data Fig. 8. Kurgans as demographic snapshots and micro-transects through time.

a, and b, show the general layout of the Marinskaya 5 and Marfa kurgans (Bronze Age only). Analysed individuals are highlighted in different colours. c, and d, PCA plots highlighting the genetic heterogeneity of individuals buried in these kurgans. e, and f, Estimates of genetic relatedness calculated as mean pairwise mismatch rates (PMR) for individuals from Marinskaya 5 and Marfa, with error bars indicating ± 2 s.e. Dashed lines represent the expected PMR for each relatedness degree.

Extended Data Fig. 9. IBD analysis results per pair of individuals with ancIBD.

a, Visualising the sum (x axis) and the number (y axis) of all IBD tracts with a length of at least 12cM. Only pairs with at least two shared chunks of 20cM length are shown. b, Map showing IBD networks of closely related individuals within and between sites of the piedmont and steppe zone. The map was generated using Base Relief: Mapzen, OpenStreetMap, and rivers, lakes and borders were added using free vector and raster map data from Natural Earth (https://www.naturalearthdata.com). OpenStreetMap is open data, licensed under the Open Data Commons Open Database Licence by the OpenStreetMap Foundation.

Extended Data Fig. 10. Assessment of runs of homozygosity and inbreeding.

a, b, and c, Inferred run of homozygosity ROH per individual in Steppe and Caucasus clusters, respectively. Results are plotted by genetic and cultural groups in relative chronological order from left to right. b, Full scale of inferred ROH for individuals in the Maykop group. The legend illustrates the expected ROH for offspring of small populations or closely related individuals. d, Karyogram of VS5001. The positions of ROH longer than 4 cM are marked on the 22 autosomes (maroon).