Sequences Evolution and Population Structure of Tunisian Date Palm (Phoenix dactylifera L.) Revealed by Chloroplast DNA Markers

Abstract

Date palm is present among the vital crops of arid and semiarid countries of North Africa and the Middle East. Chloroplast DNA is the best molecule for finding the evolutionary history of plant species. In the present study, cpDNA variation in date palm was estimated using the trnL-trnF intergenic spacer and psbZ-psbC region. The high AT values in both molecular markers may clarify the high proportion of transversions observed in this species. The neutrality tests, expansion parameter estimation (mismatch distribution), and haplotype network patterns proposed that demographic expansion had occurred in recent times. Furthermore, the taxa distribution is not related to geographical origins; neighbor-joining trees are clustered independently either from their geographic origin or from the sex of trees, suggesting a common genetic basis between different cultivars. Statistical analysis of chloroplast germplasm provides a means of assessing cytoplasmic gene flow, which occurred in Tunisian Phoenix dactylifera L. In fact, Nm was important between Tunisia and Eastern Arabic region (Nm = 2.57), which reflects high levels of connectivity between these population pairs. In conclusion, genomic studies prove date palm domestication happened in the Arabian Peninsula and showed an important gene flow with North African palm populations.

1. Introduction

Date palms (Phoenix dactylifera, Arecaceae) are of major economic and ecological importance to the oasis agriculture of arid and semiarid zones. The past distribution area of this fruit crop covers Mauritania in the west to Pakistan in the east and northwestern India [1]. It is moreover existing in sub-Saharan Africa and has been hosted in California, Peru, Australia, and other countries [2]. The presence of this crop in these areas gives it an undeniable ecological role by limiting the progression of steppe areas and the silting of farmland.

P. dactylifera was brought to Tunisia by Phoenicians before the Roman occupation and has played a vital part in establishing oases. In Tunisia, date palm cultivation covers over 46 thousand ha, with a total number of palm trees of approximately 5.4 million [3]. One hectare of land host an average of 120 date palm trees [4]. Further, the Tunisian date palm germplasm is distinct by a remarkable richness, represented by the high number of cultivars (over 250) [5, 6]. Throughout the plant domestication procedure, breeding actions, selection, migration, and admixture have given growth to cultivated populations distinct from the inherited gene pools [7]. Humans have particularly selected qualities associated with production, fruit quality, and fertility [8].

Understanding the population genetics and domestication history of cultivated species is certainly important for the genetic improvement of crops relying on the conservation and usage of the germplasm [9]. In order to contribute to the varietal improvement of date palms and to offer novel perceptions on the influence of geographic origins and human action on the genetic structure of the date palm, this study investigated the diversity of the species using chloroplast DNA.

The small, sensibly constant size, and conservative evolution of chloroplast DNA (cpDNA) make it an ultimate molecule for finding the evolutionary history of plant species. In angiosperm, the use of chloroplast DNA sequences for intra-specific phylogenetic research is now routine [10]. These investigations have been simplified by the abundant number of whole chloroplast genome sequences that are accessible from an extensive variety of angiosperms and the development of universal PCR primers in conserved coding as well as noncoding regions.

Commonly, in terms of its size, organization, and sequence, cpDNA is the most recognized conservatively evolving genome. It has been used for genetic analysis in plants and provides an accessible and well-characterized source of comparative sequence data. Moreover, identifying the footmark of positive selection is an imperative assignment in evolutionary genetic studies [11]. Otherwise, different modes of selection may result in divergent patterns of the nature and extent of genetic variation. The neutral theory affirms that all observable mutations in populations have little or no effect on an organism’s fitness, and their evolutionary dynamics are entirely measured by genetic drift [12, 13]. For such mutations, the evolution continues as equilibrium among the forces of mutation pressure, natural selection, and genetic drift. New molecular techniques based on single and combined sequences data sets have provided a vast understanding of the evolution of flowering plants [14]. Noncoding regions are usually less sensitive to natural selection than coding regions and then may be more beneficial for studying plant evolution.

