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

Abstract

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.

Fig. 1.. Spatiotemporal distribution of the samples and the population structure of Neolithic Southwest Asia.

(A) Timeline of ancient Southwest Asian individuals used in the analyses. Colored horizontal bars at the bottom represent the subperiods of the Neolithic Era in Southwest Asia. (B) The map shows EP and Neolithic populations from Southwest Asia. Shaded areas mark PPN period cultural zones (referred to as the Aceramic period in C Anatolia). (C) Çayönü building types and their approximate dates of use, considered as evidence for Çayönü’s cultural openness and ingenuity. Modified from (112). (D) The first two dimensions of the MDS plot of genetic distances. The MDS summarizes the genetic distance matrix among ancient genomes calculated as (1 − outgroup f₃) values. Outgroup f₃-statistics were calculated as f₃(Yoruba; individual₁, individual₂). The labels represent the following sites: Anatolia EP: Pınarbaşı; Anatolia PPN: Boncuklu and Aşıklı Höyük; Anatolia PN: Çatalhöyük and Barcın Höyük; Levant EP: Natufian; Levant PPN: Ain’ Ghazal, Kfar HaHoresh, Motza, and Ba’ja; C Zagros N (Central Zagros Neolithic): Ganj Dareh, Tepe Abdul, and Wezmeh Cave; S Caucasus EP (South Caucasus EP): Kotias and Satsurblia. See note S5 for a definition of “Anatolia.” PPNA, Pre-pottery Neolithic A; PPNB, Pre-pottery Neolithic B; PPNC, Pre-pottery Neolithic C.

Fig. 2.. Genetic affinities of the Çayönü population with the neighboring populations.

Formal tests computed in the form of (A) D(Yoruba, Çayönü/cay008; pop2, test) and (B) D(Yoruba, pop1; Çayönü/cay008, Anatolia EP/PPN/PN). Z scores were corrected with the Benjamini-Hochberg multiple testing correction (86). Horizontal bars represent ±2 SE. (C) qpAdm modeling of the Çayönü group and cay008. The local Çayönü group or an outlier cay008 individual was the “target”; Central Anatolia EP, Central Zagros Neolithic, and South Levant Neolithic samples were sources for both targets. The local Çayönü group was also used as the “source” for modeling of cay008. Model 1 represents the model with “local Çayönü population + C Zagros N,” whereas model 2 includes “Anatolia EP + C Zagros N” as sources. Horizontal bars represent SEs of the coefficients. All three models yielded P > 0.05. We also cannot reject a three-way model of Central Anatolia PPN, Central Zagros, and South Levant at P > 0.01, while three-way models for cay008 were not infeasible (table S5). In all analyses shown in the figure, “Çayönü” represents the nine genomes listed in Table 1, excluding relatives and cay008.

Fig. 3.. Genetic isolation by distance in Southwest Asia.

(A) Correlation between geographic (x axis) and genetic (y axis) distance for Southwest Asia Neolithic populations. The red regression line shows the linear fit with 95% confidence interval. Each point represents pairs of individuals from Southwest Asia Neolithic. Pairs from the same site and pairs separated by >1000 years of time difference were not included. All regression lines were highly significant (P < 0.001). (B) The distribution of residuals that we calculated by subtracting the observed values from the predicted values obtained from the linear regression models in (A). In all analyses shown in the figure, Çayönü represents the nine genomes listed in Table 1, excluding relatives and cay008.

Fig. 4.. Testing temporal genetic change in Çayönü.

(A) Formal tests computed in the form of D(Yoruba, pop1; test, Boncuklu) where pop1 denotes CHG/S Levant N/C Zagros N, and test denotes radiocarbon-dated Çayönü individuals (Table 1), ordered from earliest to latest in bottom to top direction, respectively. Horizontal bars represent ±2 SE. (B) qpAdm model of each radiocarbon-dated Çayönü individual. The outlier individual cay008 was excluded. The y-axis is in the same order as in (A). Black outlines around the boxes show that the P value of the model is >0.05 (indicating feasibility), whereas gray outlines indicate 0.05 > P > 0.01.

Fig. 5.. Genetic diversity in Neolithic Southwest Asia.

