Proceedings of the Royal Society B: Biological Sciences
You have accessResearch article

Archaeogenomic analysis of the first steps of Neolithization in Anatolia and the Aegean

Gülşah Merve Kılınç

Gülşah Merve Kılınç

Department of Archaeology and Classical Studies, Stockholm University, Lilla Frescativaegen 7, Stockholm 114 18, Sweden

[email protected]

Google Scholar

Find this author on PubMed

,
Dilek Koptekin

Dilek Koptekin

Department of Health Informatics, Middle East Technical University, Ankara 06800, Turkey

Google Scholar

Find this author on PubMed

,
Çiğdem Atakuman

Çiğdem Atakuman

Department of Settlement Archaeology, Middle East Technical University, Ankara 06800, Turkey

Google Scholar

Find this author on PubMed

,
Arev Pelin Sümer

Arev Pelin Sümer

Department of Biological Sciences, Middle East Technical University, Ankara 06800, Turkey

Google Scholar

Find this author on PubMed

,
Handan Melike Dönertaş

Handan Melike Dönertaş

European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Trust Genome Campus, Cambridge, Hinxton CB10 1SD, UK

Google Scholar

Find this author on PubMed

,
Reyhan Yaka

Reyhan Yaka

Department of Biological Sciences, Middle East Technical University, Ankara 06800, Turkey

Google Scholar

Find this author on PubMed

,
Cemal Can Bilgin

Cemal Can Bilgin

Department of Biological Sciences, Middle East Technical University, Ankara 06800, Turkey

Google Scholar

Find this author on PubMed

,
Ali Metin Büyükkarakaya

Ali Metin Büyükkarakaya

Department of Anthropology, Hacettepe University, Ankara, Beytepe 06800, Turkey

Google Scholar

Find this author on PubMed

,
Douglas Baird

Douglas Baird

Department of Archaeology, Classics, and Egyptology, University of Liverpool, Liverpool L69 7WZ, UK

Google Scholar

Find this author on PubMed

,
Ezgi Altınışık

Ezgi Altınışık

Department of Biology and Ecology, Faculty of Science, University of Ostrava, Ostrava, Czech Republic

Google Scholar

Find this author on PubMed

,
Pavel Flegontov

Pavel Flegontov

Department of Biology and Ecology, Faculty of Science, University of Ostrava, Ostrava, Czech Republic

A.A. Kharkevich Institute for Information Transmission Problems, Russian Academy of Sciences, Moscow, Russia

Institute of Parasitology, Biology Centre, Czech Academy of Sciences, České Budějovice, Czech Republic

Google Scholar

Find this author on PubMed

,
Anders Götherström

Anders Götherström

Department of Archaeology and Classical Studies, Stockholm University, Lilla Frescativaegen 7, Stockholm 114 18, Sweden

[email protected]

Google Scholar

Find this author on PubMed

,
İnci Togan

İnci Togan

Department of Biological Sciences, Middle East Technical University, Ankara 06800, Turkey

[email protected]

Google Scholar

Find this author on PubMed

and
Mehmet Somel

Mehmet Somel

Department of Biological Sciences, Middle East Technical University, Ankara 06800, Turkey

[email protected]

Google Scholar

Find this author on PubMed

