The genomic history of the Aegean palatial civilizations

Abstract

The Cycladic, the Minoan, and the Helladic (Mycenaean) cultures define the Bronze Age (BA) of Greece. Urbanism, complex social structures, craft and agricultural specialization, and the earliest forms of writing characterize this iconic period. We sequenced six Early to Middle BA whole genomes, along with 11 mitochondrial genomes, sampled from the three BA cultures of the Aegean Sea. The Early BA (EBA) genomes are homogeneous and derive most of their ancestry from Neolithic Aegeans, contrary to earlier hypotheses that the Neolithic-EBA cultural transition was due to massive population turnover. EBA Aegeans were shaped by relatively small-scale migration from East of the Aegean, as evidenced by the Caucasus-related ancestry also detected in Anatolians. In contrast, Middle BA (MBA) individuals of northern Greece differ from EBA populations in showing ∼50% Pontic-Caspian Steppe-related ancestry, dated at ca. 2,600-2,000 BCE. Such gene flow events during the MBA contributed toward shaping present-day Greek genomes.

Keywords: Anatolia; Bronze Age; Cycladic civilization; Greece; Helladic civilization; Minoan civilization; Mycenean civilization; ancient DNA; paleogenomics; population genetics.

Graphical abstract

Figure 1

Geographical location of archeological sites and radiocarbon dates (A) Archeological sites are indicated by square symbols, colored according to their corresponding Aegean cultural group. Asterisks indicate archeological sites for which mitochondrial DNA (mtDNA) capture data were generated (STAR Methods). The whole genomes of six individual samples (Pta08, Kou01, Kou03, Mik15, Log02, and Log04) sequenced are colored according to their radiocarbon dates. (B) Radiocarbon dates of the six whole-genome sequenced individuals from this study together with 12 individuals from present-day Greece from previous studies (Lazaridis et al., 2017, Mathieson et al., 2018, Hofmanová et al., 2016). For the two Mesolithic individuals, only mtDNA data is available (Hofmanová et al., 2016). The bar indicates the range of the Cal 1-sigma OxCal calibrated date for each individual (STAR Methods). See also Figures S1 and S3, Table 2, and Document S1.

Figure S1

Images of archaeological site Elati-Logkas (Log02, Log04), related to Table 2 and Document S1 (A) Elati-Logkas, view of the cemetery with burials covered by stones known as “periboloi.” (B) Elati-Logkas, Burial 80.1 (Log04) is a pit-grave in the circumference of an inner enclosure built from rough stones. The buried individual was in crouched position lying to the left side, with the hands bent and the palms supporting the skull. Inside the same walls, four other similar burials were excavated with no grave goods apart from only one flint stone blade in tomb 80.5. (C) Elati-Logkas, Burial 22.1 (Log02) is the main among three pithos-inhumations and one secondary burial inside the “peribolos 22.” The grave itself is bordered by rough stones, with the buried individual laid on a ceramic “stretcher.” Several vertical lines are still visible on the skeletal remains. The individual of the burial 22.1 was found in a supine position with the hands crossed on the abdomen, the legs bent in a crouched position to the left, and the skull turned to the right side. There were no grave goods found in the burial 22.1. Photo credits: Ephorate of Antiquities of Kozani, Hellenic Ministry of Culture, Greece. Courtesy of Dr. Georgia Karamitrou-Mentessidi.

Figure S2

