Use of ISSR markers to assess the genetic diversity of an endemic plant of Morocco (Euphorbia resinifera O. Berg)

Background Euphorbia resinifera is a melliferous, medicinal, and endemic plant to Morocco. Nevertheless, its ecological and genetic diversity still unknown. The objective of this study is to analyze the diversity and genetic structure of Moroccan wild populations of E. resinifera using ISSR markers. Twelve natural populations collected from its geographical range in Morocco were analyzed using 14 ISSR primers. Results A total of 125 bands were obtained, with polymorphism of 74.81%. The polymorphic information content (PIC), resolving power (Rp), Shannon’s information index (I), and total genetic diversity (Ht) were 0.33, 2.8, 0.35, and 0.21, respectively. The analysis of molecular variance showed that 75.56% of the total variability is present within populations and that 24.44% exists among populations. Also, the analysis showed a very low genetic differentiation between groups of mountain range type (FCT = 0.066), mountain versant type groups (FCT =  −0.024), and altitude groups (FCT =  −0.022). Moreover, the geographical distances between populations are correlated with their corresponding genetic distances according to the Mantel test (r = 0.507; P < 0.0001). Conclusion These results suggest that the population structuring follows a model of isolation by geographical distance. Indeed, the genetic structuring of populations into two groups obtained from PCoA and structure analysis revealed a dependence on the geographical origin of the populations. By contrast, the genetic distances are not correlated with the altitude.


Background
With nearly 2000 species, the genus Euphorbia is the largest of the Euphorbiaceae family and the second largest genus of flowering plants in the world [1,2].Despite its great vegetative diversity [3], the genus is united morphologically by the possession of a cyathium, a highly reduced inflorescence that resembles a single flower [4].In Morocco, few species exist, with a notable proportion of them being endemic [5].Euphorbia resinifera O. berg is one among them, growing naturally in the Atlas Mountains, particularly in Beni Mellal and Azilal provinces [5].The plant covers the mountain in a very discontinuous way from the High Atlas Mountain (Demnat) to the Middle Atlas Mountain (El Ksiba) [6].It presents a cushioned physiognomy constituted by a bush of juxtaposed stems.On average, the growth in height does not exceed one meter; however, their lateral growth is more important ranging from 0.5 to 2 m, which contributes to the decrease of the soil erosion.The mating system of species is primarily allogamous.Indeed, the Cyathes are arranged by 3 in axillary cyme stalks, towards the end of stems and branches; Cyathes are lateral hermaphroditic, stalked with thick peduncle; the median subsessile and male usually develop first and fall before the maturity of the capsules of the lateral cyathiums.E. resinifera employs both sexual and vegetative means of reproduction to ensure the survival and expansion of its population [7].The yellow flowers of this wild species attract and feed the bees.The beekeepers came from all regions of Morocco for take place near of the E. resinifera populations during the flowering season.Thus, with an annual honey production of about 300 tons, it is an important plant of the solidary agriculture in the region.It is well known as a melliferous plant for its high therapeutic and nutritional quality honey [8,9], which has been recognized as a local product of the region and has been labeled Protected Geographical Indication (PGI) [10].
Additionally, E. resinifera was recognized as having a vast therapeutic properties and benefits.In folk medicine, it is used to treat some types of complicated dermatoses, cancer treatment and control glycaemia in type II diabetics [11][12][13].In 1975, the resiniferatoxin major element in the latex of E. resinifera was identified [14].This molecule has many potential medical applications, analog to capsaicin but 1000 more potent than it [15].Moreover, many compounds of the latex of this species are endowed with biological activities as an anticancer [16,17], antioxidant, antibacterial, antifungal [18,19], an anti-pain and antituberculosis [20,21].
Nevertheless, the genetic resources of this species underwent important genetic erosion caused by various factors including the overgrazing, the deforestation provoked by a demographic pressure, voluntary or natural fires, and some fungal diseases attacking the plant.Consequently, emergency measures should be taken to safeguard this wild species.Thus, the safeguard and valorization of these local plant genetic resources constitute not only an imperative but also a major component of the stability of natural ecosystems and their valorization in agro-economic, therapeutic, and environmental perspectives.For that, the description of Moroccan E. resinifera populations may help to identify different genotypes and rationalize conservative treatments.Therefore, it has become imperative to establish a research program aiming at the evaluation of the genetic diversity of this species.Our first study on morphological characters showed high phenotypic diversity in twelve Moroccan populations [3].However, no information is available on its genetic diversity.For this, it becomes necessary to find more discriminating markers, which could provide information about the variability within and among populations of this species and investigate new resources of variation which might be used for specific conservation programs.
Several types of molecular markers are employed for assessment of genetic diversity and relationships in plant species, including RFLPs (restriction fragment length polymorphisms) [22], RAPDs (randomly amplified polymorphic DNAs) [23][24][25], AFLP (amplified fragment length polymorphisms) [26,27], SNPs (single-nucleotide polymorphic) [28], SSRs (simple sequence repeats) [29][30][31], and ISSRs (inter simple sequence repeats) [32][33][34][35].However, RAPDs have low reproducibility, RFLPs are time-consuming and labor-intensive, SSRs require the knowledge of the flanking regions for the development of species-specific primers, and AFLPs and SNPs have high cost [36], while the ISSR markers are hypervariable, highly reproducible, fast, inexpensive, and do not require any prior sequence information of amplified locus [36,37].ISSR is a kind of DNA sequences confined by two inverted SSR composed of the same motives, which are amplified by a unique PCR primer.ISSR-PCR detects the levels of variation in microsatellite regions and gives multi-locus schemes, which are very iterative, plentiful, and polymorphic in plant genomes [38,39].
The present work provides the first data on the genetic diversity of Moroccan E. resinifera based on molecular markers.Our objectives were to give a preliminary estimation of the genetic diversity and structure of this species.Then, twelve natural populations of E. resinifera originating from diverse altitudes and geographical area were analyzed using fourteen ISSR primers.