Published:https://doi.org/10.1098/rspb.2017.2064

    Abstract

    The Neolithic transition in west Eurasia occurred in two main steps: the gradual development of sedentism and plant cultivation in the Near East and the subsequent spread of Neolithic cultures into the Aegean and across Europe after 7000 cal BCE. Here, we use published ancient genomes to investigate gene flow events in west Eurasia during the Neolithic transition. We confirm that the Early Neolithic central Anatolians in the ninth millennium BCE were probably descendants of local hunter–gatherers, rather than immigrants from the Levant or Iran. We further study the emergence of post-7000 cal BCE north Aegean Neolithic communities. Although Aegean farmers have frequently been assumed to be colonists originating from either central Anatolia or from the Levant, our findings raise alternative possibilities: north Aegean Neolithic populations may have been the product of multiple westward migrations, including south Anatolian emigrants, or they may have been descendants of local Aegean Mesolithic groups who adopted farming. These scenarios are consistent with the diversity of material cultures among Aegean Neolithic communities and the inheritance of local forager know-how. The demographic and cultural dynamics behind the earliest spread of Neolithic culture in the Aegean could therefore be distinct from the subsequent Neolithization of mainland Europe.

    1. Introduction

    The primary zone of Neolithization in western Eurasia encompassed the Levant, Taurus-Zagros ranges of Mesopotamia, central Anatolia and Cyprus [14]. The earliest evidence for sedentary life and food storage in this region goes back to the Natufians (ca 12 500–10 800 cal BCE) [5,6]. Sedentary communities were established across this zone during the first phase of the Pre-Pottery Neolithic (PPN, or Aceramic Neolithic, ca 10 000–8500 cal BCE), and the first indications of plant cultivation appeared [79]. Between ca 8500 and 7000 cal BCE, community sizes increased, architectural elaboration intensified and a subsistence economy based on agriculture gradually became the norm [1014]. Meanwhile, portable artefacts such as figurines and stamps evolved into staples of sedentary life, and pottery production became widespread ca 7000 cal BCE [10,11]. The elements of the subsequent Pottery Neolithic culture (PN, ca 7000–5500 cal BCE), including integrated cultivation practices of domestic plants and animals, the architectural practices of sedentary life, together with portable artefacts, have been collectively described as the Near Eastern ‘Neolithic Package’ [1518].

    During the same period, there were no signs of a Neolithization in west Anatolia and the Aegean. Only after ca 7000 cal BCE did elements of the ‘Neolithic Package’ appear in these regions, eventually spreading towards Europe [1921]. Some archaeologists suggest that the emergence of the Neolithic elements in the Aegean and in Europe without a preceding PPN development period indicates the role of demic processes (i.e. migrations from the Neolithic primary zones through land and sea routes), frequently described as a leap-frog model where migrants form enclaves in new territory [15,16,2227]. Others, in contrast, favour a role for interaction between local foragers and primary zone Neolithic populations, including the adoption of Neolithic elements by locals and acculturation [16,18,28,29].

    Recent archaeogenomic data have shown that the Neolithization of central, western and northern Europe involved migration from a single eastern source, frequently termed ‘Anatolian farmer’ [3033], while in other regions, such as in the Baltic [33,34] and in southern Greece [33], acculturation may have played role. In most of Europe, there is limited genetic evidence for early admixture between farmers and local European Mesolithic (WHG) communities in the sixth millennium BCE, such that early European farmers studied to date (with few exceptions [35]) carry ancestry similar to farmers from northwest Anatolian Barcın [3133] (electronic supplementary material, figure S1). In subsequent millenia, however, WHG-like ancestry appears in Middle and Late Neolithic European populations [3638]. These observations support a leap-frog model of Neolithic spread in Europe [28]: farmers only occupied enclaves in the new territories while Mesolithic groups persisted in the same regions [3942].

    The processes behind the earliest steps of Neolithization and the Neolithic spread in the Aegean are less understood. For instance, whether Aegean Neolithic populations were recent colonists originating from areas of the primary zone of Neolithization (e.g. [27]), descendants of indigenous foragers (e.g. [20]) or admixed groups (e.g. [15,16]) is still contentious. Additionally, whether Aegeans' demographic or cultural relationships were stronger with central Anatolians [16] or with the Levantine seafaring populations [27] remains unclear. We re-analyse published ancient human genomes to answer these questions and to dissect the demographic dynamics behind the Neolithic transition in Anatolia and the Aegean.

    2. Material and methods

    (a) Compiling and mapping genomic data

    We obtained DNA sequencing data of 99 published ancient individuals (electronic supplementary material, table S1), generated using whole-genome shotgun sequencing and/or sequencing of libraries enriched by hybridization capture [3032,35,37,38,4246]. We mapped sequencing reads to the human reference genome (hs37d5) using the Burrows–Wheeler Aligner (BWA, v. 0.7.12) [47], with the parameters ‘-l 16500, -n 0.01, -o 2’. We filtered PCR duplicates using FilterUniqSAMCons.py [48]. We filtered reads shorter than 35 base pairs, with greater than 10% mismatches to the reference, and less than 30 mapping quality per read.

    (b) Preparation of population genetics analysis datasets

    We restricted our analysis to known present-day DNA variants to minimize false positives. We used two different modern reference panels, calling genotypes of ancient individuals for SNPs overlapping with (i) the Human Origins genotype dataset [42,49] and (ii) the 1000 Genomes whole-genome sequence data [50] using SAMtools mpileup (v. 1.3) [51]. For (i) we obtained a curated version of the Human Origins panel of 594 924 autosomal SNP genotype calls for 2730 present-day individuals from Lazaridis et al. [42]. We determined the SNPs of the ancient samples overlapping with this dataset. We encoded transitions as missing to avoid confounding with cytosine deamination in ancient DNA. To prepare (ii), we obtained the BAM and VCF files for the African Yoruba individuals from 1000 Genomes Project phase 3 from McVean et al. [50]. Using vcftools [52], we extracted a total of 1 938 919 transversion SNPs with minor allele frequencies of greater than or equal to 10% in the Yoruba population to avoid false positive calls [36,40]. We determined the positions in the ancient samples overlapping with this dataset. We merged ancient genotypes with these two datasets using PLINK [53] requiring base quality greater than or equal to 30 per overlapping position. We haploidized each full dataset by randomly selecting one allele per position. The Human Origins-merged dataset, which has higher number of present-day populations, was used for principal component analysis (PCA) and for calculating f3-statistics. The 1000 Genomes-merged dataset, with a higher number of SNPs, was used for D-statistics, where we require high statistical power.

    (c) Principal component analysis

    We performed PCA by calculating principal components using west Eurasian populations from the Human Origins dataset using the smartpca program of EIGENSOFT [54], with the ‘numoutlieriter:0’ parameter. We projected ancient genomes onto the reference space using the ‘lsqproject:YES’ option and plotted the results using R (v. 3.3.0).

    (d) D- and f3-statistics

    We computed D-statistics using the qpDstat program of the ADMIXTOOLS package [49]. We assessed statistical significance by calculating standard errors using a block jacknife of 0.5 Mbp. We used the Yoruba population as outgroup for the D-statistics [32]. We computed f3-statistics (i.e. genetic affinity between pairs of populations based on an estimate of shared drift between them as their divergence from an outgroup population) using the qp3Pop program of the ADMIXTOOLS package [49]. The Human Origins dataset's African Mbuti population was used as outgroup for calculating f3-statistics [42]. We performed multiple testing correction using Benjamini–Yekutieli method for all 207 D-statistics results and reported adjusted p-values together with Z scores per each test [55]. For the pairwise f3-statistics, as genetic distance measure between a pair of populations, X and Y, we used: 1-f3(Mbuti;X,Y) [30]. These pairwise distances were summarized with the multidimensional scaling (MDS) method using the cmdscale function of R. We evaluated the goodness of fit for MDS using ‘GOF’ component obtained from cmdscale function.

    (e) Heterozygosity estimates

    We calculated heterozygosity as a measure of genetic diversity in a population, using genome sequence data of (i) Bon002 (from Boncuklu, central Anatolia, pre-7000 cal BCE) [30], Tep003 (Tepecik-Çiftlik, central Anatolia, post-7000 cal BCE) [30], Bar8 (Barcın, north Aegean, post-7000 cal BCE) [31] and Rev5 (Revenia, north Aegean, post-7000 cal BCE) [31]. We calculated genome coverage per sample using GenomeCoverageBed [56]. We downsampled the genome sequences of Bon002 and Bar8 to similar levels as the other two samples using SAMtools (v. 1.3) [51]. We calculated heterozygosity per sample using ANGSD [57] as ‘angsd -GL 1 -doGlf 2 -doMajorMinor 1 -sites ReferenceSNP.pos -bam bamlist -doSaf 1 -anc referencegenome.fasta’. To minimize false positives, we only considered transversions overlapping with of Yoruba individuals from 1000 Genomes Project phase 3 from McVean et al. [50].

    (f) Modelling of admixture

    We used the qpWave/qpAdm framework [38,58] in the ADMIXTOOLS package [49] to model populations as mixtures of two or more sources. The following worldwide set of ancient and present-day outgroups, which most probably did not experience any post-split gene flow from Anatolian/Aegean populations, was used: Mbuti, Yoruba, Ust Ishim, El Miron, Goyet Q116, Villabruna, Kostenki14, Vestonice16, Papuan, Onge, Karitiana, Mixe, Chipewyan, Oroqen, Koryak, Dai, Japanese. Adding East European hunter–gatherers (EHG) as a close outgroup to increase the resolution did not change the results.

    (g) Serial coalescent simulations

    We performed serial coalescent simulations using fastsimcoal [59] under four various demographic models involving Neolithic central Anatolians, Aegeans, Iranians and WHG (not including the Levantine populations, for whom we lack whole-genome data). The simulations were designed to mimick the data with respect to tree topology, divergence times and sample sizes. We then performed D-statistics on the simulated DNA and compared these with the observed data to gain understanding into the plausibility of different models. Specifically, we generated data to represent Iranian Neolithics (10 000 BP), WHGs (Loschbour: 7200 BP), central Anatolian Neolithics (Tepecik-Çiftlik: 8500 BP; Boncuklu: 10 000 BP), the Aegean Neolithics (Revenia: 8300 BP) and present-day sub-Saharan Africans (Yoruba-YRI). We launched 100 runs for each model defined in the parameter file (input.par) for testing different population histories. For all models, we sampled 30 Mb DNA sequences for: five present-day Yoruba, two Iranian Neolithics, two WHGs, four central Anatolian Neolithics (two Tepecik-Çiftlik, two Boncuklu) and two Aegean Neolithics (Revenia). We assumed a mutation rate of 1.00 × 10−9 bp yr−1, and a recombination rate of 1.00 × 10−8 bp yr−1, and assumed 25 years per generation, again following [45]. We set the effective population size (Ne) of these populations and times of divergence between Anatolian Neolithic, WHGs and Iranian Neolithic populations based on [45]. We converted all outputs (arp file) to plink format and computed D-statistics with topology of D(YRI, Test, central Anatolian N, Aegean N) to test the relationships among populations via ADMIXTOOLS [49]. Note that the tree topology involving the Anatolian/Aegean populations, Iran, WHG and the Africans, were based on the phylogenetic analysis from Broushaki et al. [45]. The Anatolian/Aegean populations were assumed to diverge simultaneously from the same source (star shaped).

    3. Results

    (a) Early Holocene gene pools of west Eurasia and the Anatolian/Aegean gene pool

    We compiled published genome sequence data of 99 ancient individuals (sample ages: ca 11 840–4360 cal BCE) (electronic supplementary material, figure S2, table S1). Both a PCA using present-day and ancient populations (electronic supplementary material, figure S1) and an MDS analysis using only ancient genomes (figures 1 and 2; electronic supplementary material, table S2) revealed the presence of four distinct gene pools in Early Holocene west Eurasia: (a) a ‘Caucasia/Iran gene pool’, (b) a ‘Levant gene pool’, (c) a ‘European pre-Neolithic gene pool’ and (d) an ‘Anatolian/Aegean gene pool’. To objectively measure clustering in gene pools (a)–(c), we used D-statistics of the form D(Yoruba, p1; p2, p3) where ‘p’ refers to the Caucasia/Iran, the Levant or European pre-Neolithic gene pools, correcting for multiple testing. In 80% comparisons (p < 0.05; Z ≥ 3), populations belonging to the same gene pool shared more alleles with each other compared to external populations (figure 1a; electronic supplementary material, figure S3a, table S3). The only exceptions were comparisons involving a single pre-Neolithic individual from Iran for which we had relatively few SNPs and low statistical power.

    Figure 1.

    Figure 1. Genetic differentiation among ancient west Eurasians and predicted admixture events. (a,c) Results of the same multidimensional scaling (MDS) analysis, summarizing f3-statistics (shared genetic drift) between ancient population pairs (electronic supplementary material, table S2). The goodness of fit was estimated as 0.17 and 0.17 for both dimensions. Admixture events among gene pools inferred using D-statistics are represented as arrows on each MDS plot. The circles where the arrow tips touch indicate which population is involved in the inferred admixture. Tepecik-Çiftlik is labelled as Tepecik. (a) Admixture in Boncuklu (central Anatolian PPN). For clarity, the other Anatolian/Aegean populations are not plotted. Arrow ‘a’: gene flow between Boncuklu and pre-Neolithic populations of mainland Europe (relative to other gene pools). Arrow ‘b’: gene flow between Boncuklu and the Levant populations (relative to other gene pools). Arrow ‘c’: gene flow between Boncuklu and Caucasia/Iran populations (relative to other gene pools). (b) Results of D-statistics in the form of D(Yoruba, p1; p2, Boncuklu). Multiple testing correction was performed using the Benjamini–Yekutieli method [55]. (c) Arrow ‘d’: estimated gene flow between Boncuklu and the Levantine populations. This is based on testing the topology D(Outgroup, Boncuklu; Levant_preNeolithic, Levant_Neolithic), showing that the Boncuklu population showed higher genetic affinity to the Levantine Neolithics (sample ages: ca 8300–6750 cal BCE) than to the Levantine pre-Neolithics (sample ages: ca 11 840–9760 cal BCE), although the result was only marginally significant (p > 0.05, Z > 2.5) (electronic supplementary material, figure S2c, table S5). Arrow ‘e’: gene flow from the Iran PPN population into Anatolian/Aegean PN populations. Arrow ‘f’: gene flow from the Levant PPN population into Anatolian/Aegean PN populations (electronic supplementary material, table S5 and table S6).

    Figure 2.

    Figure 2. Summary of D-statistics describing population relationships within the Anatolian/Aegean gene pool and between Anatolians/Aegeans and neighbouring groups. The Yoruba genome was used as outgroup in D-statistics. The Tepecik-Çiftlik is labelled as Tepecik. All D-statistics results are reported in electronic supplementary material, tables S8–S13. (a) D-statistics results summarized as arrows on the MDS plot (same as figure 1). Each triple population compared in D-tests is framed in the same colour. If a test population has greater genetic affinity to the second population compared to a third one, an arrow with same colour as the frames is drawn from the test population to the second population (the arrows' direction or lengths are not representative of gene flow magnitudes). Arrow ‘a’ and navy frames summarize D(Yoruba, Natufian; centralAnatolian, northAegean), where Natufians had stronger genetic affinity to north Aegean PN than to central Anatolian PPN or PN groups (electronic supplementary material, table S8). Arrows ‘b’ and ‘c’ and green frames summarize D(Yoruba, CHG&Iran_Neolithic; northAegean, centralAnatolian). In six of eight comparisons CHGs and Iran PPN populations had stronger genetic affinity to the north Aegean PN than to central Anatolian PPN and PN (electronic supplementary material, table S9). Arrow ‘d’ and purple frames summarize D(Yoruba, WHG; northAegean, centralAnatolian). In all comparisons WHGs had stronger genetic affinity to the north Aegean PN than to central Anatolian PPN and PN, with the exception of D(Yoruba, WHG; Boncuklu, Barcın) being non-significant (electronic supplementary material, table S10). (b) Results of D-tests calculated as D(Outgroup, RightPopulation; BottomPopulation, LeftPopulation), where right, bottom and left refer to the positions of the populations on the matrix. For instance the top row shows that Boncuklu has significantly higher affinity to Barcın than to Tepecik-Çiftlik. The D-statistic magnitude is represented by colour, Z score by size, and significance by being filled or not. Multiple testing correction was performed using the Benjamini–Yekutieli method [55] (electronic supplementary material, table S12 and S13).

    We then investigated the relationships among ancient Anatolians and other west Eurasian gene pools, using the oldest Anatolian population yet sequenced: Boncuklu from central Anatolia (sample ages: ca 8300–7952 cal BCE), an Aceramic Neolithic population previously predicted to be the descendants of local Epipalaeolithic groups [30,60]. We computed D-statistics of the form D(Yoruba, p1; p2, Boncuklu), where ‘p1’ and ‘p2’ refer to populations belonging to different gene pools: Caucasia/Iran, the Levant or the European pre-Neolithic. In 56% of the comparisons (p < 0.05; Z ≥ 2.8), all three regional gene pools showed higher affinity to Boncuklu than to each other (figure 1a,b; electronic supplementary material, table S4). Using the qpWave/qpAdm algorithm [38,58], we further modelled the Boncuklu population as a mixture of CHG (59.1%), the Levant (31.4%) and WHG (9.5%) (electronic supplementary material, table S5).

    We next included three post-7000 cal BCE Neolithic populations from Anatolia and Aegean in the analyses: Tepecik-Çiftlik in central Anatolia [30,61], Barcın in northwest Anatolia [31,37,62] and ‘Revenia’ in Pieria of northeast Greece [31]. We computed D-statistics of the form D(Yoruba, p1; p3, p2) where ‘p1’ and ‘p2’ are Anatolian/Aegean populations, and ‘p3’ is an external population (Caucasia/Iran, the Levant or European pre-Neolithic). In 94% of the comparisons (p < 0.05; Z ≥ 2.8) all Anatolian/Aegean populations were genetically closer to each other than to any other gene pool (electronic supplementary material, figure S3b, table S5).

    Given archaeological indication that Aegean Neolithic was influenced by east Mediterranean sources [27], we further studied the genetic affinities of Aegean Neolithic people to central Anatolian Neolithics and to the Levantines. Calculating D-statistics of the form D(Yoruba, northAegean; Levant, centralAnatolia) revealed that the post-7000 cal BCE Neolithic north Aegean individuals (Barcın and Revenia) consistently share more alleles with central Anatolians compared to south Levantines, where 50% of the comparisons were significant (p < 0.05; Z ≥ 2.8) (electronic supplementary material, figure S3b, table S5).

    (b) Notable genetic diversity in the Aegean

    To assess demographic events in the Near East during the Neolithic transition, we studied signatures of regional admixture using diachronic populations from the same region (figure 1c). In 83% of the comparisons, pre-7000 cal BCE Neolithic populations of the Levant and of Iran were genetically closer to all post-7000 cal BCE Anatolian/Aegean populations (Tepecik, Barcın, Revenia) compared with the pre-7000 cal BCE Anatolian Boncuklu (p < 0.05; Z ≥ 3) (electronic supplementary material, figure S3c, table S6). Considering the radiocarbon dates of the investigated individuals, this is consistent with gene flow from both the Levant and from Iran into Anatolia, within a period ranging from the PPN to the PN (figure 1c, arrows ‘e’ and ‘f’). These results are also compatible with a regional increase in the levels of admixture during the Neolithic [30,32], although alternative explanations to gene flow remain plausible, such as population structure confounding the analysis results [63].

    Next, to gain understanding into Aegean Neolithization, we studied the population genetic characteristics of the PN Aegean groups relative to central Anatolian groups. We first compared heterozygosity estimates among these populations. If the Aegeans were recent colonists from a single origin, due to a founder effect, one might expect lower heterozygosity in the Aegean than in central Anatolia. By contrast, Barcın and Revenia individuals had higher heterozygosity levels (mean 0.25 and 0.26, respectively) than those of Boncuklu and Tepecik (0.22 and 0.19, respectively) (electronic supplementary material, table S7).

    Second, we calculated D-statistics focused on the Aegeans, which suggested higher admixture in this region than in central Anatolia:

    • (i) D(Yoruba, Natufian; centralAnatolia, northAegean) revealed that pre-Neolithic population of the Levant had stronger genetic affinity to the two north Aegean Neolithic populations (Barcın and Revenia, post-7000 cal BCE) than to the two central Anatolian Neolithic groups (Boncuklu and Tepecik, pre- and post-7000 cal BCE) (p < 0.05; Z ≥ 3) (figure 2a, arrow ‘a’; electronic supplementary material, figure S4a, table S8). Given the above-proposed gene flow event from the Levant into Anatolia during the Neolithic, this result might imply additional genetic interactions between Natufian-related populations and the ancestors of north Aegean populations that bypassed central Anatolia.

    • (ii) D(Yoruba, Caucasia/Iran; centralAnatolia, northAegean) revealed that in 50% of the comparisons CHGs and Neolithic Iran individuals shared more alleles with the two north Aegean PN populations than with the two central Anatolians (p < 0.05; Z ≥ 2.8) (figure 2a, arrows ‘b’ and ‘c’; electronic supplementary material, figure S4a, table S9).

    • (iii) Likewise, WHG individuals showed higher affinity to the two north Aegean PN populations than PN central Anatolian group groups (p < 0.05; Z ≥ 3) (figure 2a, arrow ‘d’; electronic supplementary material, figure S4a, table S10).

    • (iv) Natufians, WHGs and Iranian PPN individuals were consistently more similar to the Revenia individual than to those in Barcın (electronic supplementary material, figure S4b, table S11).

    • (v) Both the Boncuklu (PPN) and the Tepecik (PN) groups of central Anatolia had stronger affinity to the north Aegean PN populations, Barcın and Revenia, than to each other (figure 2b; electronic supplementary material, table S12). Likewise, all Anatolian groups (Boncuklu, Tepecik-Çiftlik and Barcın) were genetically closer to Revenia than they were to each other (figure 2b; electronic supplementary material, table S13).

    • (vi) All European early farmer populations examined were genetically closer to Revenia than to each other (electronic supplementary material, figure S4c, table S14)

    Observations (ii), (iii) and in particular (v) are intriguing. We asked whether these could be consistent with a number of demographic scenarios, assuming a phylogenetic topology that included Iran, WHG and Aegean/Anatolian populations, estimated by [45]. We considered the following scenarios: (a) separate extreme bottlenecks in the ancestors of the two central Anatolian populations (possibly causing differentiation between the central Anatolian populations from each other, and from all other groups, (b) independent gene flow events from external sources (WHG and Iran) into the two central Anatolian groups (possibly causing differentiation between the two), (c) independent gene flow from WHG and Iran into the Aegean, and (d) independent gene flow from WHG, Iran and the two central Anatolian lineages into the Aegean. We performed serial coalescent simulations using realistic settings and compared the results with the observed D-statistics. We could only replicate the observed results under scenario (d) that describes rampant admixture in the Aegean (electronic supplementary material, figure S5).

    4. Discussion

    The analyses presented here highlight two points regarding the process of Neolithization. First, the observation that the two central Anatolian populations cluster together to the exclusion of Neolithic populations of south Levant or of Iran restates the conclusion that farming in central Anatolia in the PPN was established by local groups instead of immigrants, which is consistent with the described cultural continuity between central Anatolian Epipalaeolithic and Aceramic communities [9,64]. This reiterates the earlier conclusion [32] that the early Neolithization in the primary zone was largely a process of cultural interaction instead of gene flow.

    The second point relates to whether Aegean Neolithization (post-7000 cal BCE) involved similar acculturation processes, or was driven by migration similar to Neolithization in mainland Europe—a long-standing debate in archaeology [16,20,22,27,28]. Here, we discuss the two scenarios based on the genetic analysis.

    (a) Model 1: migration from Anatolia to the Aegean

    A recent study reported that by the seventh millennium BCE the eastward border of the WHG gene pool extended to the Iron Gates (on the border between Romania and Serbia) [33]. Plausibly, during the Early Holocene, the WHG population could also have been present along the Aegean coastline, such that the border between central Anatolian and WHG gene pools ran along west Anatolia. If so, the Aegean Neolithization must have involved replacement of a local, WHG-related Mesolithic population by incoming easterners.

    If migration occurred, where did it originate? Because Revenia and Barcın cluster with PPN and PN central Anatolian Neolithic groups to the exclusion of the south Levant (figure 1c; electronic supplementary material, figure S3c), the latter is unlikely to be the source, leaving central Anatolia or south Anatolia (north Levant) as potential origins.

    Notably, the north Aegeans (Revenia and Barcın) show higher diversity than the central Anatolians. We had earlier shown that the highest-quality Barcın genome carries a smaller proportion of short runs of homozygosity than the highest quality Boncuklu genome [30], which also supports the notion that the ancestral effective population size of the Aegeans was larger than those of central Anatolians. Moreover, we find that the north Aegeans share more alleles with eastern, western and southern gene pools, as estimated using the D-statistic (figure 2). Although the D-statistic can be sensitive to technical biases, our result is unlikely to be a technical artefact because (a) the north Aegean data were derived from two independent studies [31,37], (b) the Barcın data were produced using two different techniques, whole-genome shotgun sequencing and SNP capture, and (c) both Barcın and Revenia display the same population genetic patterns, suggesting that the admixture signals in the Aegean individuals are reproducible. In addition, although unknown population structure can complicate interpretation of the D-statistic [63], we note that the admixture estimates are consistent with the estimated higher genetic diversity in the Aegean.

    If the Revenia and Barcın individuals studied here were descendants of Anatolian Neolithic immigrants, they must have been recent settlers, as all samples analysed here date to early stages of the Aegean Neolithic (Revenia: 6438–6264 and Barcın: 6500–6200 cal BCE). Furthermore, if the migration was directly of central Anatolian origin (represented by Boncuklu and Tepecik-Çiftlik), the putative migrants must have admixed with populations carrying alleles of distinct gene pools (the Levant, Caucasus/Iran and WHG) within a few centuries, in order to explain our observations above (figure 2a).

    Alternatively, the migration event could have originated from the Anatolian south coast or north Levant [27] (currently no genome data are available from these groups). This region could have hosted a hypothetical central Anatolian-related population exposed to admixture from CHG-, Iran- and the Levant-related gene pools in earlier millenia. A south Anatolian population could have been in contact with different central Anatolian populations from the Konya Plain (Boncuklu) and Cappadocia (Tepecik-Çiftlik), explaining the affinity of both Boncuklu and Tepecik-Çiftlik to Barcın. A seafaring population could also be in genetic contact with putative WHG-related populations of the Aegean. This hypothetical population could have initiated the Cyprus Neolithic in the eleventh millennium BCE and later Aegean Neolithic communities in the seventh millennium BCE [27].

    One surprising observation here is the apparent absence of WHG-like ancestry in Late Neolithic/Chalcolithic Aegean genomes: admixture analysis results from two individuals from northwest Anatolia (Kumtepe, approx. 5000 BCE) [65] and four individuals from south Greece (Franchti Cave and Diros, approxi. 4000 BCE) [33] all lack noticeable WHG-like ancestry components [30,31,33]. This contrasts with WHG admixture emerging in European farmer populations in the Middle and Late Neolithic [36,38], and perhaps earlier in the Balkans [33], indicating the persistence of Mesolithic populations in Europe after Neolithic migrations. Therefore, if the Mesolithic populations of the Aegean coast had indeed been WHG-related, they must have been fully replaced by the eastern migrant farmers.

    (b) Model 2: adoption of Neolithic elements by local foragers

    Alternatively, the Aegean coast Mesolithic populations may have been part of the Anatolian-related gene pool that occupied the Aegean seaboard during the Early Holocene. Under this scenario, the north Aegean PN populations would be at least partial descendants of local hunter–gatherers who adopted Neolithic lifestyle post-7000 cal BCE, triggered by contacts with central Anatolian and the Levantine populations. The following events would be conceivable. (a) during the Last Glacial Maximum (LGM), the Aegean evolved into a refuge hosting a significant human population, which is in line with climatic modelling [6669]; estimates of human population density during the Marine Isotope Stage 2 in west (but not central) Anatolia reach one of their highest levels in Europe [70,71]. The existence of an Aegean human population going back to the LGM is also consistent with mitochondrial haplogroup-based analyses [72], and that Anatolian-like mitochondrial haplogroups are found also in Mesolithic Balkan and Aegean populations [31,33]. (b) Following the LGM, Aegean emigrants dispersed into central Anatolia and established populations that eventually gave rise to the local Epipalaeolithic and later Neolithic communities, in line with the earliest direct evidence for human presence in central Anatolia ca 14 000 cal BCE [60]. This hypothetical out-of-the-Aegean event coincides with the post-LGM Near East-related migration signatures in European Mesolithic genomes [73]. (c) Between the LGM and post-7000 cal BCE Neolithization, WHG, Natufian and Caucasus/Iran-related groups admixed with north Aegeans, differentiating the latter from their central Anatolian relatives and leading to our observations in figure 2a. (d) Post-7000 cal BCE, there occurred additional, albeit limited central Anatolian gene flow back into the Aegean, giving rise to our observation in figure 2b.

    (c) The archaeological evidence

    Both the migration and acculturation models for Aegean Neolithization enjoy support from material culture investigations, but the overall evidence points to a complex process where Aegean societies were culturally influenced by diverse sources, including the central Anatolian Neolithic, the Levant Neolithic and possibly local Mesolithic traditions. In contrast to the relative homogeneity of European Neolithic cultures, such as the LBK and Cardial, the Aegean Neolithic is noted for its diversity [64]. Variation in Neolithic Package elements and primary zone traditions is notable across Aegean sites, among regions (e.g. east and west of Marmara), even between closely neighbouring villages [16,17,20,64,7479]. This diversity includes, for example, obsidian, with Greek Aegean (Melos) [80,81] or mainland Anatolian (Cappadocian) [82] sources being preferred in some settlements, and yet other settlements showing no evidence of obsidian use [64]. Cultural trait diversity involves architecture, tool types, ceramics and symbolic elements (such as figurines and intramural burial), which may show partial similarities to either central Anatolia or to the Levant, or may be unique [16,64,75]. For instance, intramural burial, a common feature among primary zone sites, is also widespread in east Marmara early Neolithic villages (including Barcın), but totally absent in settlements only 200 km west [16]. Mesolithic-like lithic industries and the prominence of seafood in some settlements further imply the continuing presence of Aegean Mesolithic traditions into the Neolithic [16,20,27,83,84]. Indeed, lively seafaring activity was prevalent in the Mediterranean and the Aegean already by the eleventh millennium BCE [85,86], as evidence from Cyprus, Crete, Franchti, Cyclops Cave, Ouirakos and other Aegean island and coastal mainland Mesolithic sites demonstrates [20,27,83,8694].

    Instead of a single-sourced colonization process, the Aegean Neolithization may thus have flourished upon already existing coastal and interior interaction networks connecting Aegean foragers with the Levantine and central Anatolian PPN populations, and involved multiple cultural interaction events from its early steps onward [16,20,64,74]. This wide diversity of cultural sources and the potential role of local populations in Neolithic development may set apart Aegean Neolithization from that in mainland Europe. While Mesolithic Aegean genetic data are awaited to fully resolve this issue, researchers should be aware of the possibility that the initial emergence of the Neolithic elements in the Aegean, at least in the north Aegean, involved cultural and demographic dynamics different than those in European Neolithization.

    Data accessibility

    This article has no additional data.

    Authors' contributions

    M.S., I.T., G.M.K., A.G. and Ç.A. designed and supervised the study; G.M.K., D.K. and A.P.S. analysed data with contributions from H.M.D., R.Y., E.A., P.F.; Ç.A., D.B., A.M.B. and C.C.B. provided archaeological and ecological interpretations; M.S., I.T., G.M.K. and Ç.A. wrote the manuscript with input from all authors.

    Competing interests

    We declare we have no competing interests.

    Funding

    This work was supported by TÜBİTAK (grant no. 114Z927 to M.S.), TÜBA GEBİP (to M.S.), METU BAP (to M.S.), Knut and Alice Wallenberg Foundation: ‘1000 ancient genomes' (to A.G.). E.A. was supported by the Institutional Development Program of the University of Ostrava and by EU structural funding Operational Programme Research and Development for Innovation, project No. CZ.1.05/2.1.00/19.0388.

    Acknowledgements

    We thank T. Günther, A. Munters, S. C. Acan, A. Birand, Y. S. Erdal, A. Fairbairn, A. Omrak, M. Jakobsson, M. Krzewinska, F. Özer, R. Özbal, L. Excoffier and J. Stora for technical support in analysis, suggestions and/or comments on the manuscript, and two anonymous reviewers for helpful suggestions.

    Footnotes

    Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3929803.

    Published by the Royal Society. All rights reserved.

    References