The genetics of an early Neolithic pastoralist from the Zagros, Iran.

The agricultural transition profoundly changed human societies. We sequenced and analysed the first genome (1.39x) of an early Neolithic woman from Ganj Dareh, in the Zagros Mountains of Iran, a site with early evidence for an economy based on goat herding, ca. 10,000 BP. We show that Western Iran was inhabited by a population genetically most similar to hunter-gatherers from the Caucasus, but distinct from the Neolithic Anatolian people who later brought food production into Europe. The inhabitants of Ganj Dareh made little direct genetic contribution to modern European populations, suggesting those of the Central Zagros were somewhat isolated from other populations of the Fertile Crescent. Runs of homozygosity are of a similar length to those from Neolithic farmers, and shorter than those of Caucasus and Western Hunter-Gatherers, suggesting that the inhabitants of Ganj Dareh did not undergo the large population bottleneck suffered by their northern neighbours. While some degree of cultural diffusion between Anatolia, Western Iran and other neighbouring regions is possible, the genetic dissimilarity between early Anatolian farmers and the inhabitants of Ganj Dareh supports a model in which Neolithic societies in these areas were distinct.

The agricultural transition started in a region comprising the Ancient Near East and Anatolia ~12,000 years ago with the first Pre-Pottery Neolithic villages and the first domestication of cereals and legumes12. Archaeological evidence suggests a complex scenario of multiple domestications in a number of areas3, coupled with examples of trade4. Ancient DNA (aDNA) has revealed that this cultural package was later brought into Europe by dispersing farmers from Anatolia (so called ‘demic’ diffusion, as opposed to non-demic cultural diffusion56) ~8,400 years ago. However a lack of aDNA from early Neolithic individuals from the Near East leaves a key question unanswered: was the agricultural transition developed by one major population group spanning the Near East, including Anatolia and the Central Zagros Mountains; or was the region inhabited by genetically diverse populations, as is suggested by the heterogeneous mode and timing of the appearance of early domesticates at different localities?

To answer this question, we sequenced the genome of an early Neolithic female from Ganj Dareh, GD13a, from the Central Zagros (Western Iran), dated to 10000-9700 cal BP7, a region located at the eastern edge of the Near East. Ganj Dareh is well known for providing the earliest evidence of herd management of goats beginning at 9,900 BP789. It is a classic mound site at an altitude of ~1400 m in the Gamas-Ab Valley of the High Zagros zone in Kermanshah Province, Western Iran. It was discovered in the 1960s during survey work and excavated over four seasons between 1967 and 1974. The mound, ~40 m in diameter, shows 7 to 8 m of early Neolithic cultural deposits. Five major levels were found, labelled A through E from top to bottom. Extended evidence showed a warren of rooms with evidence of under-floor inhumations within what may be burial chambers and/or disused houses10. The current Minimum Number of Individuals is 116, with 56 catalogued as skeletons that had four or more bones recovered11. The individual analysed here was part of burial 13, which contained three individuals, and was recovered in level C in 1971 from the floor of a brick-walled structure. The individual sampled, 13A (referred to as GD13a throughout the text), was a 30–50 year old female; the other individuals in the burial unit were a second adult (13B) and an adolescent (13).

The site has been directly dated to 9650–9950 cal BP7, and shows intense occupation over two to three centuries. The economy of the population was that of pastoralists with an emphasis on goat herding7. Archaeobotanical evidence is limited12 but the evidence present is for two-row barley with no evidence for wheat, rye or other domesticates. This implies that the overall economy was at a much earlier stage in the development of cereal agriculture than that found in the Levant, Anatolia and Northern Mesopotamian basin.

Results

The petrous bone of GD13a yielded sequencing libraries comprising 18.57% alignable human reads that were used to generate 1.39-fold genome coverage. The sequence data showed read lengths and nucleotide misincorporation patterns indicative of post-mortem damage, and had low contamination estimates (<1%, see Supplementary Fig. S3). The mitochondrion of GD13a (91.74X) was assigned to haplogroup X, most likely to the subhaplogroup X2, which has been associated with an early expansion from the Near East1314 and has been found in early Neolithic samples from Anatolia5, Hungary15 and Germany16.

