In this work we address both the question of maternal kinship relations of individuals from Çatalhöyük, and the genetic affinities of Central Anatolian populations and what follows their potential relation to the westward spread of the Neolithic. We present ten new complete mitochondrial genomes from Çatalhöyük individuals buried under the floors of adjacent buildings dated to classic levels (Mellaart Phase VI A) of its occupation.

Çatalhöyük was undoubtedly a part of a large, far-reaching exchange network [ 23 24 ] and could have potentially participated in the exchange of both goods and ideas. Elements of Çatalhöyük origin began to emerge in particular in North-western Anatolia in the middle of the 7th millennium BCE [ 25 ]. This process was unquestionably complex, presumably involving several subsequent impulses [ 15 ], as it took two millennia for the Neolithic to spread first to western Anatolia, with the earliest dates for the Neolithic being around the late 8th millennium BC [ 14 ], and then towards Marmara Region in the late 7th millennium BC [ 26 ]. The questions concerning to what extent those impulses towards west and North-western Anatolia were connected with gene flow and what was the role of the autochthonous hunter-gatherers in the spread are yet to be resolved. In some parts of the region intermediate and mixed traditions and economies have been observed [ 27 ]. At the same time, genomic data, both from the Marmara Region [ 28 29 ] and Central Anatolia [ 30 ], shows the genetic similarity of those regions and their close genetic affinity with Central European Neolithic populations. Those results support the leading role of the terrestrial route of the Neolithic spread both within and outside of the Anatolia.

The emergence and expansion of the Neolithic within and outside of Anatolia is another issue that could be addressed with ancient DNA data from Çatalhöyük. This process is thought to be a sum of several waves and trajectories of migration [ 13 15 ]. The Neolithic in Central and South-western Anatolia is thought to have developed under the influences from the upper Euphrates valley in the span of a thousand years, with the earliest evidence of some degree of Neolithic lifestyle seen in Central Anatolia in the second half of the 9th millennium BC [ 16 ]. However, the emergence of Neolithic societies in Central Anatolia was also proposed to be an autonomous process, involving local hunter-gatherers adopting the Neolithic lifestyle under the influence of farming communities from South-eastern Anatolia and Levant [ 16 18 ]. Furthermore, it has additionally been proposed that the region of Central Anatolia might not have contributed significantly to the subsequent westward movement of Neolithic tradition, as both archeological [ 19 20 ] and zooarchaeological [ 21 ] data suggest that it constituted a distinctive cultural zone. At the same time, the maritime colonization originating in the Levantine coast has been proposed as the major factor contributing to the development and the spread of the Neolithic in the Aegean coast of western Anatolia [ 22 ]. It is thought that this process did not involve any local populations, as Mesolithic occupation was sparse in the parts of the region where the Neolithic first appeared [ 14 ].

Initially it was proposed that Çatalhöyük individuals buried together in the same building were biologically related, and groupings of houses and constituting neighborhoods were defined by biological affinity [ 5 ]. Then, an uneven distribution of burials among different houses was interpreted as evidence for some of them being a burial place for larger household communities, composed of a number of houses inhabited by nuclear families [ 6 ]. A recent study based on dental phenotypes of individuals found in Çatalhöyük burials showed that individuals with close biological affinity spanned across several buildings [ 7 ]. This result was interpreted as evidence of the lack of kinship patterning in burials found at the site. However, the correlation between biological distances based on both morphological traits and genetic kinship is poorly understood, as both types of data are rarely available for the same set of samples. In the few studies where direct comparison between morphological and genetic data was available, the results were inconsistent, pointing towards weak correlation [ 8 ]. Several approaches and tools for genetic kinship estimation based on ancient DNA have been recently published [ 9 11 ], and although these tools were developed with low coverage data in mind, they still depend on significant overlap in nuclear single nucleotide polymorphisms (SNPs) between analyzed samples. However, where overall ancient DNA (aDNA) preservation between samples is poor and/or deeper sequencing data is not feasible, mitochondrial (mt) genomes can be used to exclude maternal kinship [ 12 ].