Comparison of SNP capture and WGS data, related to Figure 4 and Tables 2 and S1 (A) Number of single nucleotide polymorphisms (SNPs) in the nuclear genomic data from this study (WGS data) in comparison with previously published BA genomic data from the Aegean. The number of covered SNPs across BA Aegeans based on two SNP sets are shown. On the left: the number of SNPs based on the 1240K SNP set (Dataset I) defined by the array used in Lazaridis et al. (2017) to enrich the libraries. On the right: the number of SNPs based on the intergenic regions defined for the ABC-DL analysis below (Dataset IV, STAR Methods). The green box plots (median indicated by a horizontal line and interquartile range indicated by the box) correspond to the number of SNPs among the BA Aegean data from present-day Greece (Lazaridis et al., 2017); the blue box plots correspond to the number of SNPs among the whole genome sequence (WGS) data from this study. (B) One-dimensional Site Frequency Spectrum (SFS) for the seven whole genomes used for demographic analyses (ABC-DL). The seven genomes included here are: Mik15 and Log04 from this study, YamnayaKaragash_ EBA (3,018-2,887 BCE) (de Barros Damgaard et al., 2018), KK1 (CHG; 7,745-7,579 BCE) (Jones et al., 2015), Bar8 (Neolithic Barçın; 6,122-6,030 BCE) (Hofmanová et al., 2016), Sidelkino (EHG; 9,386-9,231 BCE) (de Barros Damgaard et al., 2018), and S_Greek-1 (SAME3302732; modern Greek from Thessaloniki) (Mallick et al., 2016).STAR Methods In blue ("WGS") the SFS for the regions included in Dataset IV (STAR Methods). In red ("1240K") the SFS for the regions in Dataset IV restricted to the sites overlapping with the SNPs included in the 1240K array.

Figure S3

Error rates, damage, and read length distributions for the WGS and nuclear capture data from this study, related to Figure 1 and Tables 2 and S1 (A) Error rate for whole genome sequencing before (lighter colors) and after (darker colors) trimming 5 bp from the extremities of the reads. Log02 was USER™-treated. (B) Error rate for nuclear capture data for different mutation types. Columns 1 and 2 show transitions and column 3 shows transversions. (C) Read length distribution for whole genome sequencing. (D) Read length distribution for nuclear capture data. (E) Post-mortem damage pattern for whole genome sequencing (C to T and G to A substitutions). Dashed lines indicate partial data removal resulting from trimming 5 bp from the extremities of the reads. The color of each curve indicates the analyzed sample according to panel A. Log02 (dark green curve) was USER™-treated. (F) Post-mortem damage pattern for nuclear capture data. Curves are colored according to panel B. See STAR Methods for details.

Figure S4

Genetic affinities between Neolithic, BA, and present-day Aegeans compared to other present-day Eurasian populations, related to Table 1 f₃-statistics of the form f₃(Yoruba; Y, X): Y corresponds to either Neolithic Anatolians or Greeks, Minoan-Petras-EBA individual from the island of Crete (Pta08), the Cycladic-Koufounisi-EBA individuals (Kou01, Kou03), the Helladic-Manika-EBA individual from the island of Euboea (Mik15), the Helladic-Logkas-MBA individuals from northern Greece (Log02, Log04), previously published BA Aegeans (Mycenaeans and Minoans), and present-day Greeks (incl. Cretans), Cypriots, while X are other present-day populations from Dataset I (Lazaridis et al., 2014, 2016, 2017) (STAR Methods). For clarity, we only show results for west Eurasian and north African populations and cap f₃ values below 0.15. For each case, we show the geographic distribution of f₃ (warmer colors represent greater sharing between populations X and Y). Beside each map, we plot the f₃ values for the 15 populations that are most closely related to each of the populations in Y (bars represent ~1.95 standard errors). In agreement with MDS and ADMIXTURE analyses, we observed that ancient and present-day Anatolians and Greeks share the most genetic drift with present-day central and southern European populations.

Figure 2

Multidimensional scaling analysis (MDS) Included are the six Aegean Bronze Age individual samples from the present study, 259 ancient genomes, and 638 modern individuals (gray shapes) of Eurasian ancestry (Table S2; STAR Methods). Population labels are given in Table 1. See also Table S5 and Document S1.

Figure 3

