Vegetation History and Archaeobotany https://doi.org/10.1007/s00334-022-00898-7 ORIGINAL ARTICLE Genomic analysis of emmer wheat shows a complex history with two distinct domestic groups and evidence of differential hybridization with wild emmer from the western Fertile Crescent Alice Iob1 · Laura Botigué1 Received: 2 March 2022 / Accepted: 15 August 2022 © The Author(s) 2022, corrected publication 2023 Abstract Triticum turgidum ssp. dicoccoides (wild emmer wheat) was one of the first plants that gave rise to domestic wheat forms in southwest Asia. The details of the domestication of emmer and its early dispersal routes out of southwest Asia remain elusive, especially with regard to its dispersal to the east. In this study, we combine whole genome data from a selection of specimens of modern wild T. turgidum ssp. dicoccoides and domestic T. turgidum ssp. dicoccum (emmer wheats) with a previously published 3,000 year old sample, to explore the phylogenetic relationships between wild and domestic populations of emmer, and especially the early dispersal routes south and eastwards to Africa and Asia, respectively. Our data confirm a marked differentiation between landraces from Europe, the Caucasus and Iran, and those from Africa, the Arabian Peninsula and India, the first group being more closely related to wild emmer from the northern and eastern Fertile Crescent. Gene flow is detected between wild emmer from the western Fertile Crescent and the second group of domestic emmer. These results support a dispersal route from southwest Asia to Africa, the Arabian Peninsula and India. We also observe a lower genetic variability in the wild emmer from the northern and eastern compared with that of the western Fertile Crescent. It is possible that the ancestors of domestic emmer that spread into Egypt still remain to be surveyed and analysed. Investigating the genetic content of ancient samples from Europe, the Caucasus or Iran would be very valuable to determine whether the two distinct types of germplasm arose during history or were already present during the early phases of dispersal. Keywords Emmer wheat · Domestication · Dispersal · Genomics · Ancient DNA Introduction even if some of the details of this process at the species level remain elusive. The domestication and spread of crop plants has aroused the Triticum turgidum ssp. dicoccoides (Asch. & Graebn.) curiosity of the scientific community for a century now, ever Thell. (wild emmer wheat) is a tetraploid mainly self-polli- since the pioneering studies by Nikolai Vavilov, detailed in nating hulled cereal. Its range of distribution covers the Fer- Studies on the origin of cultivated plants (Vavilov 1926). tile Crescent, including Israel, Jordan, southwestern Syria, Technical innovations and methodological improvements Lebanon (the western Fertile Crescent), and southeastern have allowed a re-evaluation of old theories about plant Turkey, northern Iraq and western Iran (the northern and domestication, in which evidence from both archaeobotany eastern Fertile Crescent) (Zaharieva et al. 2010; Özkan et al. and genetics shows domestication as a geographically dif- 2011). It is unclear whether domestic emmer emerged in fused and genetically varied phenomenon in southwest Asia the western or the eastern and northern Fertile Crescent or (Fuller et al. 2014; Allaby et al. 2017; Pankin et al. 2018), independently in both regions. In the site of Ohalo II in the western Fertile Crescent (Kislev et al. 1992; Weiss et al. Communicated by J. Whitlam. 2004) dated to ca. 23,000–21,000 years cal bp, as many as 36% of the emmer rachis remains carry the diagnostic * Laura Botigué scar associated with a non-shattering rachis and thus with laura.botigue@cragenomica.es domestication (Snir et al. 2015). Archaeological assem- 1 Centre for Research in Agricultural Genomics (CRAG), blages with 100% domestic emmer are first found dating to Cerdanyola del Vallès, 08193 Barcelona, Spain 10,600–10,200 years cal bp, from sites in the northern and 13 Vol.:(0123456789) Vegetation History and Archaeobotany eastern Fertile Crescent such as Çayönü, Turkey (van Zeist 5,000 bp emmer reached Ethiopia (Helbæk 1970). Whether and de Roller 2003). Ethiopian emmer is a descendant of Egyptian emmer or if it Most of the studies of genetics and genomics have been reached there through the Iranian highlands and the Arabian based on modern material. Modern domestic landraces Peninsula is not known. (Luo et al. 2007). of emmer are closer to wild emmer accessions from the Triticum turgidum ssp. dicoccum (Schrank ex Schübl.) northern and eastern Fertile Crescent (Avni et al. 2017). Thell. (fully domestic emmer), together with T. monococcum Also, early genetic models based on the similarity between (einkorn), Hordeum vulgare (barley), Pisum sativum (pea), domestic and wild emmers from the northern and eastern Lens culinaris (lentil) and Linum usitatissimum (flax) Fertile Crescent suggested that southeastern Turkey and reached western Iran by 9,000 bp (even though evidence northern Syria were the centres of emmer domestication of management of partially domestic emmer is known (Lev-Yadun et al. 2000; Özkan et al. 2002). More recently, from Chogha Golan as far back as 9,800 bp (Riehl et al. the idea that domestic emmer originated from a single 2013)) northern India (the upper Punjab plain) by the first homogeneous wild population has been questioned half of the 5th millennium bp and southern India during (Jorgensen et al. 2017), with some authors proposing a the 4th millennium bp. While emmer was important in model in which fully domestic emmer wheat emerged from India, it was not introduced into central Asia (north or a population in the northern and eastern Fertile Crescent east of Turkmenistan and to Afghanistan), where only that would carry the genetic background of multiple wild Triticum aestivum (free-threshing wheat) is found in the emmer populations, including those from the western Fertile archaeological record (Stevens et al. 2016). Crescent (Civáň et al. 2013). This would be in agreement Modern landraces of emmer from India were classified with other genetic studies, which find a greater level of by Vavilov (1926) as indostanicum group within subspecies similarity between domestic emmer populations and wild abyssinicum (Ethiopian emmer), suggesting that the populations from the northern and eastern Fertile Crescent dispersal of emmer into Ethiopia and India was somehow (Oliveira et al. 2020), but also a contribution from the connected. Two routes have been proposed to explain this. western Fertile Crescent populations to important domestic One possibility is that the introduction into northwest India haplotypes (Nave et al. 2019). occurred via Iran and Afghanistan (Mani 2004) and then These genetic findings are compatible with spread south into India. The second hypothesis proposes a archaeobotanical evidence of common finds of wild emmer first introduction to India by sea from the Arabian Peninsula, in the western Fertile Crescent from the pre-pottery Neolithic from a population characterized by low genetic diversity A (PPNA) (11.6–10.7 ka cal bp) with increasing proportions (Salunkhe et al. 2013). This theory would be compatible of domestic types from the following periods, whereas in with known trade routes from the Red Sea across the the northern and eastern Fertile Crescent other taxa were Arabian Sea or the Indian Ocean which existed since Greek preferentially consumed. It was not until the middle and and Roman times (Luo et al. 2007). Both possibilities are late pre-pottery Neolithic B (PPNB) (10.2–8.3 ka cal bp), compatible with trade routes connecting India through the that there was a change towards an increased management coast of Iran with the Arabian Peninsula and from there to and consumption of domestic wheats in the northern and Ethiopia. eastern Fertile Crescent, evident from findings from different One limitation of the genetic studies is that they are all archaeological sites (Arranz-Otaegui et al. 2016; Kabukcu based on modern data, and it is not possible to determine et al. 2021). The Iranian site of Chogha Golan represents a whether modern landraces are representative of those present good example of this, with its sequence of over 2,200 years at the time of early dispersal of cultivated emmer. The of plant management there. From the PPNA to mid pre- publication of the genomic sequence of an ancient sample pottery Neolithic B (MPPNB), only a small fraction of the of emmer from Egypt dated to 3,130–3,000 bp (Scott et al. remains from Chogha Golan is represented by wild emmer, 2019) re-opened interesting questions about its dispersal. while other cereals such as Hordeum vulgare (wild barley) This sample resembled present-day landraces from India, are predominant. However, starting from 9,800 cal bp, the Arabian Peninsula and Turkey. Interestingly, this study domestic emmer appears, soon outnumbering other large- showed signals of genetic introgression from the western seeded grasses (Riehl et al. 2013). Fertile Crescent, which possibly occurred when emmer was From southwest Asia, emmer spread towards Europe, cultivated, but before its introduction to Egypt, or during central Asia and Africa. While dispersal of emmer into later interactions, indicating a connection between early Europe has been well studied together with the spread of emmer dispersal eastwards across the Iranian plateau and agriculture (for example, Coward et al. 2008), its dispersal into the Indus valley and also to the southwest into the Nile to the south and east is still a matter of debate (Stevens valley. However, due to the lack of examples of emmer from et al. 2016). In Africa, the first settlements in Egypt date to Africa in the dataset, it was not possible to find out about 7,500–6,650 cal bp (Wendrich and Cappers 2005). Around these early emmer trade routes from this study. 13 Vegetation History and Archaeobotany Here we re-analyse this ancient sample in the context of Alignment and variant calling a more extended modern genomic dataset with representa- tion of emmer landraces from Ethiopia, to elucidate the phy- Main dataset logenetic relationships between domestic and wild emmer populations, and with the ancient sample, with the ultimate We aligned all the samples to the durum wheat genome goal to bring insight into the routes of dispersal of emmer (Triticum durum), using the program bwa v. 0.7.17 (Li (Fig. 1). and Durbin 2009). For each wheat specimen, this process takes all the sequenced DNA fragments (reads) and finds its coordinates (chromosome, position) in the reference Materials and methods genome. We then sorted and indexed the resulting *.bam files using Samtools 1.9 (Danecek et al. 2021). Variant Samples calling (the process of identification of genetic differences between the sequenced specimen and the reference genome) Main dataset was done with the genome analysis toolkit GATK v. 4.1.6 (van der Auwera and O’Connor 2020) and Picard v. 2.22.3 Previously published whole genome sequence data from (Broad Institute 2019), retaining only biallelic SNPs (single 57 emmer accessions were downloaded from https://​bigd.​ nucleotide polymorphisms, which are changes to a single big.​ac.​cn/​searc​h/?​dbId=​gsa&q=​CRA00​1951 (Zhou et al. DNA base, hereafter referred to as “variants”), and hard- 2020a). This dataset comprises 28 wild samples and 29 filtering the dataset following Zhou et al. (2020a) to keep domestic samples covering Europe, western Asia and the only those genetic variants that were identified with high Horn of Africa. confidence. We further removed variants (SNPs) with more than 10% missing information or that were found in Extended dataset one sample only using VCFtools v. 0.1.16 (Danecek et al. 2011) with the commands—max-missing 0.1 and—mac 3. For some analyses, the main dataset is extended to include The total number of variants in this dataset was 66 million an ancient sample from Egypt, 14C dated to 1130–1000 cal (66,097,433 SNPs). We further filtered our dataset for bc (Scott et al. 2019), and four domestic landraces from linkage disequilibrium using plink v. 1.9 whole genome Turkey, Oman and India, sequenced by exome capture (Avni association analysis toolset (Purcell et al. 2007), allowing et al. 2017). Details for the dataset can be found in ESM a maximum ­r2 value of 0.1 calculated in 50 kb windows Table S1, and the geographical distribution of the samples with a step size of 10 kb, reducing the dataset to 4.5 million is shown in Fig. 1. variants (4,444,631 SNPs). Increasing the value of r­ 2 to 0.4 gave the same results in both principal component analysis (PCA) and discriminant analysis of principle components (DAPC). For this reason, we used the filtering for ­r2 0.1 in all analyses. Fig. 1  Map showing the distri- bution of sites of the samples analysed in this study, coloured according to the genetic cluster- ing identified in the analyses. WWFC wild western Fertile Crescent, WNEFC wild north- ern and eastern Fertile Crescent, DNW domestic northwest route of dispersal, DSE domestic southeast route of dispersal. Samples from the WNEFC population come from only two sampling sites (7 samples, 2 sites). Sample W-LBN3 from Lebanon is not represented due to lack of geographical informa- tion 13 Vegetation History and Archaeobotany Extended dataset Population structure of emmer The modern samples from the extended dataset were also PCA aligned to the durum wheat reference genome and processed in the same way as the main dataset. The data from the In order to get a general overview of the relationships ancient samples were processed as in Scott et al. (2019) to between our samples we did a Principal Components adequately account for the characteristics of degraded DNA. Analysis (PCA) on the main dataset, using plink v. 1.9 and Variants were re-called on the whole dataset as previously plotted the results in R v. 4.1.0 (R Core Team 2021). This described and sites present in both datasets were kept, analysis led to the identification of two outliers (samples resulting in 3.6 million (3,689,770 SNPs) variants. After ISR5 and SRB3, ESM Fig. S1). In the absence of additional filtering for linkage disequilibrium (LD) (­ r2 0.1 as above) information about these samples, we decided to exclude and removing transitions to avoid errors from postmortem them from further analyses. damage, the number of variants was ca. 400,000 (433,973 SNPs). DAPC Outgroup Since emmer wheat is a self-fertilising, highly inbred taxon, we also performed a discriminant analysis of principal com- Some of the analyses required the comparison of the ponents (DAPC) on the main dataset to better discriminate genetic data observed in emmer wheat with that of an between groups. DAPC identifies clusters by minimizing the outgroup, a more distantly related group that can serve differences within groups while maximizing the differences as a reference when determining the relationships within between groups. This was done with the R package Adegenet, the ingroup. Since wild emmer is the result of an ancient in R v. 4.1.0 on variants that were at least 25,000 bases apart hybridization event between Triticum urartu (A genome) (VCFtools v. 0.1.16 command—thin 25 000 will take one and Aegilops speltoides (B genome), it is not possible to variant every 25,000 bases) resulting in 319,331 variants. In find a tetraploid outgroup in nature. We circumvented this DAPC the number of retained principle components (PCs) by following Scott et al. (2019) and used publicly available is critical, as retaining too many of them compared with the genomic data of the donors of the A and B genomes of sample size could lead to over-fitting and a subsequent distor- emmer as an outgroup. We downloaded genomic sequences tion of the results. For this reason, we performed the cross of A. speltoides (SAMEA2342530, European Nucleotide validation using the xvalDapc function, and retained the first Archive) representing the B subgenome and T. urartu ten principle components (ESM Fig. S2). (sample A082 from Zhou et al. 2020a, https://​bigd.​big.​ac.​ cn/​searc​h/?​dbId=​gsa&q=​CRA00​1951), representing the ADMIXTURE A subgenome. Paired-end reads were aligned to the durum reference genome. We retained only reads mapping to the This is a maximum likelihood based unsupervised clustering A subgenome for T. urartu and the B subgenome for A. algorithm that estimates the proportions of an established speltoides. Retaining only these reads allowed us to avoid number of ancestries for each specimen in the dataset mis-mapping biases from reads mapping to the wrong (Alexander et al. 2009). Maximum likelihood methods homologue, an issue that is common in wheat studies, due estimate the most probable model given the observed to the high levels of homology between subgenomes. data. Given a certain number of ancestral components, After re-calling variants in the whole dataset, we K, the individuals can be represented as a mixture of kept only the sites that were polymorphic in the emmer such components. In order to determine the best number dataset. This led to the identification of 56 million variants of K for the main dataset, we used cross validation error (56,200,433 SNPs), which after LD filtering reduced to 4 analysis in ADMIXTURE, with values of K from 1 to 10. million variants (4,078,964 SNPs). The best values for K proved to be 3 (cross validation error All the analyses used the dataset which had been filtered 0.58) and 4 (cross validation error 0.58), even if very close for LD, except the calculation of the nucleotide diversity. values were obtained for the values 2, 5 and 6 for K (cross validation error 0.61–0.63). All analyses were done using ADMIXTURE v. 1.3.0 and results plotted using R v 4.1.0. Phylogenetic analysis In order to discover the phylogenetic relationships between the samples in our dataset, we first converted the format 13 Vegetation History and Archaeobotany of the dataset from vcf to relaxed Phylip using the script D statistic vcf2phylip.py, downloaded from https://​github.​com/​edgar​ domor ​tiz/​vcf2p​hylip, and used the resulting phylip file In order to detect gene flow between different populations, as input for RaxML (randomised axelerated maximum we calculated Patterson’s D (ABBA-BABA test) (Green likelihood, Stamatakis 2014) and MEGA XI (molecular et al. 2010) between the two wild groups (WNEFC for evolutionary genetics analysis, Kumar et al. 2018). We northern and eastern Fertile Crescent and WWFC for constructed a phylogenetic tree using maximum likelihood western Fertile Crescent) and the two domestic groups, with RaxML v. 8.2.12 and our main dataset. We ran a DNW for the northwestern route to Europe, the Caucasus, RaxML rapid bootstrap analysis (option -f a), searching for Balkans and Iran and DSE for the southeastern route to the best tree out of 20 runs (option -# 20). For the extended Ethiopia and India in the main dataset. This analysis is dataset, since it included an ancient sample, and in order to based on the fact that, given a known phylogeny with four avoid introducing any bias due to model selection, we used populations (((P1, P2) P3) OUTGROUP) represented as neighbour joining clustering to construct a phylogenetic BBAA, in which A is the ancestral condition (allele) and tree, based on p-distance with 100 bootstrap replicates B the derived one, by analysing different genomic regions using MEGA XI; bootstrapping is a re-sampling method one can obtain a certain number of phylogenies which do for assessing the reliability of the results. We kept only not fit, such as ABBA and BABA (P2 and P3 sharing the transversions (variants that entail a change from a purine derived allele and P1 and P3 sharing the derived allele), (Adenine, Guanine) to a pyrimidine (Cytosine, Thymine) due to incomplete lineage sorting. This is random, and in or vice versa) in order to eliminate the potential mis- the absence of gene flow between populations, the number incorporations from Cytosine to Thymine and from Guanine of ABBAs and BABAs should be equal or not significantly to Adenine related to ancient DNA (aDNA) damage. The different and the D statistic should be zero. On the other resulting phylogenetic tree topology (which is the branching hand, an excess of ABBAs or BABAs and the resulting structure of the tree, indicating the patterns of relatedness deviation of D from 0 is a sign of gene flow between the among taxa) was mainly consistent with that from the main two populations. The statistic was computed using Dsuite v. dataset based on a maximum likelihood approach, which 0.4 r41 (Malinsky et al. 2021). Significant results are defined was further verified by constructing a new tree for the main by an absolute Z-score for relationship to the average larger dataset samples with neighbour-joining (ESM Fig. S3). than 3. Genetic diversity TreeMix The amount of genetic diversity in each group was calculated TreeMix v. 1.13 (Pickrell and Pritchard 2012) builds a as nucleotide diversity (pi) using VCFtools v. 0.1.16. For phylogenetic tree using maximum likelihood and allowing groups with few samples it is possible that the genetic for gene flow between populations. We used this on the variability within the group is not representative of that of modern dataset after converting the *.vcf data to TreeMix the real population. Since the northern and eastern Fertile input format, using the script vcf2treemix.sh, downloaded Crescent group consisted of only seven samples compared to from https:// ​ g ithub. ​ c om/ ​ s peci ​ a tion ​ g enom ​ i cs/ ​ s crip ​ t s/​ other group sizes of > 13, bootstrapping was used to enable blob/​master/​vcf2t​reemix.​sh, and used the output file as comparison between groups, in which four random samples input for TreeMix. We tested for 0 to 3 migration edges, were extracted from each group to calculate nucleotide which represent events of migration (gene flow) from one diversity (pi), and the process was repeated ten times. The population to another (represented as arrows in the output) averaged value of pi from all the subsets within a population (-m 0 to 3). Applying bootstrap validation with blocks is taken as its value for the whole group. of 500 SNPs (-bootstrap -k 500), we found that the only meaningful trees were the ones with 0 or 1 migration edges Analysis of the southeastern dispersal route (more edges show gene flow from the outgroup). We also and gene flow used TreeMix on the extended dataset, this time without the sample size correction (applying -noss) and, as in the We analysed the genetic make-up of populations of domestic tree using neighbour joining, we kept transversions only emmer that dispersed to the south towards Africa and the and filtered with linkage disequilibrium (LD). The results Arabian Peninsula, and also to the east towards Asia and were plotted in R using the script plotting funcs.R, which is India, in an attempt to detect hybridization events, especially included in TreeMix. with the wild western Fertile Crescent (WWFC) population. 13 Vegetation History and Archaeobotany Results The resulting clustering is used for finding out about the possible causes of these observed patterns, for example, Re‑processing the whole genome sequence data geography. The PCA led to the identification of four groups of samples, differentiating wild from domestic and grouping The samples published by Zhou et al. (2020a) were samples according to geographical patterns (ESM Fig. S1). aligned to the bread wheat reference genome. In order to Two outliers were detected, most probably due to inaccurate be able to capture as much genetic variation as possible, passport information, the meta-information such as origin, we re-processed them to align them to the durum wheat species, data of collection, that has been associated with reference genome (Maccaferri et al. 2019). The geographical this accession. After excluding these outliers, the samples distribution of the samples is shown in Fig. 1, with the were grouped into clusters for further analysis. Wild emmer samples coloured according to the results from discriminant samples from the western Fertile Crescent are shortened analysis of principal components (DAPC). to WWFC and those from the northern and eastern Fertile After re-processing this genomic sequenced data, Crescent to WNEFC in the following text. Domestic samples we had a high quality dataset of 66 million variants from Ethiopia and Oman are called DSE as in domestic (66,097,433 SNPs) for 55 samples. After filtering for linkage southeast, to refer to the southeast route, whereas samples disequilibrium (LD) the dataset was reduced to 4.4 million from Europe, the Caucasus, Balkans and Iran are called variants (4,444,631 SNPs). For those analyses that required DNW as in domestic northwest, since they broadly represent an outgroup, we kept only the sites that were polymorphic the northwest route of emmer dispersal. We discuss the in the emmer dataset. This led to the identification of 56 intriguing grouping of the domestic Iranian landraces in the million variants (56,200,433 SNPs), which after LD filtering DNW group, below. reduced to 4 million (4,078,964 SNPs) variants. For the This sample clustering was confirmed by DAPC (Fig. 2a). analysis of the southeastern migration route, adding the Interestingly, not all domestic samples are equally close to data from exome capture and the ancient samples resulted the WNEFC group, which would be expected if there had in a dataset of 3.7 million (3,689,770 SNPs) variants, which been a single domestication and dispersal event. DNW after filtering for LD and removing transitions to avoid landraces appear closer to the WNEFC group in the PCA errors related to postmortem damage resulted in ca. 400,000 and are on the same x-axis in the DAPC, whereas the DSE (433,973 SNPs) variants. For these variants the mean depth specimens are more distant from the WNEFC cluster in both of the whole genome sequenced (WGS) samples (including analyses. outgroup) calculated with VCFtools is 5.2 × , that is, on To investigate the population structure in more detail average each of the positions where these variants are found, and identify shared genetic components between samples, is covered 5.2 times by the sequenced DNA fragments. We we applied ADMIXTURE clustering analysis, which called variants on the exome capture samples and on the models shared ancestries between individuals. Values of ancient sample only on sites known to be polymorphic in K (representing the number of ancestries) between 2 and the main dataset, allowing for a minimum depth of 1 × . The 6 were tested (ESM Fig. S4), and cross-validation errors covered sites for the ancient sample have a mean depth of determined that modelling three or four ancestries provided 1.47 ×. a best fit with the data. With a K value of 2, the dataset is divided into two ancestries, one characterizing all WWFC Population structure samples (shown in pink) and the other all domestic samples (blue), with WNEFC emmer specimens having a varying To get a general overview of the relationships between proportion of the two ancestries. When an additional modern emmer landraces, we used principal component ancestral component is allowed (K = 3), the domestic DSE analysis (PCA) and discriminant analysis of principal group is assigned to this new component (yellow), the components (DAPC). Both these methods are used to WWFC group keeps its ancestry (pink) and the WNEFC and infer population structure without previous knowledge, DNW groups are largely represented by the other ancestral by determining the number of observed clusters in the component (blue) (Fig. 2b). It is interesting to note that, in dataset, based on principal components, which are the ones agreement with the PCA results, allowing a third ancestry that explain most of the genetic variability in the dataset. results in the differentiation between DSE and the DNW- In DAPC the differences between groups are maximized, WNEFC group, which now appears homogeneous. At K = 4, while the differences within them are minimized. These the WNEFC and DNW groups are further differentiated, methods are ideal for measuring the degree of differentiation with domestic specimens from DNW being largely assigned between the samples when using genetic information only. to a new ancestral component (green) (Fig. 2b). Very low levels of admixture are found between the wild populations, 13 Vegetation History and Archaeobotany Fig. 2  a Discriminant analy- sis of principal components (DAPC). Groups shown by their colours are arranged accord- ing to the PCA results and confirmed by this analysis. b ADMIXTURE analysis for the best values of K. Upper panel, K = 3; lower panel, K = 4. c Maximum likelihood tree, modern whole genome data only. The tree was created with RaxML, using a fast search to find the best tree out of 20 runs. Bootstrap values are reported at the top of the nodes and only a few samples from DNW appear to be admixed a clear geographical pattern. The DSE cluster is subdivided between WNEFC, DNW and DSE components, perhaps as in three sub-clades, the most divergent being the one with a result of the under-representation of the Mediterranean the samples from Oman. In the WWFC cluster, on the landraces. Overall, these results are in agreement with the other hand, the subclades mirror the geographical origins DAPC results. of the samples. Judging by their scattered position within We finally further explored the genetic affinities between the WWFC group, emmer from Syria seems to be the most wild and domestic modern emmer specimens by constructing diverse. This pattern was replicated with the neighbour- a phylogenetic tree using maximum likelihood (Fig. 2c). joining tree (ESM Fig. S3). Overall, the specimens are grouped into clades with the same Overall these results confirm that WWFC is the most grouping, such as that when using PCA and ADMIXTURE. differentiated population within the dataset and they rein- The first node in the tree divides WWFC from all other force the observation of the genetic differences between the populations, while the WNEFC samples cluster together DNW and DSE landraces, as well as their varying affinity and are an outgroup to all the domestic samples (with the to the WNEFC samples. It is unclear why the DSE modern exception of the only wild Iranian emmer specimen, that landraces are so different from this wild population. These appears as the closest relative to all other domestic emmers). results do not suggest that domestic emmer emerged in the Within the domestic cluster, samples maintain the DNW northeastern and western Fertile Crescent independently, and DSE groupings, even though samples from Spain and since the DSE samples are not closer to WWFC. A wide the UK appear as outliers to the DNW clade, not following range of modern wheat landraces were studied to investigate 13 Vegetation History and Archaeobotany during the early dispersal of the group. In such a situa- tion, the low genetic diversity within the group would have increased its genetic differences from other groups, delet- ing the genetic signature linking this group to the WNEFC emmer wheat ancestor, a process called genetic drift. In order to determine whether genetic drift caused this differentiation in the DSE group, we calculated nucleo- tide diversity (Fig. 3). As expected, the wild populations show higher levels of nucleotide diversity, even if, intrigu- ingly, the WWFC group shows levels of diversity (pi 0.19) almost twice those of WNEFC (pi 0.10). Domestic groups show lower levels of nucleotide diversity than WWFC and WNEFC, but similar levels of genetic diversity among themelves (DSE pi 0.06, DNW 0.07). As the genetic diver- sity of the DSE group is comparable to that of DNW, this evidence does not suggest a strong genetic bottleneck, a sharp reduction in the size of the population of the DSE ancestor, pointing to other possible explanations for its differentiation. Gene flow Another explanation for the distinctiveness of the DSE group could lie in a different level of contribution of the WWFC populations to this group. This could have occurred during the domestication process or the early dispersal of the domestic emmers into Africa. In order to investi- gate if the domestic landraces from the DSE group show signs of admixture with wild specimens from the west- ern Fertile Crescent, we calculated Patterson's D statistic (ABBA–BABA) (Green et al. 2010) for the deviation from the expected ratio of allele sharing between WWFC and the two domestic groups, plus the outgroup, using the phylo- genetic tree: (((DNW, DSE) WWFC) OUTGROUP). The results show indeed an excess of allele sharing between WWFC and DSE (tree (((DNW; DSE) WWFC) OUT- GROUP) D = 0.106, Z score = 3.15), compatible with gene flow from the wild population into the domestic one (Fig. 4 and ESM Table S2). Replacing the WWFC population by the WNEFC one yielded no evidence of gene flow between Fig. 2  (continued) any of the domestic groups (topology (((DSE; DNW) WNEFC) OUTGROUP) D = 0.074, Z score = 2.03), consist- the allelic variability in the genes responsible for a non-shat- ent with both domestic populations having more affinities tering rachis (TtBTR1A and TtBTR1B) (Nave et al. 2019). with WNEFC than with WWFC. The authors found no diversity within the domestic samples, In light of the evidence of hybridization observed suggesting that the fully domestic phenotype, to which wild between the WWFC and DSE groups, we constructed a populations WWFC and WNEFC probably contributed, had maximum likelihood phylogenetic tree with TreeMix, a single origin. The observed results could be explained if which models gene flow by introducing edges of migration the ancestors of the modern DSE group had experienced between populations (Fig. 5). As expected, under no migra- a strong founder effect, the loss of genetic variability that tion, TreeMix constructed the same tree topology previously occurs when a new population is established by a very small obtained. However, the plot of the residuals showed that the number of individuals from a larger population, possibly tree does not properly fit the data (ESM Fig. S5a). When 13 Vegetation History and Archaeobotany vs. ln(likelihood) = − 224.897 for 0 and 1 migration events, respectively) and lower standard error, ± 43.4 vs. ± 15.2 for 0 and 1 migrations events, respectively (ESM Fig. S5a, b), confirming the Patterson’s D statistic results. In order to better characterize the nature of the DSE group, we included publicly available data on emmers from Turkey, Oman and India (referred to as the Indian Ocean group) that also differed from domestic samples from Europe and the Caucasus (Avni et al. 2017). We note that the inclu- sion of modern landraces from Turkey could potentially dif- ferentiate between whether hybridization happened during the domestication process or during the early dispersals from the northern and eastern Fertile Crescent southwards. We also included the genomic data of a 3,000 year old ancient Egyptian emmer sample (Scott et al. 2019) to determine whether the genetic make-up of this domestic group is mod- Fig. 3  Nucleotide diversity (pi) in different populations, showing that the WWFC population is the most diverse, while WNEFC is less so. ern or if it already existed in the past. We first constructed a Domestic populations, DNW and DSE, have similar levels of diver- neighbour joining tree to study how the specimens from this sity extended dataset were related to those from the main one. The extended dataset does not change the arrangement of the tree obtained with respect to the main one (Fig. 6). Interest- ingly, the samples D-TUR1 and D-TUR2 from Turkey clus- ter together with the samples from Ethiopia, while D-OMN3 from Oman and D-IND1 from India cluster with the other samples from Oman. The ancient sample UC10164 from Egypt is placed as an outgroup to the DSE group and shows a quite long branch, which could be due to the low coverage of the sample or a slightly different genetic make-up. We next constructed a tree that allowed for gene flow between its branches (Fig. 7). The new emmer specimens from India, Oman and Turkey were merged into the DSE Fig. 4  Schematic representation of Patterson’s D (ABBA-BABA) group (now DSE + IO). One edge of migration in the tree statistical test results; the arrow shows the direction of gene flow between WWFC and DSE continued to support gene flow from WWFC, this time to the branch of the common ancestor of the ancient sample allowing for an edge of migration, gene flow is identified and the DSE + IO group, indicating some mixture in both from the WWFC population branch to the DSE leaf, giving of them. The likelihood of a tree with one edge of migration a tree with increased likelihood (ln(likelihood) = − 4477.06 is higher than one with no edges, while the Standard Error Fig. 5  TreeMix maximum likelihood analysis of the modern whole genome dataset allowing for 0 (left panel) and 1 (right panel) migration events, showing gene flow between WWFCSL and DSE (arrow). The colour of the arrow is proportional to the amount of migration 13 Vegetation History and Archaeobotany Fig. 6  Neighbour joining phy- logenetic tree of the extended dataset, including the Indian Ocean samples and the ancient Egyptian sample. Bootstrap 100 replicates, all nodes have boot- strap values are above 90 (SE) is lower (ESM Fig. S5c-d). The (ln(likelihood) with 0 migration events = -4.62704, SE =  ± 8.4; ln(likelihood) with 1 migration event = 159.731, SE =  ± 3.4). Discussion There are many aspects surrounding the emergence of domestic emmer and its dispersal that even after decades of study remain elusive. While archaeobotanical evidence shows that the earliest domestic emmer appeared in sites in the western Fertile Crescent, the earliest assemblages with only domestic emmer are found in the northern and eastern Fertile Crescent. Genomic data has always supported the view that modern domestic emmer is clearly closer to emmer Fig. 7  TreeMix maximum likelihood analysis of the entire data- from the northern and eastern Fertile Crescent. Regarding set including the Indian Ocean and the ancient sample, allowing for its dispersal, it is not clear how and when emmer spread east 0 (left) and 1 (right) migration events, showing gene flow between and south. How it arrived into Ethiopia (whether from Egypt WWFC and the node to modern DSE + IO and the ancient sample or through the Arabian Peninsula) or into India (considering (arrow). The colour of the arrow is proportional to the amount of migration that only free-threshing wheats are found in central Asia) is so far unknown. 13 Vegetation History and Archaeobotany Our analysis of a comprehensive collection of mod- would have needed a second adaptation, to high altitudes, ern emmer landraces from the southern and eastern route during the dispersal from Egypt to the high plains of Ethio- together with an ancient specimen from Egypt, provides pia. Introduction of emmer into India would have occurred some new insights into these events, but also raises new either through the sea routes, or through Iran but without questions. It is now clear that modern emmer landraces can introduction of the germplasm to this area. Another option is be clearly differentiated into two groups, at least in our data- that it was introduced into Iran, and that two different types set, in which landraces from other regions such as the Medi- of emmer were cultivated in this region in the past, but only terranean area are under-represented. These groups can be one persisted into modern times. In a third scenario, modern assigned to two geographical regions, loosely reflecting the landraces from the DSE + IO group might exist in Iran but proposed dispersal routes. One includes Europe and the Cau- have not been studied. Finally, we cannot exclude the possi- casus, and the other Africa and the Arabian Peninsula. Nota- bility that the genetic make-up of the DNW group, including bly, modern landraces from Iran are grouped with the ones the modern Iranian samples, is the product of a much recent from Europe and the Caucasus, which is even more puzzling dispersal into these areas, perhaps during the last millen- in light of the fact that modern landraces from India are nium. The study of ancient samples from this region would grouped with those from the Arabian Peninsula and Africa. be invaluable to investigate these different scenarios. The similarity between the ancient Egyptian sample and this The clear genetic differences between the DNW and DSE DSE + IO group from Africa, the Arabian Peninsula and groups raises another important question. In this paper we south Asia demonstrates that the genetic make-up of these have ruled out that this differentiation is driven by drastic modern landraces has not changed much, at least over the founder effects, since the genetic diversity in the two groups past 3,000 years. The clustering of the landraces from India is similar. If this differentiation reflects a past genetic struc- and Oman is in agreement with known trade routes con- ture, this would suggest that early dispersal of domestic necting southern Asia and Africa (Cuny and Mouton 2009). emmer occurred from two different germplasm sources, However, with the current data it is not possible to know opening the possibility that the full domestic phenotype was whether modern Indian landraces are descendants of emmer present in markedly different populations before its disper- which came there during early dispersal or whether they sal. This would be compatible with a diffused and protracted arrived from the Arabian Peninsula at a later stage. domestication process that may be reflected by the abun- The differentiation of these landraces from those from dance of domestic emmer remains in archaeological sites Iran poses a perhaps bigger mystery. Indeed, Iran and par- throughout southwest Asia over thousands of years. Also, ticularly the surroundings of the Zagros mountains experi- analysis of the genes diagnostic for domestication, the TtBtr enced the same change towards the consumption of increas- loci, show that all modern domestic landraces carry the same ing proportions of domestic emmer during the middle and haplotype (the same variants in a sequence of DNA) (Nave late PPNB, at the expense of other cereals (Riehl et al. 2013), et al. 2019), as do the samples studied in this dataset. The as did other sites in the northern and eastern Fertile Cres- immediate wild ancestor of one of the domestic haplotypes cent, reflecting similar patterns of early agriculture. Without is almost ubiquitous among wild emmer, both from the ancient samples from Europe, the Caucasus or Iran, it is not north-eastern and the western Fertile Crescent, so, even if possible to know if the genetic make-up of the DNW group a single origin is the simplest hypothesis, it is nonetheless is ancient or more recent. If this group represents an ancient possible that the haplotype arose independently in different genetic structure and assuming that the Iranian landraces populations. A more careful examination of the genetic com- which have been studied reflect the diversity of the whole position around this locus could help test this hypothesis. country, this would show that domestic emmer soon replaced However, if this were the case, the key question that would other wheats in the most eastern parts of the Fertile Cres- remain to be answered is how did the same fully domestic cent, reaching Iran around 9,000 bp. It then spread north and phenotype emerge in different wheat germplasms? west, but without reaching south Asia through Iran. Triticum Another interesting aspect that we have found in this study aestivum (free-threshing wheat) would have replaced T. tur- is the gene flow from WWFC to the ancestor of the ancient gidum ssp. dicoccum (hulled emmer wheat) in this route, Egyptian sample and the DSE + IO group. The presence of causing the dispersal of T. aestivum through the inner Asian two modern landraces from Turkey within the DSE + IO mountain corridor (Zhou et al. 2020b), eventually reaching group raises the possibility that gene flow from wild emmer China (Stevens et al. 2016). In parallel, domestic emmer from the western Fertile Crescent did not occur during early would also have first spread south to Egypt, then to Ethio- dispersal to the south, but rather before then. These results, pia, through the Arabian Peninsula and finally from there together with archaeological evidence, point to an early to India, where its first evidence dates to 4,700–4,500 bp contribution from the western Fertile Crescent and possibly (Stevens et al. 2016). Even though the genomic data of our even from the first domesticated emmer to at least part of dataset supports this route, it must be noted that this route the modern gene pool. Also, the clustering of these modern 13 Vegetation History and Archaeobotany landraces from Turkey with those from Ethiopia may suggest evidence, indicate the origin for the domestic gene pool from the existence of different wheat germplasms before their an admixture of wild populations, followed by a process of dispersal. However, if this was so, we would expect the this ancestral group becoming more wild type or feral, that Turkish samples to be sister clades to the DSE + IO group explains the high similarity of the WNEFC population to the on the phylogenetic tree, and not the ancient Egyptian domesticated emmer gene pool. Considering our dataset, this sample. However, the low genomic coverage obtained would mean that the extant feral population was originally with the sequencing of the ancient sample could explain its derived from the DNW, not the DSE population. However, position in the tree. On the other hand, we cannot rule out our results are insufficient to make any speculation about a recent re-introduction of DSE + IO germplasm to Turkey, this, and further analyses are needed to shed light on this perhaps during the last millennium. If gene flow did occur possible sequence of events. during emmer dispersal to the south, the wild component Supplementary Information The online version contains supplemen- of the wheat from the western Fertile Crescent would have tary material available at https://d​ oi.o​ rg/1​ 0.1​ 007/s​ 00334-0​ 22-0​ 0898-7. perhaps provided a basis for its adaptation to hot and dry environments, very different to the climate of the mountain Author contributions LB planned and designed the research. AI slopes in Turkey (Özkan et al. 2011). In any case, more conducted the research. AI and LB wrote the manuscript. ancient and modern samples from these regions would be Funding A.I. is a FPI fellow (PRE2018-083529). L. B. is a Ramón needed to accurately understand this matter. y Cajal Fellow. (RYC2018-024770-I) both fellowships funded by the Overall, our results reveal that modern landraces from Ministerio de Ciencia e Innovación—Agencia Estatal de Investigación/ Ethiopia, the Arabian Peninsula and Africa have a genetic Fondo Social Europeo. We acknowledge financial support from the Spanish Agencia Estatal de Investigación (Ministry of Science and structure that is at least 3,000 years old, very similar to Innovation-State Research Agency) (AEI), through the “Severo Ochoa emmer from Pharaonic Egypt, and most probably similar Programme for Centres of Excellence in R&D” SEV‐2015‐0533 and to the genetic component of the first domestic emmer that CEX2019-000902-S. This work was also supported by the CERCA spread from southwest Asia to Egypt. While it is not possible programme by the Generalitat de Catalunya. to establish whether the modern Indian landraces are Open Access This article is licensed under a Creative Commons descendants of the first emmer that arrived into the region, Attribution 4.0 International License, which permits use, sharing, their similarity with modern landraces from Oman supports adaptation, distribution and reproduction in any medium or format, a theory by which domestic emmer spread from southwest as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate Asia into Egypt, from there to Ethiopia, crossing the Arabian if changes were made. The images or other third party material in this Peninsula and eventually reaching India, perhaps through article are included in the article's Creative Commons licence, unless sea trade. Our data suggest that the emmer that arrived in indicated otherwise in a credit line to the material. If material is not Iran was probably from a different germplasm than the one included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted that was introduced into India. While the analysis of the use, you will need to obtain permission directly from the copyright gene responsible for the brittle rachis trait supports a single holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​ origin for all domestic emmer landraces, we do find two very org/​licen​ses/​by/4.​0/. different germplasms in modern emmer landraces. Ancient data of the DNW germplasm would be needed to confirm how old is the genetic structure of these landraces. References Finally, some of the results of this study show the need for further analyses. The low genetic variability of the Alexander DH, Novembre J, Lange K (2009) Fast model-based estima- tion of ancestry in unrelated individuals. Genome Res 19:1,655– WNEFC groups highlights the importance of adopting a 1,664. https://​doi.​org/​10.​1101/​gr.​094052.​109.​vidual more meaningful sampling strategy that shows whether this Allaby RG, Stevens C, Lucas L, Maeda O, Fuller DQ (2017) Geo- characteristic is real, or an artefact if perhaps this group is graphic mosaics and changing rates of cereal domestication. under-represented in the number of samples from the area. Philos Trans R Soc B 372:1–10. https://​doi.​org/​10.​1098/​rstb.​ 2016.​0429 It remains to be confirmed whether gene flow between wild Arranz-Otaegui A, Colledge S, Zapata L, Teira-Mayolini LC, Ibáñez JJ emmer from the western Fertile Crescent and domestic (2016) Regional diversity on the timing for the initial appearance emmer of the DSE + IO group occurred during the early part of cereal cultivation and domestication in southwest Asia. Proc of the dispersal of emmer or before then. More information Natl Acad Sci USA 113:14,001–14,006. https://​doi.​org/​10.​1073/​ pnas.​16127​97113 on these processes would allow testing whether domestic Avni R, Nave M, Barad O et al (2017) Wild emmer genome architec- emmer from the DSE + IO group could have emerged from ture and diversity elucidate wheat evolution and domestication. a now extinct wild population in the northern and eastern Science 357(6346):93–97. https://​doi.​org/​10.​1126/​scien​ce.​aan00​ Fertile Crescent and hybridized with wild emmer from the 32 Broad Institute (2019) Picard Toolkit. https://​broad​insti​tute.​github.​io/​ western Fertile Crescent when it dispersed towards the south. picard/ Moreover, several studies, according to the archaeological 13 Vegetation History and Archaeobotany Civáň P, Ivaničová Z, Brown TA (2013) Reticulated origin of domesti- Oliveira HR, Jacocks L, Czajkowska BI, Kennedy SL, Brown TA cated emmer wheat supports a dynamic model for the emergence (2020) Multiregional origins of the domesticated tetraploid of agriculture in the fertile crescent. PLoS ONE 8:e81955. https://​ wheats. PLoS ONE 15:e0227148. https://​doi.​org/​10.​1371/​journ​ doi.​org/​10.​1371/​journ​al.​pone.​00819​55 al.​pone.​02271​48 Coward F, Shennan S, Colledge S, Conolly J, Collard M (2008) The Özkan H, Brandolini A, Schäfer-Pregl R, Salamini F (2002) AFLP spread of Neolithic plant economies from the Near East to north- analysis of a collection of tetraploid wheats indicates the origin west Europe: a phylogenetic analysis. J Archaeol Sci 35:42–56. of emmer and hard wheat domestication in southeast Turkey. Mol https://​doi.​org/​10.​1016/j.​jas.​2007.​02.​022 Biol Evol 19:1,797–1,801 Cuny J, Mouton M (2009) La transition vers la période sassanide dans Özkan H, Willcox G, Graner A, Salamini F, Kilian B (2011) Geo- la péninsule d’Oman: chronologie et modes de peuplement. In: graphic distribution and domestication of wild emmer wheat Schiettecatte J, Robin CJ (eds) L’Arabie à la veille de l’Islam: (Triticum dicoccoides). Genet Resour Crop Evol 58:11–53. https://​ Bilan Clinique. De Boccard, Paris, pp 91–133 doi.​org/​10.​1007/​s10722-​010-​9581-5 Danecek P, Auton A, Abecasis G et al (2011) The variant call format Pankin A, Altmüller J, Becker C, von Korff M (2018) Targeted rese- and VCFtools. Bioinformatics 27:2,156–2,158 quencing reveals genomic signatures of barley domestication. Danecek P, Bonfield JK, Liddle J, etal, (2021) Twelve years of SAM- New Phytol 218:1,247–1,259. https://​doi.​org/​10.​1111/​nph.​15077 tools and BCFtools. GigaScience 10:1–4. https://d​ oi.o​ rg/1​ 0.1​ 093/​ Pickrell JK, Pritchard JK (2012) Inference of population splits and gigas​cience/​giab0​08 mixtures from genome-wide allele frequency data. PLoS Genet Fuller DQ, Denham T, Arroyo-Kalin M et al (2014) Convergent evolu- 8:e1002967. https://​doi.​org/​10.​1371/​journ​al.​pgen.​10029​67 tion and parallelism in plant domestication revealed by an expand- Purcell S, Neale B, Todd-Brown K et al (2007) PLINK: a tool set for ing archaeological record. Proc Natl Acad Sci USA 111:6,147– whole-genome association and population-based linkage analyses. 6,152. https://​doi.​org/​10.​1073/​pnas.​13089​37110 Am J Hum Genet 81:559–575 Green RE, Krause J, Briggs AW et al (2010) A draft sequence of the R Core Team (2021) R: a language and environment for statistical com- Neandertal genome. Science 328(5979):710–722. https://​doi.​org/​ puting. R Foundation for Statistical Computing, Vienna, https://​ 10.​1126/​scien​ce.​11880​21 www.r-​proje​ct.​org/ Helbæk H (1970) The plant husbandry of Hacilar. In: Mellaart J (ed) Riehl S, Zeidi M, Conard NJ (2013) Emergence of agriculture in the Excavations at Hacilar. University Press, Edinburgh, pp 189–244 foothills of the Zagros mountains of Iran. Science 341(6141):65– Jorgensen C, Luo M-C, Ramasamy R et al (2017) A high-density 67. https://​doi.​org/​10.​1126/​scien​ce.​12367​43 genetic map of wild emmer wheat from the Karaca Dağ Region Salunkhe A, Tamhankar S, Tetali S, Zaharieva M, Bonnett D, Tre- provides new evidence on the structure and evolution of wheat thowan R, Misra SC (2013) Molecular genetic diversity analy- chromosomes. Front Plant Sci 8:1–13. https://​doi.​org/​10.​3389/​ sis in emmer wheat (Triticum dicoccon Schrank ) from India. fpls.​2017.​01798 Genet Resour Crop Evol 60:165–174. https://​doi.​org/​10.​1007/​ Kabukcu C, Asouti E, Pöllath N, Peters J, Karul N (2021) Pathways s10722-​012-​9823-9 to plant domestication in Southeast Anatolia based on new data Scott MF, Botigué LR, Brace S et al (2019) A 3,000-year-old Egyp- from aceramic Neolithic Gusir Höyük. Sci Rep 11:1–15. https://​ tian emmer wheat genome reveals dispersal and domestica- doi.​org/​10.​1038/​s41598-​021-​81757-9 tion history. Nat Plants 5:1,120–1,128. https://​doi.​org/​10.​1038/​ Kislev M, Nadel D, Carmi I (1992) Epipalaeolithic (19,000 BP) cereal s41477-​019-​0534-5 and fruit diet at Ohalu II, Sea of Galilee, Israel. Rev Palaeobot Snir A, Nadel D, Groman-Yaroslavski I, Melamed Y, Sternberg M, Palynol 73:161–166 Bar-Yosef O, Weiss E (2015) The origin of cultivation and proto- Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: weeds, long before Neolithic farming. PLoS ONE 10:e0131422. molecular evolutionary genetics analysis across computing plat- https://​doi.​org/​10.​1371/​journ​al.​pone.​01314​22 forms. Mol Biol Evol 35:1,547–1,549 Stamatakis A (2014) RAxML version 8: a tool for phylogenetic Lev-Yadun S, Gopher A, Abbo S (2000) The cradle of agriculture. analysis and post-analysis of large phylogenies. Bioinformatics Science 288(5471):1,602–1,603 30:1,312–1,313 Li H, Durbin R (2009) Fast and accurate short read alignment with Stevens CJ, Murphy C, Roberts R, Lucas L, Silva F, Fuller DQ (2016) Burrows—Wheeler transform. Bioinformatics 25:1,754–1,760. Between China and South Asia: a middle Asian corridor of crop https://​doi.​org/​10.​1093/​bioin​forma​tics/​btp324 dispersal and agricultural innovation in the Bronze Age. Holocene Luo M-C, Yang Z-L, You FM, Kawahara T, Waines JG, Dvorak J 26:1,541–1,555. https://​doi.​org/​10.​1177/​09596​83616​650268 (2007) The structure of wild and domesticated emmer wheat pop- Van der Auwera GA, O’Connor BD (2020) Genomics in the cloud: ulations, gene flow between them, and the site of emmer domesti- using Docker, GATK, and WDL in Terra, 1st edn. O’Reilly Media cation. Theor Appl Genet 114:947–959. https://​doi.​org/​10.​1007/​ Inc, Sebastopol s00122-​006-​0474-0 Van Zeist W, de Roller GJ (2003) The Çayönü archaeobotanical record. Maccaferri M, Harris NS, Twardziok SO et al (2019) Durum wheat In: van Zeist W (ed) Reports on Archaeobotanical Studies in the genome highlights past domestication signatures and future Old World. The Groningen Institute of Archaeology. University improvement targets. Nat Genet 51:885–895. https://​doi.​org/​10.​ of Groningen, Groningen, pp 143–166 1038/​s41588-​019-​0381-3 Vavilov N (1926) Studies on the origin of cultivated plants. Institut Malinsky M, Matschiner M, Svardal H (2021) Dsuite—Fast D-statistics Botanique Appliquée et d’Amélloration des Plantes, Leningrad and related admixture evidence from VCF files. Mol Ecol Resour Weiss E, Kislev ME, Simchoni O, Nadel D (2004) Small-grained wild 21:584–595. https://​doi.​org/​10.​1111/​1755-​0998.​13265 grasses as staple food at the 23 000-year-old site of Ohalo II, Mani BR (2004) Further evidence on Kashmir Neolithic in the light Israel. Econ Bot (suppl) 58:S125–S134 of recent excavations at Kanishkapura. J Interdiscip Stud Hist Wendrich WZ, Cappers RTJ (2005) Egypt’s earliest granaries: evidence Archaeol 1:137–142 from the Fayum. Egyptian Archaeol 27:12–15 Nave M, Avni R, Çakır E et al (2019) Wheat domestication in light Zaharieva M, Ayana NG, Al Hakimi A, Misra SC, Monneveux P (2010) of haplotype analyses of the Brittle rachis 1 genes (BTR1-A and Cultivated emmer wheat (Triticum dicoccon Schrank), an old crop BTR1-B). Plant Sci 285:193–199. https://​doi.​org/​10.​1016/j.​plant​ with promising future: a review. Genet Resour Crop Evol 57:937– sci.​2019.​05.​012 962. https://​doi.​org/​10.​1007/​s10722-​010-​9572-6 13 Vegetation History and Archaeobotany Zhou X, Yu J, Spengler RN et al (2020a) 5,200-year-old cereal grains Publisher's Note Springer Nature remains neutral with regard to from the eastern Altai Mountains redate the trans-Eurasian jurisdictional claims in published maps and institutional affiliations. crop exchange. Nat Plants 6:78–87. https://​doi.​org/​10.​1038/​ s41477-​019-​0581-y Zhou Y, Zhao X, Li Y et al (2020b) Triticum population sequencing provides insights into wheat adaptation. Nat Genet 52:1,412– 1,422. https://​doi.​org/​10.​1038/​s41588-​020-​00722-w 13