(A) Small dots show pairwise genetic distance calculated as (1 − outgroup f₃) values for all pairs of individuals, whereas large dots show the median values of each population. The vertical gray line represents the total median across the eight populations. Deviation from the total median is shown with colored horizontal lines. Outgroup f₃-statistics were computed as f₃(Yoruba; ind₁, ind₂) where ind₁ and ind₂ represent individuals from the same archaeological site. (B) and (C) present ROHs in Southwest Asia. Sum of total ROH > 4 cM and number of total ROH > 4 cM are shown on the x and y axes, respectively. The baseline (red diagonal line) was computed using short ROH values (4 to 8 cM) in present-day West and Central Eurasian individuals to represent outbreed samples to determine the baseline. (C) is the zoomed version of (B) in which we draw the zoomed area with the gray rectangle. The red star denotes the only Çayönü individual, namely, cay007, that has more than 300,000 SNPs in the 1240K SNP Panel. The gray dots designate the ROH values for modern genomes.

Fig. 6.. Kinship coefficient (θ) estimates among Çayönü individuals.

Comparison of kinship coefficients inferred for 76 pairs using four different software [NGSRelate (44), READ (45), lcMLkin (46), and TKGWV2 (47)] is shown in (A). The figure only shows estimates when a pair had >2000 shared SNPs available in that analysis (which may differ for different software). (B) and (C) show autosomal and X-chromosomal estimates of θ, respectively, inferred from NGSRelate and READ. In the last two panels, NGSRelate θ estimates are shown on the x axes and READ θ estimates, calculated as (1 − normalized P₀), on the y axes. Vertical bars represent ±2 SE of P₀ values. Vertical dotted, dashed, and straight gray lines intersect with expected θ values for first-, second-, and third-degree relatives, respectively. Annotation with the gray label shows the pair cay008 and cay013. The full list of kinship coefficient estimates for all possible pairs is given in table S8.

Fig. 7.. Locations of Çayönü coburials interred in domestic buildings.

All three buildings belong to the cell building subphase. The figure shows plans of buildings coded (A) CA, (B) CL, and (C) CN. Red dots represent individuals analyzed in this study, pink dots represent individuals screened for aDNA but with insufficient preservation, and blue dots represent burials of other individuals within the same buildings. Black curved lines show the closely related pairs in each building and the estimated kinship level. The blue star in (B) indicates that the pair (cay015-cay033) is likely close genetic kin but were not reported as their number of shared SNPs was below the chosen threshold (<2000 SNPs) (Materials and Methods and table S8).

Fig. 8.. Cranial features of the cay008 toddler.

(A) Frontal flattening, post-coronal depression, bulging on the parietal tuber, and fronto-occipital grooving suggest a double-bandaged circular-type cranial deformation. (B) Cauterization with a circular depression found on the post-coronal area on the left parietal bone. The bone is very thin in the center, and the edge of the lesion is elevated. (C) An enlarged picture of a post-coronal depression and frontal flattening. (D and E) Endocranial lithic lesions similar to serpens endocrania symmetrica on the occipital bone. (F) Slightly developed cribra orbitalia on the right orbital roof. This lesion together with porotic hyperostosis is mainly related to anemia. (G) Cranial trepanation performed by drilling on the skull of Çayönü individual ÇT’78 KE 6-2/3a SK5 (not represented in our genetic sample).

Fig. 9.. Biplots of D-statistics illustrating excess allele sharing between Çayönü and post-7000 BCE populations from Central/Western Anatolia.

D-statistics were computed in the form of D(Yoruba, pop1; pop2, X), where X represents PN, Chalcolithic, and Bronze Age populations from the Anatolian Plateau. Each population is represented by a dot and error bars representing ±2 SE. The list of populations and D-statistics can be found in table S5. In both panels, pop1 corresponds to CHG on the x axes, whereas on the y axes, pop1 corresponds to the Çayönü population. pop2 is represented by Boncuklu (Central Anatolian PPN) in (A) and Pınarbaşı (Central Anatolian EP) in (B) in both axes. The slope of the diagonal dashed line is 1 showing x = y, and the intercept of both vertical and horizontal dotted lines is 0.