ADMIXTURE analysis for modern and ancient Eurasian individuals Shown are a subset of individuals for K = 3, which has the lowest cross validation error (CV = 0.974). Results for the full dataset and statistical support are shown in Figure S5. The six BA individuals whole genome sequenced in this study are highlighted with an asterisk. Abbreviations for chronological periods and population names are given in Table 1. See also Table S2, Document S1, and STAR Methods.

Figure S5

ADMIXTURE analysis using ancient and modern populations with the number of ancestry components ranging from K = 2–6 and cross-validation error, related to Figure 3, Table 1, and Document S1 (A) For this analysis we consider a total of 969 individuals (638 modern and 331 ancient) and 165,447 SNPs (Dataset II, STAR Methods). Each bar represents one individual. Individuals from the same population were grouped. For all K > 2, red represents the component mostly present in “European Neolithic-like,” light blue in “Neolithic Iran/Caucasus HG-like” and orange for “European HG-like.” (B) Cross-Validation error (CV error) for K ranging from 2 to 6. The CV-error is plotted for the ten runs for each value of K (STAR Methods).

Figure S6

Exploring differential allele-sharing in Aegean/Anatolian populations through time with D-statistics, related to Table 1 D-statistics of the form D(Anatolia_N, H2; H3, Mota) were computed, testing whether Anatolia_N (H1), or ancient/modern Anatolian and Aegeans (H2) share more alleles with CHG, Iran_N, EHG, or Steppe_EMBA (H3, STAR Methods). For this analysis, the genome of an Ethiopian individual (Mota) was used as an outgroup. Points represent D-statistics, and horizontal error bars represent ~3.3 standard errors (SE corresponding to a p-value of ~0.001 in a Z-test). Vertical bars represent upper and lower bounds of the dates available for the populations. In this figure, the populations are ordered chronologically, using either radiocarbon dates (when available) or dated archaeological context. Horizontal dashed lines indicate time periods. Vertical dashed lines mark the zero. A value of D = 0 indicates no gene flow or ancestral population structure (Durand et al., 2011), thus H1 and H2 are symmetrically related to H3 and Mota. In this case, D < 0 would indicate potential gene flow between H3 and H2, and D > 0 would indicate potential gene flow between H3 and H1. Abbreviations for chronological periods and population names are given in Table 1.

Figure 4

Model comparison in ABC-DL analysis (A) Posterior probabilities P(M|D) of different 3-leaf models (models A1–A4) calculated with ABC-DL to establish the topology of the three ancestral populations: CHG, EHG and Aegean Neolithic. (B) Seven-leaf demographic models (models B1–B6) extending the tree from (A) with the highest posterior probability, each including Neolithic, EBA and MBA Aegeans, present-day Greeks, CHG, EHG, and Pontic-Caspian Steppe_EMBA populations. Yellow (EHG-like), light blue (CHG-like), and dark blue (Steppe-like) arrows indicate a single pulse of gene flow from simulated “ghost” populations diverged from EHG, CHG, and Steppe_EMBA. Posterior probabilities are listed below each schematic topology (Document S1; STAR Methods). See also Figure S2 and Tables 1 and S4.

Figure 5

Sex-biased gene flow Comparison of X-linked and autosomal genetic ancestries associated with (A) Iran_N/CHG-like and (B) Steppe-like components in EBA and MBA Aegean ancient genomes. Violins show the distribution of point estimates across 100 replicates for the corresponding autosomal ancestry, with median indicated by a dot and interquartile ranges indicated by boxes. For the violin plots, we considered a random set of autosomal SNPs matching X-linked SNPs in number (i.e., 8,133) (STAR Methods). Orange diamonds show the point estimates and one associated Standard Error for the same ancestries on the X chromosome. See also Figure 3 and Table 1.

Figure S7

Estimated total ROH length by size category for six ancient and four modern genomes, related to Document S1 (A) ROHs estimated from 43 million imputed transitions and transversions, and (B) 13 million imputed transversions (STAR Methods).