Morphological characterization and genetic diversity analysis of Tunisian durum wheat (Triticum turgidum var. durum) accessions

Background Tunisia is considered a secondary center of diversification of durum wheat and has a large number of abandoned old local landraces. An accurate investigation and characterization of the morphological and genetic features of these landraces would allow their rehabilitation and utilization in wheat breeding programs. Here, we investigated a diverse collection of 304 local accessions of durum wheat collected from five regions and three climate stages of central and southern Tunisia. Results Durum wheat accessions were morphologically characterized using 12 spike- and grain-related traits. A mean Shannon-Weaver index (H′) of 0.80 was obtained, indicating high level of polymorphism among accessions. Based on these traits, 11 local landraces including Mahmoudi, Azizi, Jneh Khotifa, Mekki, Biskri, Taganrog, Biada, Badri, Richi, Roussia and Souri were identified. Spike length (H′ = 0.98), spike shape (H′ = 0.86), grain size (H′ = 0.94), grain shape (H′ = 0.87) and grain color (H′ = 0.86) were the most polymorphic morphological traits. The genetic diversity of these accessions was assessed using 10 simple sequence repeat (SSR) markers, with a polymorphic information content (PIC) of 0.69. Levels of genetic diversity were generally high (I = 0.62; He = 0.35). In addition, population structure analysis revealed 11 genetic groups, which were significantly correlated with the morphological characterization. Analysis of molecular variance (AMOVA) showed high genetic variation within regions (81%) and within genetic groups (41%), reflecting a considerable amount of admixture between landraces. The moderate (19%) and high (59%) levels of genetic variation detected among regions and among genetic groups, respectively, highlighted the selection practices of farmers. Furthermore, Mahmoudi accessions showed significant variation in spike density between central Tunisia (compact spikes) and southern Tunisia (loose spikes with open glume), may indicate an adaptation to high temperature in the south. Conclusion Overall, this study demonstrates the genetic richness of local durum wheat germplasm for better in situ and ex situ conservation and for the subsequent use of these accessions in wheat breeding programs. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-021-00958-3.