The chloroplast genome can be classified into three functional categories: (1) protein-coding genes, (2) introns, and (3) intergenic spacers. It has greater phylogenetic potential than nuclear DNA because it is sufficient variable but conserves to be less variable within than between species [15]. The noncoding regions are leading systematic molecular, phylogeographic, and DNA barcoding studies for plants [16, 17]. Chloroplast DNA markers are used for systematic studies of plant species [18, 19, 20, 21, 22, 23] and are particularly used to study phylogeny.

In Paris genus, Song et al. [24] described eight most variable barcodes, to discriminate between the different species, including psbC-trnS-psbZ region. Moreover, the psbC-trnS intergenic spacer was successfully used for phylogeographic study in Lolium species [25]. In 2013, Ballardini et al. confirmed that the psbZ-trnfM (CAU) region could be considered a good basis for the establishment of a DNA barcoding system in Phoenix, and is potentially useful for the identification of the female parent in Phoenix hybrids.

The date palm chloroplast genome is a typical circular double-stranded DNA molecule, and it shares a common quadri-partite structure: a pair of IRs (27.276 bp) separated by the LSC region (86.198 bp) and the SSC regions (17.712 bp) [26]. An important set of primer pairs improved for PCR amplification and covering sequencing in monocotyledons were optimized, which are distributed throughout the whole chloroplast genome, including exons, introns, and intergenic spacers (IGS).

Thus, in the present study, we undertake to find out the level of cpDNA variation in date palm using the trnL-trnF intergenic spacer and psbZ-psbC region and to accurate evolution process, which controls the pattern of polymorphism among the species.

2. Material and methods

2.1 Plant material and DNA extraction

Twenty-four accessions (20 cultivars and four male trees) of Tunisian date palm, listed in Table 1, were used in this study. Each cultivar was represented by one tree. Five varieties accommodated in Tunisian plantations (“Ghars Mettig” and “Tantabecht” from Algeria, “Berhi” and “Khadhraoui” from Iraq, and “Abou Meaan” from the United Arabic Emirates) were used in this survey. Young leaves were frozen until their use for DNA purification. Extraction of the total DNA was determined as stated by Dellaporta et al. [27] protocol. DNA concentration and integrity were checked by 0.8% agarose gel electrophoresis according to Sambrook et al. [28].

Origin	Ecotype	Label	psbC-psbZ		Combined TrnL-trnF and psbC-psbZ
			Length	%GC	Total length	%GC
Tunisia	Lagou	Lg	953	39.7	1353	37.88
	BesserHelou	Bh	951	40.1	1356	38.14
	Gasbi	Gb	967	40.1	1347	38.21
	Boufeggous	Bf	960	40.1	1360	38.21
	Hamra	Hm	971	40.0	1370	38.14
	Tazerzit Safra	Ts	952	39.8	1335	37.91
	Goundi	Gd	962	39.2	1358	38.18
	Menakher	Mk	967	40.0	1350	37.76
	Ammari	Am	969	39.9	1372	38.11
	Deglet Nour	Dn	971	40.0	1374	38.06
	Kentichi	Kt	964	39.3	1372	38.44
	Oum Laghlez	Ol	965	40.1	1359	38.14
	Arichti	Ar	973	40.1	1373	38.14
	Guelb Jemel	Gj	973	39.8	1355	37.91
	Kharroubi	Kb	967	39.7	1382	38.06
Algeria	Tantabecht	Tb	976	40.2	1357	38.18
	Ghars Mettigue	Gm	952	39.9	1352	37.99
Iraq	Berhi	Br	962	40.1	1345	38.11
	Khadhraoui	Kd	971	39.8	1373	37.99
UEA	Abou Meaan	Ab	969	40.2	1364	38.18
	Male trees
Tunisia	Borhane2	B2	967	39.9	1368	38.02
	Borhane3	B3	963	39.9	1348	38.06
	CRPh1	C1	966	40.0	1354	38.21
	CRPh5	C5	962	40.0	1362	38.14

Table 1.

Date-palm accessions studied, their origin, accession numbers, and their variation in length, GC and AT contents of the psbC-psbZ region and the combined sequences.

2.2 Amplification and DNA sequencing