Plant material and DNA extraction
The plant material used in this study included twelve natural populations representing the distribution area of E. resinifera in Morocco (Fig. 1).Six populations are originating from High Atlas Mountain and six belonging to the Middle Atlas Mountain.The geographic characteristics and meteorological conditions such as altitude slice, central latitude and longitude as well as the mean precipitation and temperature average of these populations are provided in Table 1.From each population, five bush of similar age were randomly sampled and the collected young stems were stored at −20 °C until DNA extractions.
The DNA was extracted following the method described by Doyle and Doyle [40] slightly modified [41].The quality of extracted DNA was examined by electrophoresis on 1% agarose gel, and DNA quantity was determined spectrophotometrically.Later, the samples were brought to a working concentration of 10 ng/μl.

ISSR analyses
A total of 14 ISSR primers previously displayed reliable and polymorphic band profiles [42,43] were used in this work (Table 2).The amplification reactions were performed in a volume of 12.5 μl, which contained 6.25 μL 2 × Green Taq Mix (Vazyme, Nanjing, China), 0.5 μL of each primer (0.4 μM), 5 μL of template genomic DNA (4 ng/μL), and 0.75 μL of distilled deionized water.PCRs were conducted in a DNA thermocycler (Multigene gradient, Labnet, NJ.USA).The PCR program consisted of an initial denaturation at 94 °C for 5 min, followed by 35 cycles according to the following procedure: denaturation at 94 °C for 30 s, annealing at determined temperature for 30 s, extension at 72 °C for 1 min, the last cycle was followed by a final extension for 7 min at 72 °C.Gradient PCR was used to determine the annealing temperature of each primer (Table 2).PCR products were separated by electrophoresis on 1.7% agarose gel submerged in 0.5 × TBE buffer and stained with 1 μg/μl of ethidium bromide.The DNAs were visualized under UV light using the Gel Doc system (EnduroTM GDS, Labnet).The fragment size was estimated by using a DNA marker (100 bp plus II DNA ladder, TransGen Biotech Co., Ltd).

