1. Introduction

Date palm (Phoenix dactylifera) is a keystone tree species in the oasis agrosystems of North Africa and Southwestern Asia (Tengberg, 2012; Terral et al., 2012). Date palms are a key part of the history and culture of these regions. The genus Phoenix L. includes 14 species (Govaerts and Dransfield, 2005) traditionally distributed in the Old World from the Canary and Cape Verde islands in the Atlantic Ocean, throughout Africa, Madagascar and Asia, reaching Sumatra, Taiwan (China) and the Philippines in the East. There are very few studies on the phylogenetic analysis of Phoenix and each of them provided different answers (Barrow, 1998; Pintaud et al., 2010; Pintaud et al., 2013; Torres et al., 2018). A consensus based on the chloroplastic loci support a P. dactylifera clade, which includes P. dactylifera, Phoenix sylvestris, and Phoenix atlantica (Pintaud et al., 2010). However, phylogenetic analysis based on individual nuclear genes (e.g. CYP703, LOG, cytidine deaminase) from 14 Phoenix species suggested a P. dactylifera subclade with P. dactylifera, Phoenix theophresti, and P. atlantica (Torres et al., 2018). So far, there is no comprehensive phylogenetic analysis based on either whole organellar genome assembly or whole nucleargenome assemblies of all members of the Phoenix genus.

The primary center of diversity for Phoenix is an area from India to Indochina, where eight species are found (Pintaud et al., 2010). According to archaeological data, date palm was first cultivated in ~5000-3000 BC near the Persian Gulf, and quickly spreading to the countries that are now called Iran, Iraq and the United Arab Emirates (UAE) (Tengberg, 2012; Abbo et al., 2015). The spread of date palms is believed to have occurred over many centuries along two main routes, from Iraq east to Iran, Pakistan and India and west from Egypt to North Africa (Erskine et al., 2004). However, recent analyses of date palm diversity have suggested that North Africa may have been the origin point of domesticated date palm, with subsequent dissemination east to Asia (Hazzouri et al., 2015). European exploration and colonization have spread date palm throughout much of the world. There is a growing interest in the genetics and genomics of date palm because of its exceptional contributions to the livelihoods of desert farming communities in Asia and North Africa, and also because of its potential as a source for valuable agronomic and stress-tolerance traits. Most of the wild date palm germplasm is already lost, and only a few natural populations are believed to exist (Gros-Balthazard et al., 2017; Wales and Blackman, 2017). Genetic diversity existing in the cultivated date palm is also being lost because of shifts to fewer modern varieties (El-Juhany, 2010). This will increase the vulnerability of date palm to sudden changes in climate, diseases (e.g., Bayoud, and lethal yellowing) and insect pests (e.g., red palm weevil).

Date palms have a juvenile phase of 5-8 years and life spans of over 50 years. Superior cultivars are propagated and distributed by vegetative propagation from offshoots or tissue culture. Hence, though date palm domestication is quite ancient, it is a relatively recent event on a generation time scale, in comparison to annual crops. Modern date palm improvement would be facilitated greatly by a genetic description of the diversity that is present worldwide and the correlation of particular desired traits with specific markers (Morrell et al., 2012). The value of genetic diversity to modern plant breeding is enormous. Such important traits as improved yield, disease and abiotic stress resistance, improved fruit quality, and longer shelf life have been successfully transferred from landraces and wild germplasm to elite cultivated varieties in several crop species. In spite of the many potential benefits of wild Phoenix germplasm, there is no reported effort on the phenotypic evaluation of the wild relatives and hybridization of cultivated date palm with other Phoenix species. Also, an important date palm-specific problem is the difficulty in identifying cultivars until the fruit is produced. Thus, efficient assessment of genetic composition in date palm requires markers that differentiate similar-looking cultivars. Moreover, these markers should be able to identify the genotypes of new varieties with desirable agronomic traits that may have emerged spontaneously in remote (e.g., oasis) locations through sexual reproduction. Studies of genetic diversity and gene flow in date palm can help to devise targeted approaches to conserve germplasm diversity and to investigate evolutionary processes within the genus.