Chloroplastic DNA was sequenced for the trnL-trnF (intergenic spacer) and psbC-psbZ (intergenic spacer + gene). Target regions were amplified using universal primers: (Fw1: 5′GGTTCAAGTCCCTCTATCCC3′; Rv1: 5′ATTTGAACTGGTGACACGAG3′) and (Fw2:5′CAACCTTGGCAAGAACG3′; Rv2: 5′TTGACCAACCATCAGRAGA3′) trnL-trnF spacer and for psbC-psbZ region, respectively. Amplification was carried out for 35 cycles, all comprising of a denaturation phase at 95°C for 1 min, annealing at 50°C for 1 min, and an extension step at 72°C for 2 min.

The total volume of PCR reaction was 25 μL, which contained 25 mM of MgCl₂, 2 mM of dNTP mix, 1.6 mM of each primer and 1 unit of DNA Taq polymerase, and 20 ng of DNA. Agarose-gel electrophoresis (1.5%) was used to check the PCR products.

The purified PCR products for the trnL-trnF and psbC-psbZ regions were sequenced in both strands according to the automated Sanger method [29] using automated sequencer ABI PRISM™ 310 Genetic Analyzer (Applied Biosystems). The process consists of the selective incorporation of chain-terminating dideoxynucleotides by DNA polymerases during in vitro DNA replication. In Sanger sequencing, target DNA is copied multiple times, producing fragments of varying lengths. Fluorescent “chain-terminating” nucleotides mark the ends of the fragments so that the sequence can be determined. It is the most widely used method for detecting single nucleotide variations.

2.3 Sequence analysis

The identity of both sequenced regions was confirmed through a BLASTN search in NCBI database [30]. Nucleotide sequences were aligned by Mega 5.2.2 [31]. Several genetic parameters were determined with DnaSP program [32]. Haplotype diversity (Hd) [33] and genetic diversity (Pi) [34] were calculated to evaluate genetic diversity. The average of nucleotide differences (k), the minimum number of recombination events (Rm), and the average number of nucleotide differences among cultivars were also detected. By means of selective neutrality tests, we checked the hypothesis of the mutation/drift equilibrium for a supposedly neutral polymorphism. Selection neutrality for the detected mutations was tested by both Tajima’s D [35] and Fu and Li’s D* and F* methods [36], using the DNAsp program. Demographic parameters were assessed using the distribution of pairwise sequence differences (mismatch distribution) of Rogers and Harpending [37] and site-frequency spectra of Tajima [35]. Moreover, Fu’s Fs statistics [38] was used to confirm the assumption of population growth and range expansion as revealed by the mismatch distribution of Rogers and Harpending [37]. In addition, we calculated the Harpending’s raggedness index (r) corresponding to an estimate of the fluctuation in the frequency of differences between haplotype pairs [39]. In a complementary way, the R2 statistic [40] was calculated based on the differences between the number of singleton-type mutations and the average of the nucleotide differences. These analyses were executed using coalescent simulations implemented in DnaSP software, with 1000 simulated resampling replicates.

To study the genetic relationships between the studied sequences (haplotypes), we used the reduced median network analysis available in the NETWORK software [41]. The phylogenetic relationships between the studied chloroplast haplotypes were reconstructed using the neighbor-joining (NJ) method [42]. NJ builds a tree from a matrix of pairwise evolutionary distances relating to the set of taxa being studied. Gene flow (Nm) was estimated with the mean number of migrants per generation among populations. N_m values were calculated with 1000 data permutations using the software DnaSP v 5.10.01. With reference to the standard for gene flow, we described genetic flow as low for Nm < 1, high for 1 < Nm < 4, and very high for Nm > 4 [43]. In gene flow analysis, we consider three different geographic regions of P. dactylifera L. Tunisia, Algeria, and Eastern Arabic (Iraq and UAE).

3. Results

3.1 Sequences variation in the psbC-psbZ region

The DNA sequencing of the generated bands has been successfully performed and the blast search allowed confirming the identity of the sequences. DNA sequence varied from 951 bp for “Besser Helou” cultivar to 967 bp for “Gasbi” cultivar (Table 1) with an average of 964.5 pb length. The nucleotide composition frequencies are 0.274 (A), 0.325 (T), 0.204 (C), and 0.196 (G). In addition, the GC content of the amplified sequences varied from 39.7% to 40.7%, and the AT content from 59.3% to 60.3% (Table 1).