Data analyses
The band profiles of each gel were scored visually and recorded as presence (1) or absence (0) of bands leading to the construction of data binary matrix (1,0).For each primer, the percentage of polymorphic band (PPB), the polymorphic information content (PIC), and resolving power (Rp) were determined.Also, the POPGENE software was used to measure the following parameters: numbers of alleles (Na), effective number of alleles (Ne), genetic diversity within populations (Hs), total gene diversity (Ht), coefficient of gene differentiation (Gst), and Shannon's information index (I).Partition of the observed genetic variation and calculation of the corresponding F values were carried out using different hierarchical analysis of molecular variance (AMOVA).Firstly, global analysis of AMOVA was done to apportion the total genetic variation into two hierarchical levels: among populations (FST) and within populations.Secondly, hierarchical AMOVA analysis was used to partition the variation among: -Two mountain range groups of populations: Middle Atlas (populations: KSB, TAG, AAS, OAY, AFO, and MOD) and High Atlas (populations: BIN, BZO, FMJ, IMI, WAW, and OUZ) -Three mountain versant groups of populations: North versant (populations: BIN, AAS, OAY, and AFO), Southwest versant (populations: OUZ, TAG, FMJ, and IMI), and South versant (populations: BZO, MOD, WAW, and KSB) -Two altitude groups of populations: very low (populations: BZO, OAY, TAG, AAS, FMJ, BIN, OUZ, and AFO) and low (populations: KSB, MOD, IMI, and WAW) (Table 1)   Population's specific FST indices were also calculated to investigate which population is more divergent from the remaining.The average gene diversity over loci was calculated to estimate intra-population variability.These evaluations were performed using the package ARLE-QUIN version 3.01 [44].The pairwise genetic differentiations (FST) between the twelve populations were also generated by AMOVA and the gene flow (N e m) was approximated estimated through Wright's island model: N e m = 0.25 (1/FST -1) [45].A Mantel test was used to test whether matrix of genetic distances (FST) between populations was significantly correlated with their corresponding matrices of geographic distances and difference of altitude (1000 permutations; routine MXCOMP of the NTSYS-pc; package).A Bayesian structure analysis was executed using the STRU CTU RE v.2.3.4 software to deduce the population genetic structure and define the number of groups within the studied populations [46].To identify the number of K clusters explaining the observed genetic structure, we used the STRU CTU RE Harvester website [47], which implements the Evanno method [48].Furthermore, to elucidate the populations' relationships, principal coordinate analysis (PCoA) was carried out using the DARwin software version 6.0.02 [49].

Genetic diversity
Estimates of genetic diversity of studied populations are summarized in Table 3.The results showed that the number of observed alleles (Na) was fixed in the value of 2 for all primers.The highest effective number of alleles (Ne) varied from 1.09 (UBC853) to 1.56 (UC834) with a general mean of 1.34 allele per primer, while Shannon's Information index (I) showed a lowest value (0.17) for UBC853 and the highest value (0.54) for UBC834, with an average of 0.35.In addition, the total genetic diversity (Ht) oscillated from 0.08 for UBC853 to 0.33 for UBC834 with an average of 0.21.The genetic diversity within species (Hs) ranged from 0.06 for UBC853 to 0.31 for UBC834 (mean 0.15).Moreover, Nei's coefficient of genetic differentiation between populations (G ST ) varied from 0.12 (UBC834) to 0.44 (UBC836) with an average of 0.32, indicating that 32% of total genetic variability was distributed among populations and the remaining (68%) accounted for within populations.These results were congruent with that revealed by FST value, accounting for 0. 244 (Table 4).The great level of genetic differentiation among population of this wild plant is in accordance with the low value of gene flow estimated (N e m = 0.77) which provides information on amount of number of exchanged individuals between populations studied.When AMOVA was performed, at three hierarchical levels, with two mountain range type, very low genetic differentiation was observed between groups (FCT = 0.066) even though that it is slightly significant (Table 4).This indicates that mountain range type has had little effect on populations' structuration which implies that there is no local adaptation of studied populations.Also, a very low genetic differentiation was obtained between mountain versant type groups of populations (FCT = −0.244)and among altitude groups (FCT = −0.222).