We compared GD13a with a number of other ancient genomes and modern populations615161718192021222324252627, using principal component analysis (PCA)28, ADMIXTURE29 and outgroup f3 statistics30 (Fig. 1). GD13a did not cluster with any other early Neolithic individual from Eurasia in any of the analyses. ADMIXTURE and outgroup f3 statistics identified Caucasus Hunter-Gatherers of Western Georgia, just north of the Zagros mountains, as the group genetically most similar to GD13a (Fig. 1B,C), whilst PCA also revealed some affinity with modern Central South Asian populations such as Balochi, Makrani and Brahui (Fig. 1A and Fig. S4). Also genetically close to GD13a were ancient samples from Steppe populations (Yamanya & Afanasievo) that were part of one or more Bronze age migrations into Europe, as well as early Bronze age cultures in that continent (Corded Ware)1621, in line with previous relationships observed for the Caucasus Hunter-Gatherers24.

Figure 1

GD13a appears to be related to Caucasus Hunter Gatherers and to modern South Asian populations.

(A) PCA loaded on modern populations (represented by open symbols). Ancient individuals (solid symbols) are projected onto these axes. (B) Outgroup (, GD13a; Dinka), where Caucasus Hunter Gatherers (Kotias and Satsurblia) share the most drift with GD13a. Ancient samples have filled circles whereas modern populations are represented by empty symbols. (C) ADMIXTURE using K = 17, where GD13a appears very similar to Caucasus Hunter Gatherers, and to a lesser extent to modern south Asian populations.

We further investigated the relationship between GD13a and Caucasus Hunter-Gatherers using D-statistics3031 to test whether they formed a clade to the exclusion of other ancient and modern samples (Table S4). A large number of Western Eurasian samples (both modern and ancient) showed significant excess genetic affinity to the Caucasus Hunter-Gatherers, whilst none did with GD13a. Overall, these results point to GD13a having little direct genetic input into later European populations compared to its northern neighbours.

To better understand the history of the population to which GD13a belonged, we investigated the distribution of lengths of runs of homozygosity (ROH) (Fig. 2A). A bias towards a high frequency of both long and short ROH is indicative of past strong bottlenecks followed by population re-expansion. GD13a has a distribution with few long ROH (>2 Mb), similar to that of the descendants of Anatolian early farmers (represented by the European farmers NE115 and Stuttgart17). In contrast, both Western17 and Caucasus Hunter-Gatherers24 have relatively more long as well as short ROH. Thus, GD13a is the descendant of a group that had relatively stable demography and did not suffer the bottlenecks that affected more northern populations.

Figure 2

GD13a did not undergo a recent large population bottleneck.