Background
Durum wheat (Triticum turgidum var. durum Desf.) is a tetraploid species (2n = 4x = 28, AABB) that originated and domesticated in the Fertile Crescent and spread within the Mediterranean region through different dispersal [1,2], reaching the Iberian Peninsula through Northern Africa around 7000 BC [3]. Since then, durum wheat has gained commercial importance. Today, durum wheat is cultivated worldwide, especially in the Mediterranean Basin, which is considered as the center of diversification and production of durum wheat [4,5]. The Mediterranean Basin is characterized by highly variable environments, ranging from warm and dry to cool and wet climates [6]. Durum wheat accessions collected from the Mediterranean region exhibit higher genetic diversity than those collected from other regions of the world [7].
Within the Mediterranean region, Tunisia is one of the main centers of diversity of durum wheat [8,9]. Old Tunisian durum wheat cultivars are known by their high level of genetic diversity and their specific adaptability to North African drylands [10]. Despite their notable genetic diversity, Tunisian landraces have been progressively abandoned since the first decade of the twentieth century and replaced by improved, high-yielding and genetically uniform semi-dwarf cultivars (known as "modern varieties") developed through international breeding programs [11,12]. This has led to a significant reduction in the genetic diversity of local durum wheat [13,14]. Nonetheless, the genetic diversity of durum wheat could be preserved by: (1) characterizing the remaining durum wheat landraces; (2) re-introducing these landraces into breeding programs; and (3) protecting these landraces through effective conservation strategies. Therefore, the genetic and morpho-phenological characterization of landraces, which are either sparsely cultivated under the current farming system or stored in gene banks, would allow the identification of unexplored sources of genetic diversity that may be important for adaptation to several biotic and abiotic stresses [7,15,16]. The availability of landraces for breeding programs may also have particular relevance for breeding cultivars suitable for suboptimal and marginal environments such as the Mediterranean Basin, where durum wheat and other crop species are largely cultivated under unstable and limited water conditions, which cause considerable yield fluctuations [17,18].
Previously, the agro-morphological evaluation of Tunisian durum wheat accessions using quantitative and qualitative spike-related traits, mostly concerning the grains, revealed high morphological diversity within the Tunisian durum wheat landraces [19,20], and more than 35 durum wheat landraces were recorded [13]. However, few studies have been conducted to analyze the morphological and genetic features of durum wheat simultaneously. Moreover, the correlation between genetic population structure and morphological aspects of durum wheat has not been investigated to date. Previously, analysis of the level of genetic diversity in Tunisian durum wheat germplasm using amplified fragment length polymorphism (AFLP) and simple sequence repeat (SSR) markers revealed an important polymorphism within cultivars [10]. More recently, investigation of the genetic diversity and population structure of 196 durum wheat landrace accessions (including Tunisian and North African accessions) using diversity array technology sequencing (DArTseq)-based markers showed that genetic variation was higher among landraces than within landraces, and the Tunisian and North African landraces showed remarkable genetic similarity [21]. Furthermore, Slim et al. [22] evaluated the genetic structure of Tunisian durum wheat germplasm, and suggested the existence of five subpopulations with a strong genetic stratification from the north to the south of Tunisia, probably due to the prevalence of modern cultivars in the north. By tracing the history of cultivation, Tunisian durum wheat germplasm collections have been divided into three distinct categories: traditional varieties or old landraces, old cultivars (cultivated up to 1970s) and modern cultivars (cultivated up to present) [10,13,22]. Since traditional local landraces have been derived either from artificial selection of traditional farming systems or from natural adaptation to adverse growing conditions, these landraces might harbor key traits for breeding programs.
Taking into account the value of traditional Tunisian durum wheat landraces, we aimed to: (i) evaluate the genetic diversity and population structure of 304 Tunisian durum wheat accessions collected from central and southern Tunisia using SSR markers; (ii) study the phenotypic diversity of these accessions, based on the morphological characterization of spike-and grainrelated traits; and (iii) analyze the relationship between genetic and phenotypic variation.
The 304 durum wheat accessions investigated in this study were grouped into 11 landraces, namely Azizi, Jneh Khotifa, Taganrog, Mekki, Richi, Souri, Roussia, Badri, Biskri, Biada and Mahmoudi, recorded in the catalog of durum wheat landraces cultivated in Tunisia [13]. These landraces were characterized by 12 specific morphological traits, based on the International Plant Genetic Resources Institute (IPGRI) [23] and International Union for the Protection of New Varieties of Plants (UPOV) [24] (Table S1, Table S2). All 12 spike and grain characteristics were almost homogeneous among accessions of the same landrace. This was supported by the Shannon-Weaver index (H′), which was relatively low for each landrace, ranging from 0.00 (Badri and Jneh Khotifa) to 0.23 (Richi), with an overall mean of 0.11 (Table S3). For instance, Mahmoudi accessions had particularly large spikes with sub-pyramidal shape, very long awns and large grains, whereas spikes of Azizi accessions were rectangular and very flat. Biskri accessions had fusiform and large spikes. The spike color, length and shape varied among the studied accessions from dark to light and from short to long spikes. For example, Badri spikes were very short and thick with a greyish color, whereas Biada spikes and awns were very light (white) in color. Souri and Roussia were both characterized by tight, red-colored spikes with a distinct spike shape, i.e., either rectangular (Souri) or cylindrical (Roussia). Souri and Roussia landraces were also characterized by a distinct orange colored grain. Interestingly, Richi accessions showed a unique feathery spike, while Mekki accessions were characterized by short and dense spikes with parallel edges. Finally, Taganrog accessions were characterized by white colored spikes with black stains, while Jneh Khotifa accessions showed very dark (black to purple), long and dense spikes and awns.