Microsatellites, otherwise known as Simple Sequence Repeats (SSRs), have been used to assess genetic diversity and relatedness of date palm varieties in Algeria (Akkak et al., 2009), Iran (Arabnezhad et al., 2012), Iraq (Jubrael et al., 2005), Tunisia (Hamza et al., 2011, 2012; Zehdi-Azouzi et al., 2015), Saudi Arabia (Al-Abdoulhadi et al., 2011), Sudan (Elshibli and Korpelainen, 2008) Oman (Al-Ruqaishi et al., 2008), Qatar (Ahmed and Al-Qaradawi, 2009; Elmeer and Mattat, 2015), UAE (Chaluvadi et al., 2014), Libya (Racchi et al., 2014), and Morocco (Sedra, 2010). Recent studies have looked into the genetic diversity of worldwide date palm germplasms, using either SSRs (Zehdi-Azouzi et al., 2015), single nucleotide polymorphisms (SNPs) (Hazzouri et al., 2015; Mathew et al., 2015) or comparisons of whole genomes (Hazzouri et al., 2015). Each of these studies investigated different sets of accessions (with some common genotypes) and they used different tools to assess diversity. In general, these broad germplasm investigations identified separate germplasms for Asian and Africa accessions, and some structure within these major subgroups. Further analyses are needed to investigate the complete germplasm collection for this important desert crop and its wild relatives.

The principal aim of this work was to resolve the phylogeny of the genus Phoenix using whole plastome assemblies of all Phoenix species and to analyze a broad distribution of cultivated date palm and related species collected from 13 countries to evaluate the genetic diversity and structure of date palm germplasm. These results were analyzed to understand better the relationship between date palm accessions collected from North Africa and Asia, and to assess the roles of biogeographical history and human activity on the current diversity and distribution of date palm germplasm.

2. Methods 2.1. Chloroplast genome assembly, annotation and phylogenetic analysis