The polymorphism pattern of the psbC-psbZ sequence in P. dactylifera L. reveals a high level of mutation with 55 polymorphic sites (Table 2). We observed 33 singleton variable sites and 22 parsimony-informative ones. Moreover, a low transitional/transversional ratio (ti/tv = 0.68) occurs in the psbC-psbZ sequence. In addition, 22 haplotypes were detected among 24 cultivars analyzed, yielding a haplotype diversity of 0.989 (Table 2). Sequence variation observed between cultivar groups and nucleotide diversity was estimated (Pi = 0.00862) (Figure 1a). The average number of nucleotide differences “k” was estimated to be 8.264. Furthermore, seven regions of conserved DNA were detected using DnaSP program yielding a value of sequence conservation C = 0.913. Mismatch distribution within the species was unimodal for the psbC-psbZ marker (Figure 2a), suggesting demographic expansion had occurred in recent times.

	psbC-psbZ sequence	Combined sequence
Number of accessions (N)	24	24
Lengths (pb) (Average)	964	1360
GC (%) (Average)	40	37.38
AT(%) (Average)	60	62.66
Number of polymorphic sites (S)	55	69
Parsimony informative sites (Pis)	22	31
Number of haplotypes (H)	22	24
Nucleotide diversity (Pi)	0.00862	0.00861
Haplotype diversity (Hd)	0.989	1
Average of pairwise nucleotide differences (k)	8.264	11.471
Tajima’s D	−1.96519; (0.10 > P > 0.05)	−1.74847; (0.10 > P > 0.05)
Fu and Li’s D*	−2.345430; (0.10 > P > 0.05)	−1.97901; (0.10 > P > 0.05)
Fu and Li’s F*	−2.61261; (0.10 > P > 0.05)	−2.23846; (0.10 > P > 0.05)
Fu’s Fs	−12.927; (0.10 > P > 0.05)	−15.474; (0.10 > P > 0.05)
R2 (Ramos-Onsins and Rozas’s test)	0.000; P < 0.05	0.000; P < 0.05
Raggedness (r)	0.000; P < 0.05	0.000; P < 0.05

Table 2.

Sequences polymorphism and divergence within date palm cultivars.

Figure 1.

Variability of nucleotide diversity (Pi) and segregating sites (S) for the psbC-psbZ chloroplastic DNA (a) and the combined sequences (b).

Figure 2.

Mismatch distribution of chloroplastic DNA sequences of the date palm cultivars based on pairwise nucleotide differences in the psbC-psbZ region (a) and the combined trnL-trnF spacer and psbC-psbZ (b). In addition, graphical representation of the site frequency spectrum of chloroplast DNA sequences. Solid lines in the site-frequency spectra indicate the expected distributions under neutrality and at equilibrium.

A phylogenetic tree was constructed using NJ methods (Figure 3a). The evolutionary distances were estimated by the maximum composite likelihood method [44]. Haplotypes of the psbC-psbZ regions made it possible to make moderately sustained phylogenetic trees (Bootstrap values ≤74%). This dendrogram supported the varieties’ organization into two main clusters, which were subdivided into different subclusters. Clearly, the obtained clustering is made independently from the geographical origin of cultivars since the foreign varieties, newly introduced in the Tunisian growing areas, did not remarkably separate from the autochthonous ones. In addition, the male and female trees did not diverge from each other.

Figure 3.

Summary of phylogenetic analysis within Phoenix dactylifera based on psbC-psbZ sequences: (a) neighbor-joining tree of 24 Tunisian date-palm cultivars. (b) Median-joining network of the haplotypes inferred from the analyzed sequences. Nodes are proportional to haplotype frequencies and branch length is proportional to the number of mutations. The circles represent haplotypes; circle diameters are proportional to the frequencies. Black bands correspond to mutational steps.

In a complementary way, a genetic network based on the psbC-psbZ sequences was reconstructed (Figure 3b). All haplotypes are connected to the most frequent haplotype H1. It represented the ancestor sequence of the other ecotypes, which mutate during evolution. Several putative haplotypes, corresponding to intermediate evolutionary steps, were detected in our network (mv in Figure 3b). The haplotype patterns (star-like pattern) reflect a signature of a recent expansion in Tunisian date palm.

