A genomic snapshot of demographic and cultural dynamism in Upper Mesopotamia during the Neolithic Transition

Upper Mesopotamia played a key role in the Neolithic Transition in Southwest Asia through marked innovations in symbolism, technology, and diet. We present 13 ancient genomes (c. 8500 to 7500 cal BCE) from Pre-Pottery Neolithic Çayönü in the Tigris basin together with bioarchaeological and material culture data. Our findings reveal that Çayönü was a genetically diverse population, carrying mixed ancestry from western and eastern Fertile Crescent, and that the community received immigrants. Our results further suggest that the community was organized along biological family lines. We document bodily interventions such as head shaping and cauterization among the individuals examined, reflecting Çayönü’s cultural ingenuity. Last, we identify Upper Mesopotamia as the likely source of eastern gene flow into Neolithic Anatolia, in line with material culture evidence. We hypothesize that Upper Mesopotamia’s cultural dynamism during the Neolithic Transition was the product not only of its fertile lands but also of its interregional demographic connections.

Çayönü retains a long and uninterrupted sequence from the beginning of the PPN to the emergence of pottery during the latter part of the Neolithic. Even though the archaeological remains exhibit significant dynamism and external interaction within this lengthy occupation, there is no archaeological evidence for complete population turnover(s) that might have resulted in these changes. Instead, the archaeological record shows a continuous development through time, likely local responses to both local and non-local contingencies. In this regard, it is possible to surmise that the genetic composition of the PPNA occupants of Çayönü was not radically different from the main profile identified among the sampled PPNB individuals from the site. In light of the similarities in material culture and the spatial proximity, such a genetic profile might also apply to the broader Upper Tigris region in earlier time periods. Yet it becomes harder to evoke such arguments as one moves further east following the Tigris and its tributaries and further west towards the Euphrates. Even though parallel archaeological patterns can be detected in Upper/Middle Euphrates and Northwest Zagros, communities in these regions were also susceptible to further genetic interchanges with the Levant, Anatolia, and Central Zagros, most readily visible in the networks of obsidian transmission that they were entangled in (30,137).

Supplementary Note 3
Among the four Fertile Crescent populations studied (Levant, C Anatolia, U Mesopotamia, Zagros), the differentiation between Zagros and the other groups was found to be higher than what would be expected simply from isolation-by-distance, which we interpret as "lowest effective migration" between Zagros and Upper Mesopotamia. This east-west differentiation is estimated using the outgroup-f3 statistic calculated between pairs of individuals from different groups in the form of f3(outgroup; popA, popB). This is theoretically not affected by strong drift (bottlenecks) within groups, but reflects average split time between groups (26). Our f3 analysis thus suggests earlier split times between Zagros populations and other "Fertile Crescent" groups, relative to split times among the latter populations (i.e. Anatolia, Levant, as well as Çayönü). Moreover, the difference cannot be fully explained by geographic distance between Zagros and the other populations. We refer to this pattern as east-west differentiation. Importantly, the observed differentiation pattern does not require population isolation during the Neolithic. In fact, the Çayönü genomes indicate the contrary, such that in Çayönü we find both east and west Fertile Crescent ancestry. Instead, the east-west differentiation is likely the remnant of earlier population isolation, e.g., isolation during the LGM. Although post-LGM population admixture may have partly erased this differentiation signature, it can still be observed in our early Holocene genomic data.

Supplementary Note 4
The genetic diversity differences between two populations, say higher diversity in population A than population B, can arise through different means: (1) population A may carry higher interindividual variability in its ancestry proportions than B (e.g. as ref. (138) has described for ancient Rome); (2) within-population ancestry proportions may be equally homogenous between A and B, and (2a) population A may be composed of admixture between genetically more distinct groups, and thus each individual will carry higher heterozygosity, or (2b) population A may have higher effective population size, and thus higher heterozygosity than B. The qpAdm analyses suggest that among these possible mechanisms, not (1) but (2a) and (2b) may contribute to higher diversity estimates in Çayönü compared to Boncuklu or Aşıklı (Fig. 5A).