Genetic relationship
The pairwise F ST values and geographic distances between the 12 populations are presented in Table 5.According to the Mantel test (r = 0.507; P < 0.0001), the geographical distances between populations are correlated with their corresponding genetic distances.These results suggest that the populations' structure follows a model of isolation by geographical distance.By contrast, no significant association was obtained between the genetic distance and difference of altitude between populations (r = −0.16,P = 0.19).Thus, the populations "AAS" and "FMJ" with a low difference of altitude (2 m) have a high value of genetic distance (0.376), and the populations "BZO" and "WAW" with a high difference of altitude (743 m) have obtained the low value of genetic distance (0.195).Between 66 pairwise F ST values, 54 values are significant, meaning that the populations are widely different from each other.The significant values varied from 0.084 (FMJ/BZO; 12) to 0.466 (OUZ/ KSB; 87), which mean that populations FMJ/BZO are the most genetically similar and OUZ/KSB are the most divergent.The result indicated that the Bin El Ouidane (BIN) and Bzou (BZO) populations have the lowest value of specific FST index (respectively 0.209 and 0.212), and Elksiba (KSB) population has the highest value of this index (FST = 0.289), implying that this latter population is the most divergent from the others studied.Regarding the intrapopulation variability assessed by gene diversity over loci (data not shown) which oscillated from 0.148 for Elksiba (KSB) to 0.261 for Bin El Ouidane (BIN) and 0.257 for Bzou (BZO), reflecting consequently for any in situ and/or ex situ conservation strategy should aim to include those later populations.
The structure analysis based on the ΔK method showed that the best number of genetic clusters (K) was 2, suggesting that all individuals fell into two clusters (Fig. 2).Moreover, based on the permuted average Q-matrix generated by Clumpak, the highest H′ was observed for K = 2 (H′ = 0.947), indicating the stability of the result for this model.Considering the genotypes as pure when the membership coefficient was greater than 0.80 and as a hybrid or admixture when the membership coefficient was lower than 0.80, 41 individuals among the 60 analyzed (68.33%) were assigned to one of the model's defined groups.The first group (red) was formed by the individuals from High Atlas Mountain populations, namely Bzou (BZO4), Ouled Ayyad (OAY2, OAY3, and OAY5), Bin El Ouidane (BO3 and BO4), all individuals   from populations of Ouzoud (OUZ1, 2, 3, 4, and 5), Imi n'Ifri (IMI5), and four from Wawla (WAW1, WAW2, WAW4, and WAW5), with a membership coefficient oscillated between 0.815 to 0.990, while the rest of this group were originating from Middle Atlas Mountain populations: Afourer (AFO2 and AFO4), Modj (MOD1), with a membership coefficient ranging from 0.953 to 0.980.Nevertheless, the other seven individuals, namely (TAG3, TAG5) from Tagzirt population, (AFO3, AFO5) coming from Afourer population, MOD3 from Modj population, IMI1 belonging to Imi n'Ifri population, and WAW3 from Wawla population could be considered as admixed (coefficients ranged from 0.507 to 0.644).The second group (green) contained individuals coming from High Atlas Mountain populations: (BO1, BO2, and BO5) of North versant, with a membership coefficient ranging from 0.954 to 0.975, all individuals of Foum Jemaa population (FMJ1, 2, 3, 4, and 5), and one individual (IMI2) from Imi n'Ifri population belonging to Southwest versant (coefficients between 0.840 and 0.959) and (BZO1, BZO2, and BZO3) collected from South versant having a membership coefficient between 0.922 and 0.986, while the rest of group include bushes coming from Middle Atlas Mountain populations: (AAS3 and AAS4) of North versant with an assignment coefficient of 0.888 and 0.86 respectively, (TAG1, TAG2 and TAG4) of Southwest versant with a membership coefficient from 0.932 to 0.969, and all bushes of Elksiba population (KSB1, 2, 3, 4, and 5) from South versant mountain, having a coefficient of assignment oscillating between 0.929 and 0.976.Finally, 12 remaining bushes which could be considered as admixture belonged to the OAY, AAS, AFO, MOD, IMI, and BZO populations (coefficients between 0.508 and 0.775).
The genetic structure of Moroccan E. resinifera populations was further reconstructed by using the PCoA.Indeed, about 24% of total variance was explained by the first two components and the plot of PCoA divided studied populations in two groups (I and II, Fig. 3), which corroborates the populations' structure obtained by Bayesian analysis (Fig. 2).Consequently, the genetic structure of all studied populations bushes in two main groups was operated independently of altitude and mountain versant type.