(A) GD13a has similar runs of homozygosity (ROH) lengths to Neolithic individuals, while Caucasus Hunter Gatherers (Kotias and Satsurblia), like European Hunter Gatherers (Loschbour and Bichon), underwent recent large population bottlenecks potentially associated withthe LGM. (B) Map showing geographical location of Anatolian Neolithic samples, Caucasus Hunter Gatherers (CHG) and GD13a. Background colours indicate mean temperature (°C) of coldest quarter during the LGM (data from the worldclim database60 generated by the CCSM4 model)60, with LGM sea levels. Map of populations was generated with MATLAB R2015b (Mathworks, http://www.mathworks.com/)61.

The phenotypic attributes of GD13a are similar to the neighbouring Anatolian early farmers and Caucasus Hunter-Gatherers. Based on diagnostic SNPs, she had dark, black hair and brown eyes (see Supplementary). She lacked the derived variant (rs16891982) of the SLC45A2 gene associated with light skin pigmentation but likely had at least one copy of the derived SLC24A5 allele (rs1426654) associated with the same trait. The derived SLC24A5 variant has been found in both Neolithic farmer and Caucasus Hunter-Gatherer groups51524 suggesting that it was already at appreciable frequency before these populations diverged. Finally, she did not have the most common European variant of the LCT gene (rs4988235) associated with the ability to digest raw milk, consistent with the later emergence of this adaptation51521.

It is possible that farmers related to GD13a contributed to the eastern diffusion of agriculture from the Near East that reached Turkmenistan32 by the 6th millennium BP, and continued further east to the Indus Valley33. However, detecting such a contribution is complicated by a later influx from Steppe populations with Caucasus Hunter-Gatherer ancestry during the Bronze Age. We tested whether the Western Eurasian component found in Indian populations can be better attributed to either of these two sources, GD13a and Kotias (a Caucasus Hunter Gatherer), using D-statistics to detect gene flow into an ancestral Indian component (represented by the Onge). Overall, for all tests where a difference could be detected, Kotias and GD13a were equally likely sources (Fig. S9 and Table S6). Whilst the attribution of part of the Western Eurasia component seen in India to Bronze Age migrations is supported by dating of last contact based on patterns of Linkage Disequilibrium34, our analysis highlights the possibility that part of that component might derive from earlier contact during the eastern diffusion of agriculture.

Discussion

GD13a had little direct genetic input into later European populations compared to the Caucasus Hunter-Gatherers (its northern neighbours) as demonstrated using D-statistics. This lack of connectivity with neighbouring regions might have arisen early on, since we also find that hunter-gatherers from the Caucasus show higher affinity to Western Hunter-Gatherers and early Anatolian farmers; this result suggests the possibility of gene flow between the former and these two latter groups to the exclusion of GD13a. An alternative, but not mutually exclusive, explanation for this pattern is that GD13a might have received genetic input from a source equally distant from all other European populations, and thus basal to them.

The Last Glacial Maximum (LGM) made entire regions in northern Eurasia uninhabitable, and therefore a number of hunter-gatherer populations likely moved to the south. For Europe there may be a separation between Western and Eastern populations with minimal occupation of the Central European plains22. For Eastern Europe, Central Asia and the northern Near East, glaciation itself was less a factor. In these areas, our understanding of how populations weathered the LGM is still vague at best. It has previously been suggested that differences in the frequency of long and short runs of homozygosity in ancient samples may be associated with different demographic experiences during the LGM1524. Neolithic farmers, with their lower frequency of short ROH, have been argued to have been relatively little affected by the LGM compared to Western and Caucasus Hunter-Gatherers1524 which are characterised by more long ROH (>2 Mb). GD13a has a profile similar to that of the descendants of Anatolian farmers (i.e. early European farmers), suggesting that her ancestors also faced more benign conditions compared to populations further north. Superimposing the sampling locations of these groups onto climatic reconstructions from the LGM (Fig. 2B), however, does not reveal clear climatic differences among the regions. It is possible that the ancestors of the Anatolian and Ganj Dareh farmers spent the LGM in areas further south or east, which experienced milder climate. But it is also possible that they exploited local pockets of favourable climate (refugia). Whilst high elevation sites in the Zagros were abandoned during the LGM35, there are a number of sites in the valleys that were occupied during that period and might have experienced more favourable conditions36.

The archaeological record indicates an eastward Neolithic expansion from the eastern regions of the Near East into Central and South Asia3237. Our analysis shows that both the Caucasus Hunter Gatherer Kotias and GD13a are plausible sources for the Eurasian Ancestry found in that part of Asia. Even though part of the Western Eurasian component found in India can be linked to Bronze Age migrations by dating the last contact using Linkage Disequilibrium (thus coming from the Kotias lineage), our results highlight the possibility of an older contribution from a source genetically close to GD13a (which would be hard to disentangle from the later gene flow from the Steppe). Eventually, ancient DNA from the Indus Valley will be needed to detect conclusively whether any genetic traces were left by the eastward Neolithic expansion from the Near East, or whether this process was mostly cultural.

The presence of two distinct lineages (Anatolian-like agriculturalists and Zagros mountain herders) in the Near East at the beginning of the Neolithic transition raises an interesting question regarding the independence of innovations arising at different locations. Even within the Central Zagros, economies vary greatly in their rate and pathway towards Neolithisation35. Ganj Dareh, in the high Zagros, has the earliest known evidence for goat domestication789, and the foothills of the Zagros mountains have also been argued to have been the site of early farming3. In addition, early sites such as Sheikh-e Abad (11.650-9,600 cal BP) provide evidence of early stages of barley cultivation38. Were these innovations independent of similar achievements that made up the Neolithic package that North West Anatolians brought into Europe? Or were they exchanged culturally? If the latter, it would imply a cultural diffusion in the absence of much genetic interchange.

Methods

DNA extraction and library preparation

Sample preparation, DNA extraction and library preparation were carried out in dedicated ancient DNA facilities at University College Dublin. The dense part of the petrous bone was isolated, cleaned and sequenced following experimental procedures outlined in Gamba et al.15. DNA was extracted from 310 mg of ground bone powder using a double-digestion and silica column method as described in Gamba et al.39. Indexed Illumina sequencing libraries were constructed with a protocol based on Meyer et al.40 with modifications including blunt end repair using NEBNext End Repair Module (New England BioLabs Inc) and heat inactivation of Bst DNA polymerase15.

Sequence processing and alignment

Libraries were sequenced over a flow cell on a HiSeq 2000 at the TheragenEtex (South Korea) using 100 bp single-end sequencing. Adapter sequences were trimmed from the 3′ ends of sequences using cutadapt version 1.341, conservatively requiring an overlap of 1 base pair (bp) between the adapter and the read. Reads were aligned using BWA42, with the seed region disabled, to the GRCh37 build of the human genome with the mitochondrial sequence replaced by the revised Cambridge reference sequence (NCBI accession number NC_012920.1). Data from separate lanes were merged using Picard MergeSamFiles (http://picard.sourceforge.net/) and duplicate reads from the same library amplification were filtered using SAMtools rmdup43. Sequences were further filtered to remove those with mapping quality <30 and length <30 bp. Indels were realigned using RealignerTargetCreator and IndelRealigner from the Genome Analysis Toolkit44. The first and last 2 bp of each read were soft-clipped to a base quality of 2. The average genome-wide depth of coverage was calculated using the genomecov function of bedtools45. A summary of alignment statistics can be found in Table S1.

Authenticity of results

The data were assessed for the presence of typical signatures of post-mortem DNA damage4647. The sequence length distribution of molecules was examined as outlined in Gallego Llorente et al.48 (Fig. S2) while the prevalence of nucleotide misincorporation sites at the ends of reads was evaluated using mapDamage 2.0 and a random subsample of 1 million reads49 (Fig. S3).

The mitochondrial contamination rate was assessed by evaluating the proportion of non-consensus bases at haplogroup defining positions in the mitochondrial genome1550. Only bases with quality ≥20 were used. The X chromosome contamination rate could not be evaluated as the sample was determined to be female, using the script described in ref. 51.

Mitochondrial Haplogroup Determination

To determine to which haplogroup the mitochondrion of GD13a belonged, a consensus sequence was generated using ANGSD52. Called positions were required to have a depth of coverage ≥3 and only bases with quality ≥20 were considered. The resulting FASTA files were uploaded to HAPLOFIND53 for haplogroup determination. Coverage was calculated using GATK DepthOfCoverage44.

Dataset preparation for population genetic analyses

Genotypes were called in GD13a at sites which overlapped those in the Human Origins dataset (Lazaridis et al.17, filtered as described in Jones et al.24) using GATK Pileup44. Triallelic SNPs were discarded and bases were required to have quality ≥30. For positions with more than one base call, one allele was randomly chosen with a probability equal to the frequency of the base at that position. This allele was duplicated to form a homozygous diploid genotype for each position called in GD13a. This method of SNP calling was also used for selected ancient samples described in Jones et al.24 Cassidy et al.25, Gunther et al.26, Omrak et al.6 and Olalde et al.27. Genotype calls for these ancient samples were merged with calls from modern samples found in the Human Origins dataset and ancient samples provided in the Mathieson et al.5 dataset which also included genotype calls for previously published ancient samples151617192021222327. To avoid biases caused by post-mortem DNA damage, only transversion sites were used for PCA, ADMIXTURE, f3-statistics and D-statistics.

Principal component analysis

To explore GD13a and other ancient samples in the context of modern variation in Eurasia, we performed PCA with a panel of contemporary populations (196 contemporary populations, 145,004 transversion SNPs). The analysis was carried out using SmartPCA28; the components were loaded on the contemporary populations, and the ancient individuals were projected onto these dimensions (Fig. 1 and Fig. S4).

ADMIXTURE

A clustering analysis was performed using ADMIXTURE version 1.2329, using the full panel of modern and ancient samples described above. SNPs in linkage disequilibrium were thinned using PLINK (v1.07)54 with parameters –indep-pairwise 200 25 0.516, resulting in a set of 116,834 SNPs for analysis. Clusters (K) (2–20) were explored using 3 runs with fivefold cross-validation at each K with different random seeds. The minimal cross-validation error was found at K = 17, but the error already starts plateauing from roughly K = 10, implying little improvement from this point onwards (Fig. S5). The ADMIXTURE proportions are shown in Fig. S6 for all samples and in Fig. 1C for GD13a and selected modern and ancient populations harbouring the component dominant in GD13a.

Outgroup f 3-statistics and D-statistics

Outgroup f3-statistics and D-statistics were performed using the qp3Pop and qpDstat programs from the ADMIXTOOLS package30.

Neighbour-joining tree

We used a custom Matlab script to calculate pairwise pi from genome-wide genotype data using a panel of 22 individuals (from the dataset described above), including GD13a, representative ancient samples, and different modern populations from the same geographic area as GD13a, and generated an unweighted pair group method with arithmetic mean (UPGMA) tree using the seqlinkage function in Matlab’s Bioinformatics Toolbox55.

Runs of homozygosity

In order to examine runs of homozygosity (ROH) we used imputation to infer diploid genotypes in our sample following the method described in Gamba et al.15. We used GATK Unified genotyper44 to call genotype likelihoods at SNP sites in Phase 3 of 1,000 genomes project56 (version 5a downloaded from the BEAGLE website, https://faculty.washington.edu/browning/beagle/beagle.html). Genotype likelihoods were called for alleles observed in the 1,000 Genomes Project and equal likelihoods were set for positions with no spanning sequence data as well as positions where the observed genotype could be explained by deamination. Genotypes were imputed using Beagle 4.0 with default parameters in intervals of 1 Mb57. We imposed a genotype probability threshold of 0.99 (any SNP without a genotype exceeding this threshold had a missing genotype assigned) while converting to PLINK-format genotype data. These data were merged with the dataset used in Jones et al.24 and ROH analysis was carried out as outlined in Gamba et al.15 and Jones et al.24.

Phenotypes of interest

Genes associated with a particular phenotype in modern populations were explored in GD13a. Observed genotypes were called using GATK Unified genotyper44, calling alleles present in Phase 1 of 1,000 genomes dataset58 with base quality ≥20. As many diagnostic markers had 1-fold coverage or less, we also used imputation to infer genotypes at these positions. Imputation was performed as described in section S11, imputing at least 1 Mb upstream and downstream of the SNP position (this interval was reduced for those variants within the first 1 Mb of the chromosome). The Hirisplex prediction model59 was used to explore hair and eye colour (Table S5). For the observed data, if the sample had 1x coverage, the variant was called as homozygous for that allele. Hair and eye colour predictions were confirmed using imputed data.

Additional Information

Accession codes: Raw reads from Ganj Dareh 13a are available for download through the EBI European Nucleotide Archive (ENA) accession number PRJEB13189.

How to cite this article: Gallego-Llorente, M. et al. The genetics of an early Neolithic pastoralist from the Zagros, Iran. Sci. Rep. 6, 31326; doi: 10.1038/srep31326 (2016).