Principal component analysis (PCA)
PCA of 12 spike and grain morphological traits of 304 durum wheat accessions showed that PC1 and PC2 axes accounted for 25.73 and 22.34% of the total genetic variation in these traits, respectively (Fig. 1). PC1 was mostly associated with SS, SL, number of spikelet per spike (NS), grain color (GC) and awn length (AL), whereas PC2 was mainly associated with GSp, GSz and grain number per spikelet (GN) (Fig. 1a). The colorcoding of accessions in the two-dimensional PCA plot (PC1 vs. PC2) showed a good correspondence between the morphological trait-based grouping and landrace denomination (Fig. 1b), and accessions belonging to the same landrace were included in the same PCA subgroup. Biskri, Jneh Khotifa and Taganrog accessions grouped together, showed positive correlation with both PC1 and PC2 and shared similar spike characteristics, such as SL (mostly long spikes), high GN (> 3), black awn color (AC) and AL longer than the spike. Azizi accessions were grouped into a distinct subgroup, mainly characterized by rectangular medium-sized spikes with a tan color. Mahmoudi accessions also formed a distinct subgroup, mainly characterized by unique pyramid-shaped spikes. Accessions of Souri and Roussia formed almost a single subgroup characterized by red-colored loose and long spikes as well as red colored glumes and awns. Landraces Badri and Mekki formed distinct subgroups negatively correlated to PC1 and PC2, and both subgroups were mainly characterized by short spikes with a low to intermediate GN. Biada and Richi accessions were grouped mainly in the center of the plot and were particularly characterized by white-colored spikes, glumes and awns (Table S2). Overall, PC1 and PC2 could separate all landraces, based on 12 spike-and grain-related morphological traits; the only exceptions were the groups of Roussia and Souri landraces and Biskri, Jneh Khotifa and Taganrog landraces, which could not be distinguished based on SL and SC. Thus, additional morphological traits, such as glume form, were considered to classify the latter landraces into distinct subgroups (Table S2).
Genetic diversity and population structure of Tunisian durum wheat accessions SSR polymorphism Ten SSR markers were used in this study to analyze the genetic diversity and population structure of Tunisian durum wheat accessions. These SSR markers were mapped onto eight different chromosomes and therefore were considered largely independent ( Table 2, Table S4). The percentage of missing data was low (< 10%) for each locus. All 10 SSR markers amplified a total of 99 alleles and from 302 accessions, 188 multilocus genotypes (MLGs) were identified. The accumulation curve ( Figure  S1), showed that these SSR markers were able to reach the maximal range of differentiation among the MLGs.

Analysis of population structure and relationship with morphological characterization
We analyzed the population structure on 188 MLGs. The maximum likelihood (LnP(K)) and delta K (ΔK) methods indicated that the most likely number of genetic groups (K) was 11 ( Fig. 2a, b). The estimated genetic group membership coefficient of each accession at K = 11 is shown in the population structure plot (Fig. 2c).
Overall, each genetic group corresponded to a single landrace. The genetic groups G2, G3, G4, G5, G7, G9, G10 and G11 corresponded to Jneh Khotifa, Taganrog, Mekki, Richi, Badri, Beskri, Biada and Mahmoudi landraces, respectively. Moreover, a significant correlation was detected between the genetic distance matrix and morphological distance matrix (P = 0.01; R xy = 0.435). However, a discrepancy between the genetic grouping and the morphological characterization was observed for Azizi, Souri and Roussia landraces; Azizi landraces were grouped by STRUCTURE into two different genetic groups G1 and G8, while Souri and Roussia landraces were grouped together into one genetic group (G6), despite their distinct morphological characteristics.