Genomic shotgun sequence data were obtained from 27 Phoenix accessions representing 14 Phoenix species, available in NCBI Genbank as a part of a previous study (Torres et al., 2018) (Supplementary Table 1). The reference plastomes assemblies of P. dactylifera cv. Khalas (GenBank: NC_013991.2) and cv. Aseel (GenBank: FJ212316) were also included in this analysis. This study included three important accessions of P. dactylifera, Khalas, Deglet Noor and Aseel, which are major cultivars in Saudi Arabia, North Africa and Pakistan, respectively. Plastid-homologous sequences were selected from raw Illumina sequence data of Phoenix accessions and then de novo assembled using Velvet and further scaffolding in Geneious 10.1.2, as previously described (Frailey et al., 2018; Vaughn et al., 2014). Plastome assemblies were annotated using the program DOGMA (Wyman et al., 2004). Protein-coding regions, rRNAs, tRNAs, introns, and intergenic regions were all annotated and extracted from DOGMA. A total of 52 gene sequences shared by all Phoenix accessions and one Elaeis guineensis (oil palm) accession were individually aligned using ClustalW (Larkin et al., 2007) in Geneious 10.1.2 (https://www.geneious.com). The individual genes used in the phylogenetic analyses were accD, atpB, atpF, ccsA, infA, ndhA, ndhC, ndhD, ndhF, ndhG, ndhH, ndhI, ndhK, orf188, petA, petB, petD, petG, petL, petN, psaA, psaB, psaC, psaI, psaJ, psbE, psbG, psbH, psbI, psbJ, psbK, psbL, psbN, psbZ, rbcL, rpl14, rpl16, rpl20, rpl22, rpl32, rpl33, rpl36, rpoA, rpoC1, rpoC2, rps2, rps3, rps4, rps8, rps11, rps15, and rps16. Identically aligned regions of each gene sequence were extracted, concatenated, and realigned using ClustalW in Geneious 10.1.2. The E. guineensis plastome was utilized as the outgroup. The Bayesian phylogenetic analysis was done on MrBayes 3.2.1 (Ronquist et al., 2012), as reported in a previous study (Frailey et al., 2018). A separate phylogenetic analysis was also carried out using the alignment of whole plastome assemblies of all the Phoenix species. The whole plastome assemblies were aligned using progressiveMauve aligner, which accurately aligns colinear sequences even if they have undergone large numbers of nucleotide substitutions, indels and rearrangements (Darling et al., 2004, 2010).

2.2. Sampling of date palm accessions for population genetic analysis

The initial date palm collection for our population genetic analyses consisted of DNA samples from 210 cultivated accessions of P. dactylifera and 16 accessions of other Phoenix species from 13 countries. These samples include 112 accessions collected and maintained by the United States Department of Agriculture - National Clonal Germplasm Repository for Citrus & Dates (USDA-NCGRCD) in California, 81 samples collected and maintained at research stations in UAE, 12 accessions collected from Tunisia, 19 accessions collected from India and two accessions collected from Pakistan (Supplementary Table 2). DNAs were isolated from leaf samples using DNeasy plant mini kits as per manufacturer's (Qiagen, USA) instructions. The DNA sample representing each accession was extracted from a single leaf of a single plant of that cultivar. Multiple (commonly three) trees were independently sampled for each accession grown in the California USDA collection (Supplementary Table 2).

2.3. SSR analysis

Nineteen polymorphic nuclear SSR loci (Elmeer et al., 2011, Zhao et al., 2012) (Supplementary Table 3) were amplified with polymerase chain reaction (PCR) performed using a three-primer system with an M13 universal fluorescent-labeled primer (FAM, HEX, NED), as described previously (Chaluvadi et al., 2014). The PCR products were detected and sized on an ABI 3730 sequencer (Applied Biosystems, Foster City, CA). Resultant chromatograms were scored using Soft Genetics Genemarker Version 2.4.

2.4. Population genetic analyses

Of the 226 accessions investigated, 195 accessions (which included 182 cultivated varieties/landraces of P. dactylifera and 13 accessions of other Phoenix species) were chosen for full analysis. The accessions were chosen because they yielded less than 11% null alleles with 19 microsatellite loci, and they also reduced overrepresentation of accessions originating from Iraq and UAE in our initial sampling. Population genetic statistics were calculated using Identity (Wagner and Sefc, 1999), Microsat 2.0 (Minch et al., 1995) and GenAlEx (Peakall and Smouse, 2012). These analyses included characterization of allelic diversity, taxon-specific alleles, Shannon diversity index, expected heterozygosity (gene diversity), observed heterozygosity, and Wright's F_ST (Wright, 1950). We also tested deviations from HardyeWeinberg Equilibrium with GENEPOP version 4.5 (Rousset, 2008).

To partition the genetic variance based on the continent and country of origin, we analyzed molecular variance (AMOVA) (Excoffier et al., 1992). The fixation index F_ST was estimated with the overall dataset within each country, and among all pairs of populations. The statistical significance of F_ST estimates was tested by 1000 random permutations of individuals across populations using GeneAlEx. Nei's genetic distance (Nei and Chakraborty, 1973) and Chord genetic distance (Cavalli-Sforza and Edwards, 1967) were calculated between pairs of genotypes for use in cluster analysis and Principal Coordinates Analysis (PCoA).

The genetic structure of the population was studied using the model-based (Bayesian) clustering method implemented in software package STRUCTURE Ver 2.4.1 (Pritchard et al., 2000). Predefined numbers of populations (k) ranged from 2 to 15, with an initial burn-in period of 50, 000 replicates and 50, 000 Markov Chain Monte Carlo (MCMC) iterations. An accession was assigned to a cluster if the admixture coefficient was >80% (Qi > 0.8) for that group. Accessions with membership probabilities less than 0.8 were assigned to an admixture group (Diez et al., 2015; Stich et al., 2005). Ten independent simulations were run for each K value. We did not use prior information to define the clusters. Because these analyses require codominant alleles and are sensitive to missing data, only 19 microsatellites (each with fewer than 11% missing data) were used. The average K value was calculated from the ten runs and Delta K was calculated by a web-based program, Structure Harvester (Earl, 2012), to identify the number of populations that best reflect the population structure of our samples (Evanno et al., 2005). The STRUCTURE patterns chosen for display (Fig. 5 and Figure S1) were those with the highest statistical support from the ten independent runs. The multiple replicate runs from STRUCTURE were integrated with CLUMPP software (Jakobsson and Rosenberg, 2007).

3. Results 3.1. Plastome assembly and phylogenetic analysis

Full chloroplast genomes were assembled for most Phoenix species. Wherever we could not assemble, the gaps were filled with Ns. The missing regions in the plastomes of P. atlantica_female, Phoenix caespitosa_female, Phoenix canariensis_female, Phoenix paludosa_female, Phoenix pusilla_male, and Phoenix reclinata_female are likely less than 100 bp in total based on comparison to other sequences in this study. The plastome sizes of Phoenix species ranged from 156, 496 bp in P. reclinata_male to 160, 758 bp in Phoenix pauludosa_female. All the Phoenix genomes have 68-70 single copy genes, 18-20 duplicated genes, 31 tRNA genes and four rRNA genes. The whole plastome sequences of all the accessions are included in Supplementary Data. Sizes of all sequenced plastomes and their composition can be found in Supplementary Table 1. There were no major rearrangements in the sequenced plastomes relative to the P. dactylifera assembly.

Phylogenetic trees calculated by the maximum likelihood and Bayesian analyses were congruent in the overall topology, hence only trees calculated by Bayesian analysis using MrBayes are shown (Fig. 1A and B). The phylogenetic tree based on the alignments of only genes (Fig. 1A) has posterior probabilities (PP) of 0.95-1.0 for all nodes except the P. atlantica node. The plastome sequences of male and female accessions of each species were not identical. However, the male and female accessions of each species clustered as pairs in all the Phoenix species, except a few cases (Phoenix acaulis, Phoenix rupicola and Phoenix theophrasti), perhaps indicating a degree of sequence change since their divergence or mislabeling of some cultivars. Fig. 1A shows that the cultivated date palm accessions shared the same clade with male and female accessions of P. sylvestris, Phoenix caespitosa, P. atlantica and male accessions of P. acaulis and P. pusilla. The posterior probability of the P. atlantica branch was only 0.55 while all the other nodes in the phylogenetic tree were from 0.95 to 1.0. The date palm cultivar Aseel appeared at the base of the P. dactylifera clade. Fig. 1B provides the phylogenetic tree based on the whole plastome alignments. The branching structure of this phylogenetic tree is mostly in agreement with the chloroplast gene-based phylogenetic tree (Fig. 1A) with a few differences. The posterior probabilities of all the nodes are 1.0, except the node where P. reclinata (PP = 0.93) and the node where P. theophresti (PP = 0.93) diverged. All the male and female accessions of each species, except with P. acaulis, clustered together.

Based on the phylogenetic analyses of plastome sequences, we separated the phylogenetic trees into three groups. We call Group A the dactylifera group, which includes P. dactylifera, P. sylvestris, Phoenix acualis, P. caespitosa, P. pusilla and P. atlantica. Group B has P. rupicola, P. theophresti, P. canariensis, P. pauludosa, Phoenix roenbereni and P. reclinata. The third group, the most distantly related to P. dactylifera, consists of Phoenix hanceana and Phoenix andamanensis.

3.2. Genetic diversity analysis

An average of 42.1 alleles per locus was detected within the sample data (Table 1). The genetic variation in each accession and each locus, as estimated by the number of alleles, observed heterozygosity (Ho), expected heterozygosity (He) and fixation index (F_ST), are presented in Tables 1 and 2 and Supplementary Table 2. The values of expected heterozygosity ranged from 0.64 for the locus DP172 to 0.94 for the locus JLBDP20. The values of observed heterozygosity were lower than the value of expected heterozygosity at most of the loci (0.14-0.85). Values for the mean allelic diversity (richness) per locus ranged from 1.1 in date palm genotypes AbuMan, Dabbasi, Medjool and Shahil to 1.8 in date palm genotypes Ashrasi2121, BlackSphinx3361, Dayri7490, Halawy7419, Haziz7434, Haziz7435, Sayer7402, Zahidi749, Ashrasi7461, Ashrasi7462, Khadrawy2254, Khadrawy7454 and Sayer7502.

Unique alleles specific to each accession were identified for most of the accessions (Supplementary Table 2). Values for the mean heterozygosity ranged from 0.05 in date palm genotypes AbuMan, Dabbasi, Medjool and Shahil to 0.42 in date palm genotypes Ashrasi7461, Ashrasi7462, Khadrawy2254, Khadrawy7454 and Sayer7502 (Supplementary Table 1). The accessions of Iraqi origin have the highest allelic diversity and the highest number of private alleles (Table 2), but this may be largely because we had more accessions from Iraq (42) than from any other country (ranging from 2 to 20). We did not find any accessions that were 100% identical at all tested loci, even among an independent sampling of accessions with the same cultivar name.

Principal Coordinates Analysis was performed to visualize the relationships between the date palm accessions (Fig. 2A). The first two axes explained 20.4% and 15.4% of the variability, respectively. These two components separate the studied date palm accessions into two loose clusters with the remaining accessions scattered in between. One group contains accessions that were collected predominantly from Asia (UAE, Oman, Iraq, Saudi Arabia, Pakistan, and India), whereas the second group included accessions largely from North Africa and accessions of unknown origin from the USA. The data points scattered in between were mostly composed of accessions from Tunisia, thereby suggesting a particularly diverse and intermixed germplasm pool in this nation. Non-parametric AMOVA to find significant differences between continents produced a pvalue of 0.001. A separate analysis was conducted to include 16 accessions from India. The results were mostly in agreement with Fig. 2A except that the accessions collected from India formed a third cluster that was mostly related to the African germplasm (Fig. 2B).

From the accession genotyping data, we calculated pairwise chord genetic distances between all genotypes to generate a dissimilarity matrix from which a Weighted Neighbor-Joining tree was calculated. P. reclinata was used to root the phylogenetic tree. The relatedness tree based on chord genetic distances depicted in Fig. 3 showed 12 well-resolved clusters. The wild Phoenix species formed two clusters. One cluster consists of P. reclinata, P. sylvestris, P. hanceana and P. canariensis and one P. dactylifera accession from Iran. The other cluster contains Phoenix roebelenii, P. paludosa, Phoenix acaualis and some accessions from Algeria, Egypt, USA, and Iraq. Indian accessions formed a distinct cluster. The neighborjoining tree indicated that date palm accessions are not clustered primarily in a manner related to their country of origin. There are, however, country-specific clusters of accessions originating from India, Iraq, UAE, Tunisia and Morocco.

Most often, accessions with the same name, although collected from different locations, clustered together. However, some accessions with the same name but collected from different locations did not cluster together. For example, multiple accessions of Medjool and Deglet Noor were found in 2 and four different clusters, respectively. Moreover, we were able to identify alleles specific to each cultivar collected from a given location.

We observed few loci/country combinations that were in HWE (Supplementary Table 4), as expected given the history of natural and human selection acting on this crop. The mean F_ST value in the AMOVA analysis based on the average across continents and countries are 0.15 and 0.44 respectively. The hierarchical AMOVA revealed that the majority of total genetic variance (92%) was due to variation among accessions, while 8% was due to variation among countries. An analysis of genetic variation between Asia and Africa showed that 98% of the total genetic variance was due to the variation among accessions within each continent and only 2% of total variance was due to variation among continents (Supplementary Figure 1). We tested the significance of variance components by pairwise population comparisons of F_ST values using a non-parametric permutation approach (Excoffier et al., 1992). Supplementary Table 5 shows a matrix of the results for all the populations. Fifteen of 21 pairwise country comparisons found that the accessions originating from different countries were significantly different (p < 0.001).

3.3. Population structure analysis

Bayesian structure analysis followed by the delta K measure was used to estimate the number of sub-populations (Evanno et al., 2005). We observed the highest delta k peak at k = 2, which is followed by smaller peaks at k = 3, k = 5, k = 7, and k = 12. Fig. 4, derived from ten replicate runs from STRUCTURE integrated by CLUMPP software (Jakobsson and Rosenberg, 2007), shows the calculated accession ancestry at different k values. The samples are depicted by the approximate geography of origin, with African accessions on the left and Asian accessions on the right. At k = 2, each accession is color-coded with either African alleles (blue) or Asian alleles (orange). Fourteen of the 54 accessions originating from Africa have Asian alleles. The majority of these accessions with Asian alleles came from Tunisia. Similarly, 28 Asian accessions out of 99 have African alleles. The majority of these Asian accessions are from the collections originating from Iraq. At k = 3, we observed a distinct group with 17 accessions from Africa and five accessions from Asia.

At k = 12, when the STRUCTURE results are organized by degree of similarity, individual clusters are more clearly defined (Fig. 5). For instance, Iraqi accessions from 5 out of seven groups are admixtures with at least some portion of the genome sharing Iraqspecific alleles (Fig. 5). A new distinct group appeared in UAE accessions. Indian accessions were more like African accessions in all STRUCTURE analyses with k > 2, and Pakistani accessions clustered with Indian accessions (Figs. 4 and 5).

4. Discussion

This study presents an intrageneric phylogeny of Phoenix derived from 52 shared chloroplast genes and whole plastomes of 29 accessions representing 14 Phoenix species. This analysis is a substantial improvement relative to prior studies of the genus phylogeny, which used very few loci and/or few species. The phylogenetic trees based on the shared plastid gene sequences and the whole plastome sequences agree in that the cultivated date palm accessions are most closely related to the male and female accessions of P. sylvestris, Phoenix cespitosa, P. atlantica, and male accessions of P. acaulis, and P. pusilla. Some recent literature considers that P. atlantica, which is a native to Cape Verde islands, may not be a separate species but rather a feral population of P. dactylifera that was naturally vectored to these isolated islands from cultivated date palms in Africa (Gros-Balthazard et al., 2017; Pintaud et al., 2010). Of the other four species in the dactylifera clade (Group 1), three species (P. sylvestris, P. acaulis, P. pusilla) are native to the Indian subcontinent, while P. cespitosa is native to the Arabian Peninsula and Somalia. The plastome-based phylogenetic trees showed marked cytonuclear discordance when compared with the nuclear gene-based (Cyp703) phylogenetic tree reported in Torres et al. (2018). The Cyp703-based phylogenetic tree suggested that P. dactylifera is more closely related to P. theophresti, P. atlantica and P. reclinata than to Phoenix sylverstris. However, both these studies agree that P. sylvestris may not be the direct progenitor of cultivated date palm. It could be that the interbreeding of more than one Phoenix species followed by natural and human selection might have resulted in the origin of P. dactylifera. Given that our chloroplast results only track maternal inheritance, one possible explanation for the different conclusions between our study and that of Torres et al. (2018) is that the female ancestor was closely related to P. sylvestris, while the male ancestor of P. dactylifera may have been a species (perhaps now extinct) that was no more closely related to P. sylvestris than it was to P. acaulis, P. cespitosa or P. pusilla.

Our microsatellite data analysis did not provide any evidence regarding the most closely-related species to domesticated P. dactylifera among the wild Phoenix species. P. canariensis, P. hanceana, P. reclinata and P. sylvestris were all placed in a separate cluster and thus were equidistant from P. dactylifera. This is in contrast to the previous finding that P. sylvestris is the most likely progenitor to the cultivated date palm (Gros-Balthazard et al., 2017). That we did not detect an unusually close relatedness of cultivated date palm with P. sylvestris could be a function of the relatively small dataset employed in our study, compared with the whole genome comparison employed earlier (Gros-Balthazard et al., 2017). One unexpected observation, however, was that Iranian accession Zamordi of P. dactylifera clustered with P. canariensis. We expect that this is caused by a mislabeling of the Iranian material, but a further investigation of this issue is warranted.

The SSR genotyping data in this current study and in several earlier studies (Chaluvadi et al., 2014; Elmeer et al., 2011; Moussouni et al., 2017) have shown that observed heterozygosity was lower than the expected heterozygosity in date palms. Expected heterozygosity increases with an increase in the number of alleles and with an even distribution of alleles. Deeper SNP marker analyses of date palm accessions also observed the higher frequency of expected then observed heterozygosity in data palm and further indicated that long runs of homozygosity (up to 500 kb) could be found within otherwise heterozygous genotypes of date palm (Hazzouri et al., 2015a). This result could be an outcome of occasional inbreeding due to farmer selection for desirable traits (Hazzouri et al., 2015a), but could also be an outcome of mitotic recombination (Rovcanin et al., 2014). These two models can be evaluated when locations of centromeres are determined on the scaffolds, because mitotic recombination is expected to yield homozygosity from the site of recombination (e.g., double strand DNA breakage repair) that extends to the end of the chromosome arm (i.e., the telomere). The fact that a high level of heterozygosity remains in most of these accessions indicates that farmer selection or natural selection can still be acting on heterozygote versus homozygote fitness even in this vegetatively propagated crop.

Though date palms are predominantly propagated vegetatively, most of the accessions used in our study as well as in other studies (Elmeer et al., 2011; Racchi et al., 2014) are highly heterozygous and rich in allelic diversity. Although commercial groves are often exclusively female, purchased pollen can exhibit genotypedependent variability in its effects on fruit size, quality and maturity, otherwise known as 'Metaxenia' (Swingle, 1928; Crawford, 1936; Nixon, 1936). Thus, chance propagation of resultant seed, and clonal propagation appear to have shaped the evolutionary dynamics of date palm even after domestication.

The SSR genotypes classified all date palm germplasm into two major groups with one group predominantly enriched with African accessions and the other group enriched with Asian accessions. This agrees with earlier studies and also suggests that date palm may have been independently domesticated in Asia and North Africa (Hazzouri et al., 2015a; Mathew et al., 2015a). The oasis agrosystem is common for date palms in the deserts of West Asia and North Africa, and tends to generate population divergence because of the often great distances between oases. Date palm domestication along the Persian Gulf has been documented (Beech and Shepherd, 2001; Hazzouri et al., 2015a), but may not have been detected in Africa because of a lack of either early written evidence or archaeobotanical studies (Terral et al., 2012).

It is particularly interesting that the Indian accessions were genetically narrow, well differentiated from all other accessions (except Pakistani accessions) and most similar to Asian accessions at k = 2. However, with the greater differentiation ability at k = 3, Indian and Pakistani accessions showed more alleles that are present predominantly in African accessions than in Asian accessions. A greater similarity of African and Pakistani accessions was also predicted in a previous study (Mathew et al., 2015). These results suggest that India and Pakistan received their date palm germplasm primarily from Africa, and not from the nearer germplasm sources in the Middle East. Whether this reflects specific trade patterns, specific cultural relationships or shared environmental demands of transplanted date palms is not clear. However, these results suggest that investigation of date palms in Somalia and Ethiopia might be particularly informative.

Date palms were first introduced into the New World by the Spanish during the colonial period. The industrial planting of date palm began in the late 1800s, mostly in the low desert areas of California and Arizona. To support this industry, the United States Department of Agriculture imported germplasm beginning in the late 1800s (Johnson, 2010). Although the main cultivars grown in the US are of North African origin, the germplasm collection includes varieties of Asian origin as well as locally developed varieties (Wright, 2016). Most of the US accessions used in our study have records regarding their donors and place of origin, but some do not. Our analyses of the USDA accessions without prior information on origin showed that they are primarily from North Africa. These USDA collections are already seeing a reverse migration, as several Arab countries undertook large-scale programs to increase date palm acreage but found a shortage of suitable offshoots. Primarily because of the absence of major pests and pathogens, California and Arizona have become highly desirable sources of offshoots (Johnson, 2010).

We observed that independent samples of genotypes with the same name, collected from different locations, usually clustered together in phylogenetic analyses. Starting from an offshoot is expected to dramatically decrease the likelihood of genetic divergence between samples of the same cultivar, and the farmer has a secondary check on accession validity because each cultivar has very distinctive fruit traits. However, differences were observed in all accessions with the same name, confirming that mutation is ongoing, especially for highly polymorphic markers like SSRs. In a few cases (e.g., for some Medjool sources), dramatic variation was observed, suggesting unintentional outcrosses. According to Devanand and Chao (2003) and Elhoumaizi et al. (2005), Medjool can be considered an ancient landrace as well as a modern cultivar. The SSR analysis described herein is an efficient and definitive technology for discerning such disparities from the expected genotypic constitution, and thus may become a tool for routine germplasm assessment and verification.

5. Conclusions

Our results show that the plastomes of cultivated date palm accessions are most closely and about equally related to the plastomes of the male and female accessions of P. sylvestris, P. caespitosa and male accessions of P. acualis, and P. pusilla. Future comparisons of Phoenix nuclear and organellar genomes, using our results as a baseline, should be able to identify the level of intercrossing between other Phoenix species and P. dactylifera, with the conflicting potential to both erode natural variation in wild species and to provide new allelic variation for domesticated date palm improvement. The date palm accessions display high genetic diversity and relatively low observed heterozygosity in the analyzed SSR loci. The accessions are genetically structured according to their geographic origin and form two main groups, African and Asian. Most of the US accessions, both with known and with unknown origins, are closely related to African genotypes. Indian and Pakistani date palms appear to have a distant African origin rather than an Asian origin, while Tunisia is unusual for its robust mixture of both African and Asian genotypes. Future studies on date palm germplasm should be targeted on providing information regarding the importance of heterozygosity versus homozygosity in particular genomic regions of the date palm genome. These analyses could indicate mechanism(s) for the origin and possible agronomic value of the homozygous regions that are consciously or unconsciously selected in breeding programs. This study also suggests future direction regarding the sources of alleles related to geographical adaptation and future breeding for improved cultivars.

Conflict of interest

The authors declare no conflict of interests.

Acknowledgments

This research was funded by endowment funds from the University of Georgia Giles Professorship and the Georgia Research Alliance. The authors also thank Jeff Wagner at the Georgia Genomics and Bioinformatics Core for assistance with SSR genotyping.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2018.11.005.