3.2 The combined region trnL-trnF and psbC-psbZ spacers

The combined sequence varied from 1347 bp for “Gasbi” cultivar to 1382 bp for “Kharroubi” cultivar (Table 1). Current sequences do not contain insertions or deletions (indels). These sequences revealed 69 polymorphic sites and defined 24 haplotypes. Among the 69 variable sites, 31 were parsimoniously informative and 38 were singletons sites. The means of the haplotype (Hd) and nucleotide (Pi) diversity (Figure 1b) were higher than for trnL-trnF spacer [16] and the psbC-psbZ region taken separately. In fact, these values were 1 ± 0.014 and 0.00861 ± 0.001, respectively. The average of pairwise nucleotide differences (k) was equal to 11.471 (Table 2). In addition, the GC content of the combined sequences varied from 37.88% to 38.44% (Table 1). The polymorphism pattern in P. dactylifera L. revealed an elevated level of mutation points with a low transitional/transversional ratio (ti/tv = 0.63).

Mismatch distributions within the species were generally unimodal for markers, suggesting population expansions (or past selective sweeps) (Figure 2b). Numerous significant negative values of Fs and Tajima’s D in trnL-trnF [16], psbC-psbZ spacers, and the combined region (Table 2) rule out the hypothesis of a constant population and propose either selective sweeps or past demographic expansion (Table 2). In addition, low values of R2 and Harpending’s raggedness index (r) will be characteristic of a recent expansion.

For gene flow estimation, N_m was high (N_m > 1.0) between Tunisia and Algeria (N_m = 1.8) and Tunisia/Eastern Arabic accessions (N_m = 2.57) indicating that few differentiations could be established among them. On the other hand, Nm was weak between Algeria and Eastern Arabic (N_m = 0.66), this reflects low levels of connectivity between these population pairs.

Genetic relationships among date palm cultivars are investigated, also, using the variation observed in the combined (trnL-trnF and psbC-psbZ) sequences. As illustrated by Figure 4, the phylogeographical patterns were very similar between the psbC-psbZ gene and the combined sequences. The NJ dendrogram supported the variety organization into three monophyletic branches and one cluster. The first branch is composed of “Lagou” (Lg), the second one CRPh5 (C5) Tunisian male tree, while the third is made of Iraquian variety “Berhi” (Br).

Figure 4.

Neighbor-joining tree based on trnL-trnF spacer and psbC-psbZ gene sequences.

The remaining cultivars are ranged in the unique cluster, which is divided into four subgroups. In fact, genotypes’ clustering is independent of the sex of trees, and it is not structured according to geographical origin. Indeed, with the exception of the foreign cultivar “Berhi” (Br) divergence, the other foreign cultivars did not significantly diverge from the Tunisian plantations (Figure 4).

4. Discussion

In the present survey, we tested the reliability of the noncoding chloroplast markers to identify date palm cultivars. Therewith, we sought an indication of the level of genetic variation and genetic classification within the date palm cultivars grown in Tunisia.

In this object, we represented the evolution of the psbC-psbZ gene and the combined sequences (trnL-trnF spacer/psbC-psbZ) to reveal polymorphisms in P. dactylifera L. The high AT values in both molecular markers may clarify the high proportion of transversions (ti/tv = 0.67 and ti/tv = 0.63) observed in this species. This result corroborates the discovery in angiosperm chloroplast noncoding regions with a ratio ti/tv not exceeding one for any of the examined taxa [45, 46, 47, 48, 49, 50]. Base content explains the relatively high proportion of transversions [45].

The use of coalescent theory [51, 52] for the P. dactylifera L. sequences has permitted inferences about past and present population size. It has been confirmed that demographic expansions conduct to star-shaped genealogies [53], an excess of rare mutations [54], and a unimodal mismatch distribution [37]. Haplotype network patterns, the neutrality tests, and the expansion parameter estimations (mismatch distribution) proposed that demographic expansion had occurred in recent times as indicated by the excess of mutations type singletons in the P. dactylifera sample.