Neolithic Çatalhöyük (7100–5950 BC) is a world-renowned Neolithic settlement. Its size, remarkable preservation, presence of numerous works of Neolithic art, and large amounts of archeological data obtained through meticulous excavation have consolidated its unquestioned importance in the identification of a wide range of constituent elements of the Neolithic [ 1 ]. The settlement was composed of a conglomeration of clustered neighborhoods with clearly defined modular house units [ 2 ]. All houses were apparently occupied and used for domestic purposes [ 3 ]. Burials were located under the floors of most buildings, especially under elevated platforms in northern and eastern parts of the living rooms. However, some of the buildings, notably the ones with more elaborate art installations, contained more burials (up to almost 70 individuals, more than one would expect from the estimated number of their inhabitants), implying their special status [ 4 ]. Those buildings are thought to have been “history houses” that provided or controlled ancestors and rituals for a larger kin or other group [ 3 ].

2. Materials and Methods

Four adjacent, roughly contemporary buildings from Çatalhöyük South Area, dated to Mellaart Phase VI A (6450–6380 cal. BC [ 31 ]), were selected for the study ( Figure 1 C). We assumed that the selected buildings represented ordinary houses as neither of them was recognized as a “history house” by the researchers of the site, however, an above-average number of art installations were found in building 80, and all of the buildings, with the exception of building 89, contained more than 10 burials. All available individuals excavated from those buildings were sampled. Where possible, petrous part of temporal bones were collected. The samples were all taken from the Çatalhöyük Research Project depot with the use of disposable gloves and facemasks. In total, 47 bone samples were acquired from 37 skeletons, including ten individuals from building 96, six from building 97, five from building 89, and 16 from building 80. The detailed information on the selected samples can be found in the Supplementary Information Text and Supplementary dataset S1

DNA was extracted from teeth and petrous parts of temporal bones in laboratories dedicated to working with human aDNA. The surface of the samples was decontaminated with the use of 2% bleach and UV light, and only the inner part of both the teeth and petrous bones were drilled for extraction. DNA isolation was performed both at the Middle East Technical University in Ankara (METU), Turkey, and at the Faculty of Biology, Adam Mickiewicz University in Poznan (AMU), Poland. In the METU laboratory, the DNA was extracted using a silica-powder-based method and modified binding buffer, as described by Allentoft et al. [ 32 ]. In the AMU facility, the DNA was extracted using a silica-based method as in [ 33 ], with modification by [ 34 ]. Total genomic DNA libraries for all samples were constructed at the AMU laboratory using the methods described previously by Juras et al. [ 12 ].

The genomic libraries underwent Illumina sequencing (150 bp PE, ca. 1.4 mln reds/sample) at the Science for Life Laboratory (SciLifeLab) facility in Stockholm (NGI Stockholm), Sweden. All preliminary pipeline computations of the sequencing data were undertaken on resources provided by the Swedish National Infrastructure for Computing (SNIC) through the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) [ 35 ].

The RNA baits for capture enrichment of complete mitochondrial genomes were prepared following the protocols described in Juras et al. [ 36 ]. Two rounds of mitochondrial DNA (mtDNA) enrichment were carried out on 22 libraries that showed either a proportion of reads mapping to human reference genome (version hs37d5) equal to at least 0.4%, or mtDNA coverage ranging from 0.02 × to 5 × after initial Illumina shotgun screening ( Supplementary dataset S1 ). The libraries enriched in mtDNA were then sequenced on Illumina HiSeq2500 (125 bp, paired end, each library 1/80 lane) in the SciLifeLab facility in Stockholm (NGI Stockholm), Sweden.

DNA sequencing data from both shotgun screening and mtDNA capture were processed with the use of a customizable analytical pipeline, described in [ 37 ]. The adapters were removed and read pairs were merged, requiring an overlap of 11 bp and summing up base qualities using MergeReads-FastQ_cc.py script, according to Meyer & Kircher [ 38 ]. BWA software package version 0.7.8 [ 39 ] with the parameters -n 0.01 -o 2 and disabled seeding, was used to map merged reads as single-end reads against both the revised Cambridge Reference Sequence (rCRS) [ 40 41 ] (GenBank: NC_012920) and human reference genome (version hs37d5). To collapse duplicate sequence reads with identical start and end coordinates, FilterUniqueSAMCons.py was used [ 38 ]. The ratio of reads mapping to Y and X chromosomes (Ry) (with mapping quality greater than 30) was calculated to assign molecular sex of individuals [ 42 ].

