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

1. Introduction

The primary zone of Neolithization in western Eurasia encompassed the Levant, Taurus-Zagros ranges of Mesopotamia, central Anatolia and Cyprus [1–4]. 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 [7–9]. Between ca 8500 and 7000 cal BCE, community sizes increased, architectural elaboration intensified and a subsistence economy based on agriculture gradually became the norm [10–14]. 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’ [15–18].

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 [19–21]. 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,22–27]. 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’ [30–33], 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 [31–33] (electronic supplementary material, figure S1). In subsequent millenia, however, WHG-like ancestry appears in Middle and Late Neolithic European populations [36–38]. 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 [39–42].

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 [30–32,35,37,38,42–46]. 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 f₃-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 f₃-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 f₃-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 f₃-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 f₃-statistics, as genetic distance measure between a pair of populations, X and Y, we used: 1-f₃(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⁻⁹ bp yr⁻¹, and a recombination rate of 1.00 × 10⁻⁸ bp yr⁻¹, and assumed 25 years per generation, again following [45]. We set the effective population size (N_e) 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.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

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 [66–69]; 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,74–79]. 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,86–94].

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.

Footnotes

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