Furthermore, there was no geographical structure in the relationships among the haplotypes; neighbor-joining trees are clustered independently either from their geographic origin or from the sex of trees, indicating a common genetic basis between different cultivars. A low genetic structure is usually related to natural habitat change or to human activities that increase gene flow between populations [55]. Human impact on these regions may be the reason for these findings. Actuality, in Tunisian localities, cultivars are manipulated by farmers after continuous selection, cloning, and exchange of varieties. Similar results have been observed using other molecular markers, where P. dactylifera ecotypes clustering made independently from either the tree’s sex or their geographic origin in spite of their great distribution [16, 56, 57, 58, 59, 60]. This result is contrasted with Zehdi-Azzouzi et al. (2015), where NJ classification was clearly coherent with a geographical structuring into two clusters: Eastern pool (Djibouti, Oman, Iraq, the UAE) and the Western pool (North Africa accessions). Also, using three chloroplast regions of rbcL, matK, and trnH-psbA, Iranian date palm cultivars were separated into clades corresponding to their geographical distribution [61].

Statistic analysis of chloroplast germplasm provides a means of assessing cytoplasmic gene flow, which occurred in Tunisian P. dactylifera L. The pattern observed confirmed that gene flow is designed to be a significant evolutionary factor in date palm history, allowing a variability of the genetic diversity at different scales (local and foreign). In fact, N_m was important between Tunisia and Eastern Arabic region (N_m = 2.57), which reflects high levels of connectivity between these population pairs. Zehdi-Azzouzi et al. [9] proved the movement of gene flows between Eastern and Western origins, mostly from east to west, following a human-mediated diffusion of the species.

In addition, genomic studies prove date palm domestication has happened in the Arabian Peninsula and that the North African population has mixed ancestry with components from Middle Eastern P. dactylifera [62]. In this subject, Gros-Balthazard et al. [63] noted that a first domestication event was shown to be improbable, given the very similar form of seeds in cultivars from Africa and the Middle East. In fact, travelers and pilgrims brought many cultivars of date palms grown in Tunisian oases, particularly from the east [64]. The origins of domestication, the direction of germplasm flows, and the breeding history within the traditional cultivation are sincerely the reason for the observed exchange.

In addition, the mean number of migrants per generation between Tunisia and Algeria (N_m) was important; this could be the result of the relative proximity of geographical sites. In fact, Deglet Nour cv. was introduced into Tunisia plantations four centuries ago from western Algeria [64, 65], and became the most valuable cultivar. After this period, the human migrations between Algeria and the south of Tunisia were perturbed by the French occupation of Algeria [64]. Despite that, other works prove that in Algeria, chloroplast diversity presents about 70% of the eastern Arabic chloroplast. In Tunisia and Morocco, this proportion only ranges from 11–42%, although Algerian nuclear diversity is similar to that obtained in Tunisia [66].

In conclusion, environmental change, over-exploitation of water reserves, and an increasing tendency toward Deglet Nour cv. Monoculture consists of some of the difficulties of date palm culture in Tunisia. Those are significant issues, which need in-depth research. Considering these, it is decisive to comprehend the genetic diversity and variation among and within accessions, and extensive programs must be prepared, to include collection, estimation, and preservation of plant genetic resources. Determining the gene flow and the demographic expansion of this population should be possible with our two molecular markers, nevertheless, the exploitation of additional nuclear markers clearly provides more accurate results. Actually, with the availability of the date palm nuclear and chloroplast genome sequences, more molecular markers are accessible to explore systematic, cultivar relatedness, and genetic map structure. Furthermore, a webpage regrouping all information about P. dactylifera germplasm conservation and utilization has an obligation to be planned in collaboration with all producer countries in order to engender bases for inhibiting its genetic erosion.

Conflict of interest

The authors declare that they have no conflict of interest.

Abbreviations

cpDNA	chloroplast DNA
dNTP	deoxynucleoside triphosphate
Hd	Haplotype diversity
IR	inverted repeats
IGS	Intergenic Spacers
k	average of nucleotide differences
LSC	large-single-copy region
MgCl2	magnesium chloride
Nm	Gene flow
PCR	polymerase chain reaction
Pi	genetic diversity
Rm	minimum number of recombination events
SSC	small-single-copy regions
ti/tv	transition/ transversion