Discussion
DNA markers have become a useful tool for evaluating the genetic diversity of many different plant species.In this work, ISSR markers were used to assess the genetic diversity of E. resinifera populations in Morocco.The 14 tested ISSRs primers revealed a high percentage of polymorphism with an average of 74.81%.This high percentage implies that there is an important genetic diversity in this endemic species.This result is higher than that reported with ISSR markers for other Euphorbia species (E.khabrica, E. buhsei, E. osyridea, and E. austro-iranica) in Iran by Pahlevani et al. (13.02%, 20.71%, 24.85%, and 14.20%, respectively) [50], in Saudi Arabia by Moustafa et al. (Euphorbia prostrata Aiton, Euphorbia peplus L., and Euphorbia terracina L.) (20.28%, 14.08%, and 11.44%, respectively) [51], and for 15 Euphorbia species by El-Hawary et al. (total polymorphism = 57.7%)[52].However, our finding was lower than that observed by Reginaldo et al. for Brazilian Euphorbiaceae (Croton urucurana Baill) (89%) based on ISSR markers [53] and by Dorset et al. for American E. telephioides (80.7%) revealed by allozyme markers [54].In addition, the high values of PIC (0.33) and Rp (2.8) parameters show that the ISSR markers are very informative and efficient for analyzing the diversity and genetic structure of E. resinifera.These values are comparable with these obtained by Reginaldo et al. for Brazilian Euphorbiaceae (Croton urucurana Baill) using ISSRs markers (PIC = 0.29 and Rp = 3.4) [53].Moreover, the high multi-locus value of Ht (= 0.21) suggests the presence of a high level of polymorphism of this endemic species.Indeed, this high polymorphism is confirmed by the Shannon index value (I = 0.35).This value is higher than that obtained by Reginaldo et al. (I = 0.26) for Brazilian Euphorbiaceae using ISSRs markers [53].Besides, the gene diversity within the species (Hs) was  E. resinifera populations is in agreement with the general trend for allogamous and long-lived woody perennial species (Ht = 0.28) and for angiosperm species (Ht = 0.28) [55].The wild populations of E. resinifera were largely differentiated (Gst = 0.32, FST = 0.244) which could be due to restricted gene flow (N e m = 0.77) between populations.This finding is in concordance with that observed for other herbaceous outcrossing perennial plant species (Gst > 0.20) [56].The Gst value (0.32) detected in this investigation is higher than that recorded by Dorset et al. for the endemic North American species Euphorbia telephioides (0.10) using allozyme markers [54].Similar results were revealed by Ki-Ryong for Euphorbia fauriei in Korea (Fst = 0.237; N e m = 0.80) and Euphorbia jolkinii in Taiwan (Fst = 0.245; N e m = 0.77) [41].The indirect estimate of gene flow on the basis of FST was low (Nem = 0.77) and might be due to the presence of geographical barriers between populations such mountains.Indeed, many factors such as discontinuous distribution of populations, limited pollinator movements, and low rate of seed migration could be an efficient obstacle to the gene flow and the origin of the high value of differentiation among E. resinifera populations.
The analysis of molecular variance (AMOVA) showed that 75.56% of the diversity accounted for within population leaving 24.44% for among populations.The existence of high genetic variability within population should help the population to cope with local environmental changes.This result implies that sampling from a small number of populations, particularly those with high intrapopulation variability, is sufficient for species in situ and/or ex situ conservation.Bin El Ouidane (BIN) and Bzou (BZO) populations with high intrapopulation variability are more convenient for this purpose.
Moreover, hierarchical AMOVA revealed a very low genetic differentiation between the two groups of mountain range (FCT = 0.066), even though that it is slightly significant, which suggests that mountain range type have had little influence on structuration's populations.These findings give choice to sampling from populations of Middle or High Atlas Mountains for conservation and breeding species.Likewise, when assembling the populations according to their mountain versant type and altitude, a very low genetic differentiation was obtained between respective groups (respectively, FCT = −0.024and FCT = −0.022),indicating that mountain versant type and altitude did not have an effect on population Strangely enough, genetic distances were not correlated to the difference of altitude between the populations (r = −0.16,P = 0.19).A similar result was reported for Moroccan walnut using SSR marker [57].Also, other studies described no genetic differentiation among populations at low and high altitudes [58,59], probably due to the overlap of flowering phenology in populations at different altitudes, species' extensive pollen flow, and long distance seed dispersal among different altitudes by animals especially birds.This result is confirmed by the biased model and PCoA which showed that the twelve populations are gathered in two groups undependably to mountain range, mountain versant type, and altitude.In contrast, geographic distances have explained the genetic differentiation between populations according to the Mantel test (r = 0.507; P < 0.0001).This result suggests that the population structure follows a model of isolation by geographic distance.Consequently, our finding suggests that more closely situated populations tend to be more genetically similar to one another.