The mtDNA capture binary alignment map (bam) files were merged with shotgun screening data, using the merge tool from SAMtools package v1.8 [ 43 ]. Misincorporation patterns for merged files were assessed using mapDamage v2.0.5 with the default parameters [ 44 ]. Contamination estimates for mtDNA sequences were then preformed using both contamMix_1.0-10 script [ 45 ] and Schmutzi package v1.5.4 [ 46 ] with the default parameters. Any sample that failed at least one of those tests, showing more than 18% contaminating sequences, was discarded from further analysis. Consensus sequences were built using ANGSD v0.910 [ 47 ]. Only reads with a mapping score of 30, a minimum base quality of 20, and positions with a minimum coverage of 3 were accepted [ 48 ]. All the computations were performed using resources provided by The Polish Grid Infrastructure (Pl-Grid). Mitochondrial haplogroups (hgs) were assigned for each individual utilizing HAPLOFIND [ 49 ], and Haplogrep [ 50 ] both based on the PhyloTree phylogenetic tree build 17 [ 51 ]. Biomatters IGV software v2.3.66 [ 52 ] was used to visualize final sequences, as well as mutations reported as unexpected or missing in the original binary alignment map (BAM).

ST ) are described in detail in For comparative analyses, complete ancient mtDNA genomes were obtained from the literature and the ancient human mitochondrial genomes database (AmtDB) [ 53 ]. Where available reconstructed fasta files were acquired, and in cases where only whole genome data was available, the mitochondrial genomes were reconstructed from the available bam files using the pipeline described above. All samples were then grouped into sets of at least 10 individuals based on their dating, geographical location, and/or attribution to particular archeological cultures. Only one sample from the confirmed kin pairs and groups was selected for population analyses. Additionally, READ [ 10 ] was used on 67 reference Neolithic Anatolian and Near Eastern samples in order to detect potential genetic kinship relations missed in previous studies. All comparative populations and samples used for principal component analysis (PCA), t-Distributed Stochastic Neighbor Embedding (t-SNE), and pairwise genetic distances (F) are described in detail in Supplementary dataset S2

For the purpose of this study, due to the limited number of samples available for Iran and Turkmenistan Neolithic and Chalcolithic, the samples were grouped together into Neolithic Middle East group (NME). Similarly, Pre-Pottery Neolithic samples from Jordan were grouped together with epipaleolithic Natufian samples from the same region into the Natufian and Neolithic Levant group (NNL). Furthermore, the Bronze Age samples from Turkey were grouped together as Bronze Age Near East group (BNE) with Jordan samples from the same period, and Neolithic, Chalcolithic, and Early Bronze Age samples from Armenia were merged into the Neolithic to Bronze Age Caucasus group (NBC).

55, The Çatalhöyük samples were analyzed as part of the Central Anatolia Neolithic group containing additionally three individuals from Boncuklu Höyük and three from Tepecik-Çiftlik sites. It has been shown that the European gene pool was shaped by three major ancestral populations, including autochthonous hunter-gatherers, and two migrant groups from Near East in the Early Neolithic and Eurasian steppe in the Late Neolithic period [ 54 56 ]. However, since this work is mostly focused on genetic relationships within the early farming populations, all the analyses were performed using both: (i) only the Neolithic to Bronze Age populations from Near and Middle East with the addition of initial farming populations from Europe (as seen on ( Figure 1 A)), and (ii) all the above with the addition of the Yamnaya steppe groups and hunter-gatherers from Europe (divided into Western, Eastern, and Balkan populations), added as proxies of the other two major components of the European gene pool. Only one individual in each pair of known first degree relatives was used in the analyses. Additionally, individuals for which less than 85% of the mitochondrial genome was recovered were excluded from the analyses.

The map in Figure 1 was generated using QGIS 2.12.2. PCA of mtDNA hgs frequencies was calculated using Python 3.5 and the Scikit-learn v. 0.18.1 [ 57 ] package. The PCA results and mtDNA hgs loadings were plotted with the use of Matplotlib 1.5.1 Python package [ 58 ].