Supplementary Note 5
We use "Anatolia" here following the traditional geographic definition, referring to the west of the Anatolian Diagonal, which is defined as the line running from southwest to northeast Turkey.   Principal Component Analysis (PCA) of Southwest Asian early Holocene populations. The ancient genomes (coloured squares and circles, with Çayönü genomes shown with red stars) were projected onto principal components calculated using 55 present-day West Eurasian populations (grey dots) (Table S4). We computed 95% confidence intervals using the ellipse function implemented in smartpca (90) for all Çayönü individuals and highlighted ellipses of two individuals, cay008 and cay015, also labelled in the figure. We note three points: (a) a Southern Levant Neolithic individual (KFH2_KFH002) falls into the Anatolian cluster, which is consistent with earlier observations (7). (b) The cay015 genome appears in the Southern Levant cluster. However, it has the lowest coverage in our sample (only 6,986 SNPs available to this analysis) and its 95% confidence range overlaps with most of Çayönü individuals. (c) The cay008 genome stands out as an outlier, such that its range overlaps with nearly none of the individuals belonging to the main Çayönü cluster.

Fig. S4.
Modelling Çayönü with Zagros or Caucasus sources. The figure shows the results of feasible qpAdm models of the Çayönü core group (9 individuals) and cay008 genomes. Each line represents a model. These include either Zagros_N (Central Zagros Neolithic) or CHG (South Caucasus pre-Neolithic) genomes as source, in addition to Anatolia EP (Epipaleolithic) or Anatolia PPN, and South Levant Neolithic. All models are feasible with p>0.01. The full results are presented in Table  S6.

Fig. S5. ADMIXTURE analysis of modern-day West Eurasian (A) and ancient Southwest Asian (B)
genomes. We note that we excluded one of the individuals from each close kin pair in ADMIXTURE, to avoid creating a bias in allele frequencies (e.g. (91)). K=3 yields the lowest cross validation (CV) error.

Fig. S6.
ADMIXTURE analysis of Southwest Asian ancient genomes at K=2,3,4,5. The x-axis was ordered temporally from earliest to latest.

Fig. S7.
Modelling ancestry components of each individual Çayönü genome. Each bar represents a qpAdm model result, with grey bars indicating that the relevant models did not fit the data (p < 0.01). In panels (A) and (B), the Anatolian Epipaleolithic genome from Pınarbaşı (Anatolia EP) and in panels (C) and (D), Pre-Pottery Neolithic populations (Aşıklı and Boncuklu) from Central Anatolia were used as Anatolian sources. We modelled each individual as two-or three-way models shown in A, C and B, D, respectively. While in two-way models we used Anatolian and Central Zagros sources, three-way models additionally include South Levant Neolithic sources. Genomes are ordered on the y-axis from earliest (bottom) to latest (top). Individuals with >0.05X genome coverage are marked with asterisks. Notably, models with Levant-related components appear plausible only for these individuals with these relatively higher quality genomes in our sample. We did not include models of the cay008 genome, the results which are presented in Figure 2C. West Anatolia. Small dots represent genetic distance between each pair, whereas larger dots show median values of each pair set. We observe a significant difference between genetic distances between co-burial pairs and pairs found in different buildings in Çayönü (permutation test p<0.001), and a similar but non-significant trend in Aşıklı (p>0.05). (B) Frequencies of co-buried individuals with or without identified close genetic kin.   and (B) maternal relatives of cay008 to resolve her relationship with cay013 who is an adult female. Blue circles with dots correspond to the possible female third-degree relatives of cay008. Given that cay008 is an infant all relatives were simulated as related through the mother or father but not as descendants. The shaded area shows potential values similar to real data where X-chromosomal θ is higher than autosomal θ. Red dots represent θ values 1 -Normalised P0 between cay008 and cay013 individuals, whereas horizontal and vertical bars infer 95% confidence intervals. Triangles show the potential female relatives falling in this range. Annotated numbers correspond to relatives in the two pedigrees shown in Supplementary Figure 10. Jitter was added to the points to visualise all the data.

Fig. S13.
Phylogenetic tree of Y-chromosome computed with PathPhynder (80). All male individuals were placed onto the tree using the best path method. Hematoma on the left femoral diaphysis. Tables   Table S2. Material culture elements from PPN occupations at given sites. Table S3. Sequencing statistics of Çayönü individuals. Green rows represent merged libraries of each individual. Table S4. Published ancient and modern genomes used in this study. Table S5. D-statistics showing population affinities in Southwest Asia. Table S6. Two-and three-way qpAdm models. Table S7. Permutation test results of with-in population diversity in Southwest Asia. Upper triangle shows p-values and bottom triangle shows the effect size. Table S8. Kinship analyses on autosomes performed using READ and NGSRelate v2. Table S9. Kinship analyses on X-chromosome performed using READ and NGSRelate v2.