Conclusion
The present study is the first work aiming to evaluate the genetic diversity and structure of E. resinifera populations in Morocco using ISSR markers.The results of this study confirmed that ISSR markers could be powerful tools for detecting genetic diversity among and within E. resinifera populations.The level of genetic diversity was high, and the genetic variation mainly existed within populations.The results led to structure of populations in two gene pools independently of mountain range type, mountain versant type, and altitude.Based on these results, it could be sufficient to sample from a few populations, particularly those most genetically diversified, for any in situ and/or ex situ conservation strategy.

Fig. 2
Fig. 2 Bar plots represent STRU CTU RE inferences of individual assignments (N = 60, K = 2, ΔK = 13.71,H' = 0.947).Each individual is represented by a vertical bar showing membership coefficients to each genetic cluster

Fig. 3
Fig. 3 PCoA plot of the E. resinifera individuals showing relationships based on ISSR data

Table 1
Geographic and ecological characteristics of E. resinifera populations used in the study

Table 2
Properties of ISSR markers used in the study and statistical parameters: polymorphism information content and resolving powers AT°C Annealing temperature in C°, Y (C,T), R (A,G), PPB Percentage of polymorphic bands, PIC Polymorphic information content, Rp Resolving power ISSR

Table 3
Genetic diversity revelated by the 14 ISSR primers Na Number of observed alleles, Ne Effective number of alleles, I Shannon's Information index, Ht Total genetic diversity, Hs Genetic diversity within group, Gst Genetic differentiation among group

Table 5
Pairwise F ST values (below diagonal) and corresponding geographic distances in km (above diagonal) and corresponding difference of altitude (in m, above diagonal in bold) for twelve populations of E. resinifera analyzed by ISSRs.Abbreviations as in Table1