Analysis of diversity indices and molecular variance
The 11 groups identified by STRUCTURE analysis presented different levels of genetic diversity (Table 3). , based on 10 replicates per K, for K ranging from 1 to 20, with a burn-in period of 100,000 and Monte Carlo Markov Chain replicates of 100,000 iterations; (c) Membership coefficient bar plot displaying population structure at K = 11 for 302 Tunisian durum wheat accessions genotyped with 10 SSR markers inferred from STRUCTURE [26]. Each MLG is represented by a vertical line and they are ordered by color-coded genetic group (G1 to G11). For each genetic group, corresponding durum wheat landrace is mentioned. * Azizi landrace was divided into two genetic group G1 and G8 Table 3 Diversity indexes of 302 Tunisian durum wheat accessions sorted by genetic groups as defined by STRUCTURE [26], by regions and by climate stages : Frequence (0.7-1). b : a DA is a rare allele with a frequence > 70% for a population or region and < 30% for the others while G2 and G7 contained no private alleles (PA = 0) (Table S5). Groups G10 and G4 contained two diagnostic alleles (DA) each, while G3, G5 and G7 contained one DA each, with frequency > 70%. The fixation index (Fis ranged from 0.698 (G4) to 1.0 (G7), where observed heterozygosity (Ho) was 0.100 and null, respectively. Furthermore, analysis of variance (ANOVA) showed that 59% of the total genetic diversity was observed between distinct genetic groups, while 41% of the genetic diversity was explained by differences within each group (Table 4).
Diversity analysis of Tunisian durum wheat accessions based on regions and climate stages Analysis of morphological diversity among different regions and climate stages The Shannon-Weaver index (H′) of 12 spike and grain related traits was compared among five regions (Sousse, Mahdia, Kairouan, Gabes and Medenine) and three different climate stages (low semi-arid, high-arid and mid-arid) ( Table 1) . Morphological traits were also variable from one climate stage to another. Values of SL and glume color were the highest in high-arid climate (0.48 and 0.96, respectively) and lowest in low semi-arid climate (0.0 and 0.41, respectively). By contrast, AC was the lowest in higharid climate (0.12) and the highest in mid-arid climate (0.71). However, no variation was observed among regions for GC and among climates for GN.
In addition, a dominant phenotypic class of some morphological traits was observed among regions (within more than 70% of accessions), except Sousse, which did not show any variation in morphological traits (Table S9). In Gabes, the majority of accessions showed long spikes (SL > 9 cm; 84%), with light color (92%) and cylindrical shape (79%), awns shorter than the spike (84%), moderately elongated grain shape (82%), small grains (GSz < 0.3 cm) (82%) and an intermediate number of grains per Admixed genotypes were excluded from the analysis spikelet (GN = 2-3; 79%), whereas accessions with medium length spikes (SL: 6-9 cm) were dominant in Medenine (73%). In Mahdia, the majority of accessions showed spikes with AL equal to the spike (72%) and small GSz (78%). However, most of the accessions in Kairouan had spikes with AL longer than the spike (72%). Among different climate stages, the mid-arid was dominated by accessions with small grains (GSz < 0.3 cm; 72%), whereas the high-arid climate stage was rich in accessions with dark colored spikes (72%) and black awns (96%). No particular phenotypic class was observed within the low semi-arid climate (Table S9).

Analysis of genetic diversity among different regions and climate stages
The results of ANOVA showed that 19 and 10% of the total genetic diversity was observed among regions and among climate stages, respectively, while 81 and 90% of the genetic variability was explained by differences within regions and within climate stages, respectively (Table 4). Genetic diversity among regions showed Ne ranging from 1.366 (Sousse) to 3.031 (Gabes) ( Table 3). Overall and among all investigated regions, Sousse showed the lowest genetic diversity indexes, while Gabes showed the highest genetic diversity indexes; the number of MLGs was the highest in Gabes (31) and lowest in Sousse and Medenine (7), and the Shannon's information index was also the highest in Gabes (H′ = 1.296) and lowest in Sousse (H′ = 0.305). Moreover, the percentage of polymorphic loci (P) was 100% for all regions, except Sousse   (50%). Moreover, the number of private alleles was also the highest in Gabes (PA = 17) and lowest in Sousse and Medenine (PA = 1). The value of Fis was greater than 0.800 in each region, except Sousse (Fis = 0.691). Interestingly, the DA number and heterozygosity index were the highest in Sousse. In fact, three diagnostic alleles (frequency > 70%; Ho = 0.100) were registered in Sousse, whereas only one such allele was identified in Gabes. Analysis of SSR data obtained from different climate stages revealed that the mid-arid climate was the most outstanding, with the highest number of effective alleles (Ne = 3.174), the highest Shannon's information index (I = 1.318) and the highest number of private alleles (PA = 19). By contrast, the high-arid climate stage showed the lowest number of effective alleles (Ne = 2.707), the lowest Shannon's information index (I = 1.050) and the lowest number of private alleles (PA = 2). However, the value of Fis was similar (> 0.800) among all studied climate stages (Table 3).

Correlation between genetic distance and geographic distance
The Mantel test showed a significant correlation (P = 0.010; R xy = 0.286) between genetic and geographic distances among durum wheat accessions, suggesting that accessions growing in close geographical proximity were genetically related. Azizi and Mahmoudi landraces showed the most widespread geographical distribution in central and southern Tunisia, except Sousse, and all climate stages. However, Azizi was more frequent in Gabes (25 accessions out of 38), while Mahmoudi was mostly found in Medenine (13 accessions out of 22) and Mahdia (11 accessions out of 27) (Fig. 4). In addition, all G5 genotypes, corresponding to the Richi landrace, were found in Sousse; all G7 and G2 genotypes, corresponding to Badri and Jneh Khotifa landraces, respectively, were found in Kairouan; and the landrace Taganrog, representative of the genetic group G3, was exclusively found in Mahdia.
Furthermore, we compared morphological traits between Azizi and Mahmoudi accessions collected from central and southern Tunisia. None of the traits, showed significant differences (P > 0.05), except for spike density (SD) which showed significant differences within Mahmoudi (P = 0.00). Mahmoudi accessions collected from central Tunisia had compact spikes (SD = 7), whereas those collected from southern Tunisia were characterized by loose spikes (SD = 5) ( Table 5).

Genetic and morphological diversity within the Tunisian durum wheat germplasm
In the present study, we investigated the genetic diversity of 302 Tunisian durum wheat accessions using 10 SSR markers, which enabled maximal differentiation among MLGs, suggesting that these markers have a good resolution power. Overall, the studied collection was characterized by high genetic diversity (overall Na = 9.9; PIC = 0.690; He = 0.346). Similar level of polymorphism (Na = 8; PIC = 0.68) was previously reported using 15 SSR markers in a Tunisian durum wheat collection composed of 7 modern cultivars and 27 old cultivars [10]. More recently, Slim et al. [22] reported genetic diversity indexes (PIC = 0.57; He: 0.28-0.82; Na: 2-13) of Tunisian durum wheat germplasm, consisting of 41 traditional varieties and 13 cultivars, using 16 SSR markers. A higher level of polymorphism (Na = 10; He = 0.71) was reported in a wider geographical collection of 172 durum wheat landraces (collected from 21 Mediterranean countries) and 20 modern cultivars genotyped by 44 SSR markers [27]. However, lower genetic diversity was observed in 33 Anatolian, 136 south Italian and 40 North-West African durum wheat landraces using 14, 44 and 29 SSR markers, respectively [7,15,28]. Robbana et al. [21] also reported low genetic diversity (PIC = 0.1; He = 0.25) in 196 Tunisian durum wheat accessions; this was possibly due to (i) the use of bi-allelic DArTseq markers, which are less informative than the multiallelic SSR markers, and (ii) to limited number of landraces (5). This variability between studies suggests that the ability to capture the maximum genetic diversity depends on the type of markers, number of landraces, their origin and geographical distribution. In this study, the level of phenotypic diversity detected on the basis of 12 morphological traits was consistent with that of genetic diversity, with a Shannon-Weaver index (H′) of 0.80. The morphological diversity observed in this study was higher than that described previously for 930 Tunisian durum wheat accessions (H′ = 0.53), collected from fewer sites in southern Tunisia, using 22 qualitative and 3 quantitative traits [19]. Lower phenotypic diversity was also observed for Moroccan durum wheat populations composed of 101 landraces (H′ = 0.62) [29] and 59 traditional durum wheat accessions (H ′ = 0.78) [30] using nine agro-morphological traits. Ethiopian durum wheat populations composed of 32 landraces showed an H′ value of 0.71 using eight qualitative traits [31], while Oman populations composed of 161 accessions showed H′ value of 0.52 and 0.66 using 15 qualitative and 17 quantitative traits, respectively [32].

Population structure, network analysis and relationships with morphological characterization
In this study, the Tunisian durum wheat germplasm collection was genetically structured into 11 groups, which were consistent with morphological characteristics. Indeed, 8 out of 11 landraces corresponded to distinct genetic groups. This result highlights the effectiveness of SSR markers in distinguishing wheat varieties; thus, these markers continued to be widely exploited for DNA fingerprinting in plants [35]. DArTseq markers also grouped the Tunisian old durum wheat accessions according to the corresponding landrace nomenclature [21], confirming that genetic structure is modulated by farmer-directed selection pressure. Moreover, SSR markers divided Azizi accessions into two different genetic groups, which were initially collectively considered as a single landrace. However, 10 SSR markers were not sufficient to differentiate accessions of Souri and Roussia, as these were clustered together in a single genetic group. Therefore, we speculate that increasing the number of SSR markers would most likely improve the genetic differentiation between these landraces. Identification of 11 landraces, based on the spike and grain characteristics and in accordance with SSR fingerprinting data, suggests that morphological characterization is an efficient tool for varietal discrimination. Based on a landrace collection named by farmers, Mahmoudi and Biskri landraces were grouped into a single genetic group, as described by Robbana et al. [21]. This result highlights the importance of landrace identification based on a precise characterization of spike related traits and not solely on farmer-determined nomenclature. In fact, mixtures of different spike morphologies are often observed in a single field [19]. In addition, the mixture of landraces favors hybridization between different genetic groups, which explains the observed level of admixture (14%) herein. Most of the genetic groups displayed high level of genetic diversity, which is related to the high frequency of different MLGs. A predominance of a single MLG was found in Badri (G3) and Mahmoudi (G11), thereby reducing their genetic diversity. A predominant MLG of Mahmoudi landrace can be explained by the selection and multiplication, since 1908-1909, of a high-yielding accession aimed at increasing farmer income. In fact, the Mahmoudi landrace is preferred for its straw and grain yield as well as its ability to produce a high yield under drought and heat stress conditions prevalent in southern Tunisia. By contrast, Karim is the most popular modern variety in northern and central Tunisia, where heat and drought are not major problems. Badri is an old cultivar (released in 1969) obtained from a cross between two old landraces, Zenati Bouteille and Mahmoudi [11,13]. Therefore, the predominant MLG of Badri would correspond to previously released lines.
A high gene flow was observed between regions (Nm = 1.037) and climate stages (Nm = 3.813). In fact, Mahmoudi (G11), Beskri (G9) and Azizi (G1 and G8), together accounted for 47% of the entire durum wheat collection, were distributed across different geographical locations. The widespread distribution of these landraces is explained by their earlier introduction (since 1893 or 1894) from Algeria and southern Europe, followed by their spread through local seed commercial trade and seed exchange between farmers. Thus, the exchange of seeds of different landraces between farmers from distant regions, followed by the introgression of these landraces into the pre-existing germplasm, could explain the level of high genetic variation observed within regions (81%) and climate stages (90%). Gabes and Mahdia showed the highest diversity indexes (including the number of MLGs), where 80 and 77% of the accessions, respectively had a unique MLG; and the highest number of private alleles (at a frequency < 0.4). In Gabes, 72% of the accessions belonged to cluster C1, whereas in Mahdia, 22 and 60% of the accessions belonged to clusters C1 and C2, respectively. This result suggests that a more diverse germplasm resource was available to farmers in Mahdia than to farmers in Gabes for breeding purposes. In fact, in Mahdia, 27 accessions were identified as belonging to five different landraces with a frequency less than 41%, whereas, 38 accessions were described in Gabes, with a predominance of the Azizi landrace (47%). Moreover, a high morphological diversity was observed among regions (H′ = 0.55) due to the presence of different phenotypic classes for all the studied phenotypic traits within regions. According to Chentoufi et al. [30], the presence of wide morphological variability in wheat in different traditional agroecosystems could be explained by different seed exchange practices between farmers from different regions. Indeed, traditional management systems contributed effectively to the conservation of diversity of local durum wheat populations [36,37] and the maintenance of landrace varietal characteristics in Tunisia. Nevertheless, gene flow could be counter-balanced by farmer selection for preferred landraces, which could result in locally adapted accessions. In fact, a moderate genetic variability (19%) was observed among regions. This ascertainment highlights the effect of selection pressure directed by farmers, based on their preferences for specific wheat types in the preparation of local traditional dishes, which may have led to the adaptation and predominance of landraces in certain eco-geographical zones [7,22,28,30,38,39]. Notably, Kairouan, Sousse and Mahdia were characterized by specific landraces. For example, landraces Badri and Jneh Khotifa were only found in Kairouan, whereas landraces Taganrog and Richi were only found in Mahdia and Sousse, respectively. In addition, distinct phenotypic classes were detected within regions and climate stages (frequency > 70%), indicating that these classes might be characteristic of certain geographical areas, and environmental conditions may play a role in shaping the phenotypic diversity of durum wheat landraces. A local genetic adaptation pattern was also revealed in this study. For example, Gabes and Sousse were characterized by the presence of diagnostic alleles. Farmers in these regions have been selecting seeds and cultivating their old traditional landraces over many generations. This practice would result in the local adaptation of germplasm for a given eco-geographical environment [22]. Landraces under cultivation might undergo evolutionary changes if farmers keep using their own seed stock [38]. In fact, Fayaz et al. [40] defined landraces as locally adapted genotypes resulting from different environmental conditions and agricultural practices among ancient communities.
In addition to farmer's selection pressure, natural selection was found to affect morphological characteristics within a single landrace, Mahmoudi. Mahmoudi accessions collected from southern Tunisia showed significantly looser spikes than those collected from central Tunisia (compact spikes). We speculate that the loose spike, characterized by an open glume in southern Tunisian Mahmoudi accessions, could provide tolerance to high temperature by maintaining fertility, as shown in rice germplasm [41]. The loose spike trait of southern Tunisian Mahmoudi accessions could be used in breeding programs for developing heat stress tolerance, and for the identification of genes and mechanisms involved in flower development useful, thus improving wheat adaptation to arid and marginal environments.
MSN analysis grouped the accessions into two major clusters, C1 and C2. However, neither one of these clusters correlated with the geographical origin of landraces. Notably, both Mahmoudi and Biskri were introduced in Tunisia from Algeria, while Jneh Khotifa, Azizi, Mekki, Biada and Roussia were considered local landraces cultivated mainly in northern and central Tunisia. Although various origins have been reported for landraces Azizi and Mekki, no origin has been reported for Richi and Taganrog, which are very old, but non-local, landraces [12,13]. According to Deghais et al. [13], the landrace Jneh Khotifa was also known as Jneh Zarzoura and/or Kahla; the denomination of landrace Souri was extended in 1915 to Sarebouza obtained from Armenia. Soriano et al. [27], used 44 SSRs to show that Tunisian durum wheat landraces have four geographical origins, including East Mediterranean, East Balkan and Turkey, West Balkan and Egypt, and West Mediterranean, with dominance (> 50%) of the West Mediterranean genetic group. In addition, Sorriano et al. [27] demonstrated that western Mediterranean landraces are characterized by the heaviest grain weight compared with the other three genetic groups. In the current study, grain size did not significantly differ between C1 and C2 accessions, suggesting that both clusters contain accessions of the western Mediterranean origin. Moreover, Robbana et al. [21] reported that most of Tunisian landraces were introduced from the early Carthage trade maritime activity in the Mediterranean Sea through Lebanon, Greece and Italy.

Conclusions
Tunisian old durum wheat, characterized here by both high genetic and morphological diversity, represents an important and valuable genetic resource that should be included in breeding and well-established conservation programs. In this study, we showed that Tunisian old durum wheat is structured into landraces, revealing the effect of selection pressure directed by farmers for specific wheat types and agro-morphologies. Nevertheless, the morphogeographical spike density trait, apparent specifically in Mahmoudi accessions, suggests that environmental selection also acted on Tunisian durum wheat. Thus, our results provide important data for improving the adaptation of wheat to extreme or fluctuating Mediterranean conditions. Further physiological and agronomic analyses are needed to ascertain whether the spike density trait could be exploited in durum wheat breeding programs for tolerance to heat and drought.

Collection and multiplication of local durum wheat accessions
A collection of 304 old durum wheat accessions provided by the National Gene Bank of Tunisia (BNG) were used in this study. Accessions were collected from five regions in three distinct climate stages: Sousse and Mahdia (low semi-arid climate) and Kairouan (higharid climate) located in central Tunisia, and Gabes and Medenine (mid-arid climate) located in southern Tunisia. The Global Positioning System (GPS) coordinates of 163 out of 304 accessions were registered (Table S10). Seeds of each accession were sown and increased from a single spike-derived lineage by the BNG team, and a BNG code was assigned to each accession. All accessions were further multiplied for spike characterization. All accessions used in this study have been preserved at the BNG of Tunisia and are available upon request.

DNA extraction and SSR marker-based genotyping
Five seeds collected from one spike of each accession were germinated and grown under controlled conditions (20°C day/16°C night temperature, 16-h light/8-h dark photoperiod and 70% relative humidity) at Bioger research unit, INRAE, France. At the seedling stage (Zadok scale: [13][14], one leaf of each accession was sampled and placed in an extraction plate. The plates were placed at − 80°C for 12 h before DNA extraction. DNA was extracted from the leaf samples of each of the 304 accessions using the DNeasy PowerPlant Pro HTP 96 Kit (Qiagen, France). DNA concentrations were quantified using a Nanodrop spectrophotometer (ND-1000) and stored at − 20°C until needed for subsequent processing. The DNA of each accession was adjusted to a concentration of 15 ng·μl − 1 and genotyped using 10 SSR markers (Table S4), which were selected from a collection of 15 SSR markers used previously [29]. Forward primers were labeled with fluorescent dyes, and SSR markers were multiplexed, as described by Gautier et al. [42]. Each multiplex PCR, followed by gel electrophoresis, was performed according to the protocol established by Eurofin (https://www.eurofins.fr). Briefly, PCR was performed by preheating the DNA at 95°C for 5 min, followed by 35 cycles of 95°C for 30 s, 60°C for 90 s and 72°C for 30 s, with a final extension step of 60°C for 30 min. PCR products were analyzed by electrophoresis on a 2% agarose gel, and fragments were separated according to their size on an ABI Prism Genetic Analyzer (Applied Biosystems). Data was checked again using the Peak scanner software (version 1.0) (https://www.thermofisher.com). Two accessions with missing data for all SSRs were excluded from the study.

Morphological characterization of durum wheat accessions
Morphological characterization was carried out using five spikes per accession (total 1520 spikes). Accessions were evaluated using 12 quantitative and qualitative spike-and grain-related morphological traits. Spike evaluation was based on durum wheat descriptor standards of the IPGRI [23] and UPOV [24] (Table S11). Spike and grain morphological traits, defined by distinct phenotypic classes, were visually and numerically estimated. Traits including SC, GlC, AC, GC, SD, SS, GSp, GSz, and AL relative to SL were visually assessed, whereas other traits such as GSz, SL, NS and GN were quantified and then converted into codes. Subsequently, all accessions were named based on the catalog of cereal varieties cultivated in Tunisia [13]. This catalog serves as a reference for reporting and describing typical varietal characteristics of more than 40 old local durum wheat landraces cultivated in Tunisia.

Data analysis Polymorphism of SSR markers
Based on the SSR data generated on 302 accessions, the number of MLGs were identified with using the GIMLET software (version 1.3.2) [43]. To check the resolution of the 10 SSR markers used in this study, a genotype accumulation curve, was generated under R software [44] using the package 'pegas' package and 'genotype_curve' function in the R 3.3.2 [44]. This analysis determines the minimum number of loci necessary to discriminate between the Tunisian durum wheat genotypes by randomly sampling up to n − 1 loci (without replacement) and counting the number of MLGs observed.
To assess the informativeness of SSR markers, the average PIC value of each marker was calculated by determining the frequency of alleles per locus using the following equation [45]: where fi is the frequency of the i th allele in the set of 302 genotypes. SSR markers with PIC ≥0.5 were considered informative.

Polymorphism of morphological traits
Frequencies of different phenotypic classes were calculated for each of the 12 spike-and grain-related traits in the entire collection, by landraces (Table S1), regions and climate stages (Table S9). Based on these frequencies, the Shannon-Weaver index (H′) was calculated for each trait using the Past software [46]. H was estimated for the entire durum wheat collection, for accessions in different regions and climate stages and for each landrace. Each value of H was standardized by conversion to a relative phenotypic diversity index (H′) to express the values of H′ within the range of 0-1. The Shannon-Weaver index (H′) was calculated as follows: where H max Ln (S), S = the number of phenotypic classes.

Morphological and genetic structures
To investigate the morphological structure of the 304 accessions, PCA was performed based on 12 spike and grain related traits using R 3.3.2 [44]. The population genetic structure of these accessions was analyzed based on MLGs using STRUCTURE software (version 2.3.4) [26]. The STRUCTURE program was run with K values ranging from 1 to 20 in an admixture ancestry model by applying 10 independent runs for each K value, 100,000 burn-in phase iterations and 100,000 Markov Chain Monte Carlo (MCMC) iterations. The run with maximum likelihood was used to assign individual genotypes to different genetic groups. Genotypes with affiliation probabilities (inferred ancestry) > 75% were assigned to a distinct genetic group, and those with inferred ancestry < 75% were treated as admixed. The optimal number of genetic groups was determined using the mean posterior probability (ln P(D)) value per cluster (K) and the delta-K method of ln P(D) STRUCTURE harvester (version 0.6.94) [47]. To classify the 302 accessions according to their genetic relationship, MSN analysis was conducted based on Bruvo's distance [48] using 'poppr' and 'adegenet' packages in R 3.3.2 [44]. Furthermore, the average value of each of the 12 traits was calculated for accessions belonging to the different clusters, as defined by the MSN analysis, using the following equation: where N is the number of genotypes per genetic cluster, as defined by the MSN analysis; n is the number of individuals per phenotypic class; and C i is the i th phenotypic class per morphological trait.
To determine significant differences between genetic clusters for each morphological trait, ANOVA of calculated means was carried out using R 3.3.2 [44].

Analysis of population structure based on different regions and climate stages
Values of Na, Ne, PA (alleles specific to a single population), I, He, Ho, Fis, P, and DA (rare alleles with frequency > 70% for a genetic group or region and < 30% for others) were calculated within each genetic group, region and climate stage using GenAlEx (version 6.501) [25]. In addition, the correlation between genetic distance and log-transformed geographic distance (1 + geographic distance) of accessions was analyzed using a Mantel test [49] for the entire collection with GenAlEx (version 6.501) [25]. Correlations between the genetic distance matrix and morphological distance matrix were also assessed using a Mantel test. Furthermore, AMOVA was performed using GenAlEx (version 6.501) [25] to investigate the significance of genetic group differentiation (as defined by STRUCTURE) and genetic variability explained by regions and climate stages.
Moreover, mean values of all 12 spike and grain related traits were estimated for Azizi and Mahmoudi accessions located in different climate stages of central and southern Tunisia. To test potential regional effects on morphological traits, ANOVA of mean values was conducted in R 3.3.2 [44].