Introduction
The genus Dracocephalum L. (Lamiaceae) contains about 60-70 species. These species are mostly perennial herbs, and rarely annual taxa growing in alpine and semi-dry regions mainly in temperate Asia, with a few species occurring in Europe, and one species in North America (Brach and Song, 2006). Eight Dracocephalum species have been reported in Iran that grow in north and central parts of the country (Rechinger, 1982). These species have medicinal values including anticancer, antioxidant, anti-hypoxic and immunomodulatory activities (Zeng et al., 2010).

Population genetics study is an important step in the way of planning genetic and breeding programs for crop plants and medicinal plant species. It provides data on the genetic variability, gene flow versus population genetic isolation, population genetic fragmentation, and the role of genetic drift, the bottleneck and any other evolutionary forces acting on population divergence (Sheidai et al., 2012; 2013).

Dracocephalum thymiflorum L. grows in limited areas in Iran and forms a few local populations. We didn’t have any information available about genetic diversity and population genetic structure of these medicinal species in the country. Therefore, the present population genetic study was performed to provide data on population genetic structure, potential gene pools and gene flow in Dracocephalum thymiflorum in Iran. The information obtained can be used in conservation of this medicinal plant species. Different molecular markers have been used in population genetics studies. We used ISSR (Inter-simple sequence repeats) molecular markers, as they are reproducible, simplified the work and not expensive (Sheidai et al., 2014).

1 Results
1.1 Genetic diversity analysis
Genetic diversity parameters determined in 5 studied populations of Dracocephalum thymiflorumis are presented in Table 2. The highest value for He (0.243), polymorphism percentage (P% = 73.91) and genetic diversity due to population (Hs = 0.294) occurred in population 1 (Amol population). The lowest value for the same parameters (0.120, 47.83, and 0.195, respectively) occurred in the population 5 (Ramsar population).
 

AMOVA test revealed significant molecular difference (PhiPT = 0.27, P = 0.01) among the studied populations. It showed that 25% of total genetic variability occurred due to among populations genetic difference, while 75% due to within population genetic variability. Pairwise Fst values obtained for the studied populations were significant for most of the studied populations (P = 0.01). This indicated genetic divergence of all studied populations (Table 3). High Hickory theta B value (0.35) obtained supported AMOVA result.
 

The population genetic differentiation analysis revealed G'st_est = 0.283, P = 0.001, and D_est = 0.010, P = 0.001 and CVA plot separated each population based on its genetic variance (Figure1). These results indicate that the populations studied are genetically differentiated.
 

1.2 Populations’ genetic affinity
The grouping of the populations by NJ tree and Ward clustering produced similar results. Therefore, the Ward dendrogram is only presented here (Figure 2). Almost all plants of each population were grouped together in a distinct cluster. This result is in agreement with AMOVA result and revealed the populations, genetic divergence in Dracocephalum thymiflorum. Moreover, it showed the use of ISSR molecular markers in population genetic studies of this species.
 

Two major clusters were formed. The populations 1 and 2 showed higher degree of genetic affinity and were placed close to each other in the first major cluster. The populations 3-5 were placed in the second major cluster and showed some degrees of intermixture.

The PCoA plot is presented in Figure 3 The plants of population 2 have been distributed in the left corner of this plot, separated from the other studied populations. This showed genetic difference of this population from the other studied populations. Moreover, this population showed a high level of within-population genetic variability, and its plants are distributed from top to down of the PCoA plot.
 

The plants of population 1 were placed close to the population 2 as also revealed by Ward dendrogram. However, some of its plants were placed far due to genetic difference. The plants of population 1 also revealed a good level of within-population genetic variability.

Plants of population 3-5 were placed close to each other, with plants of population 3 intermixed with the plants of populations 4 and 5. This population also showed high level of genetic variability as they were scattered from left side to right side of the PCoA plot. In general, PCoA plot revealed high within and between population genetic variability in Dracocephalum thymiflorum population, that is good from conservation point of view.

1.3 Population genetic structure and gene flow
The Evanno test and K-Means clustering produced the best number of genetic groups as K = 2. Therefore, STRUCTURE analysis was performed for K = 2. The STRUCTURE plot obtained (Figure 4) revealed close genetic affinity between the studied populations due to ancestral shared alleles. However, it also revealed that population 2 is genetically differentiated in its genetic structure. The plot also showed a higher degree of genetic admixture among populations 3-5, followed by population 1. This Bayesian approach analysis is in agreement with PCoA plot result.
 

Population assignment test revealed more detailed information on genetic admixture and gene flow among these populations (Table 4). This test is based on maximum likelihood of plants membership to their own or other populations. It revealed that out of 55 studied plants, 16 plants were inferred to be from other populations. The population 2 was comparatively less involved in gene exchange with the other populations and the main stream of gene flow occurred between populations 1, 3, 4, and 5.
 

The Mantel test produced significant correlation (r = 0.41, P = 0.01) between genetic distance and geographical distance of the studied populations. Therefore, the populations that are in closer vicinity had the chance for gene flow between each other. LFMM analysis identified 12 out of 46 ISSR loci are adaptive.

2 Discussion
With increase in the size of human population, the crop plants and medicinally important plant taxa are consumed and destroyed faster that before. This treat is greater for those plant species that are rare and grow in limited number and therefore, their conservation become an important task (Sheidai et al., 2013, 2014). Medicinal plants such as Dracocephalum thymiflorum are extensively used by locals and therefore are subject to be reduced in number or elimination from the natural habitat. The disappearance and fragmentation of natural populations could lead to reductions in the rate of gene flow among populations. This in turn increases in genetic differentiation among populations and reductions in genetic variation within populations due to genetic drift (Setsuko et al., 2007; Hou and Lou, 2011).

To plan a properly oriented conservation plan, the knowledge of genetic diversity available in the target species becomes important. The genetic variability can help plant taxa to adapt to changing environments they are growing in (Freeland et al., 2011).

The AMOVA test indicated that out of total genetic variation, 75% was due to within population genetic variability in the studied Dracocephalum thymiflorum populations. This should be related to outcrossing nature of this plant species. The presence of high within population genetic variability helps the population to cope with local environmental changes.

STRUCTURE analysis and K-Means clustering revealed some degree of population genetic fragmentation in the studied populations. This was mainly due to the plants in population 2 that differed genetically from the others.

Moreover, AMOVA, Gst and differentiation parameters revealed significant genetic difference among the studied populations. Among population genetic differentiation is the product of absence and or limited between population gene flow, genetic drift, inbreeding, and local adaptation (Hou and Lou, 2011; Sheidai et al., 2014).

The assessments of the levels of within- and among-population genetic variation have been used to prioritize for conservation efforts (Petit et al.,1998) with, all else being equal, more weight given to those populations exhibiting higher levels of within-population variation, and to those that are more genetically divergent from others. In our study the population 2 was genetically differentiated from the other studied populations and also contained a high degree of within population genetic variability.

These populations may have increased likelihood of persistence over less variable population and hence the ability of a population to contribute demographically to the species through time, and have increased adaptability in the face of future environmental change.

The STRUCTURE plot and population assignment revealed some degree of genetic admixture among the studied Dracocephalum thymiflorum populations. Therefore, in spite of population differentiation, limited gene flow occurred among these populations.

Gene flow is also important in conservation contexts, particularly for the species with few local populations. In these species, the genetic characteristics are strongly influenced by genetic drift and inbreeding (Frankham et al., 2002). Fortunately, although Dracocephalum thymiflorum populations are few in number and are confined to some ecological places, they showed good within-population genetic variability and limited amount of among population gene flow. Gene flow among local populations could mitigate losses of genetic variation caused by genetic drift in local populations and thus save them from extinction (Richards, 2000).

Mantel test revealed isolation by distance in the studied Dracocephalum thymiflorum populations. The plant species that form geographical populations, as geographical isolation increases, a reduction in both seed dispersal and pollen flow will result in decreased gene flow between populations. The resulting genetic isolation may lead to pronounced geographical structuring in genetic variation within a species as population differentiation increases (Jump et al., 2003). High degree of genetic variability observed in Dracocephalum thymiflorum populations can be used in programming conservation plan of this medicinal plant species in Iran.

We conclude that, genetic differentiation, genetic drift, limited gene flow and local adaptation have played role in Dracocephalum thymiflorum population divergence.
 
3 Materials and methods
3.1 Plant materials
Fifty-five plant specimens were randomly collected from 5 geographical populations of Dracocephalum thymiflorum L. Details of localities are provided in Table 1. Voucher specimens are deposited in Herbarium of Shahid Beheshti University (HSBU). Fresh leaves were collected and used for DNA extraction and molecular study.
 
3.2 DNA extraction and ISSR assay
Fresh leaves were collected randomly in each of the studied populations and dried in silica gel powder. Genomic DNA was extracted using CTAB activated charcoal protocol (Sheidai et al., 2013). The quality of extracted DNA was examined by running on 0.8% agarose gel.

Ten ISSR primers; (AGC)5GT, (CA)7GT, (AGC)5GG, UBC810, (CA)7AT, (GA)9C, UBC807, UBC811, (GA)9T and (GT)7CA commercialized by UBC (the University of British Columbia) were used. PCR reactions were performed in a 25 μl volume containing 10 mM Tris-HCl buffer at pH 8; 50 mM KCl; 1.5 mM MgCl2; 0.2 mM of each dNTP (Bioron, Germany); 0.2 μM of a single primer; 20 ng genomic DNA and 3 U of Taq DNA polymerase (Bioron, Germany). The amplifications, reactions were performed in Techne thermocycler (Germany) with the following program: 10 Min initial denaturation step 94°C, 30 S at 94°C; 1 Min at 57°C and 1 Min at 72°C. The reaction was completed by final extension step of 7 Min at 72°C.

The amplification products were visualized by running on 2% agarose gel, followed by the ethidium bromide staining. The fragment size was estimated by using a 100 bp molecular size ladder (Fermentas, Germany).

The ISSR analysis performed was tested for marker reproducibility/genotyping errors, by repeating the experiment 3 times and scoring only the sharp and consistent bands obtained (Pompanon et al., 2005). Moreover, the initial DNA material was obtained from 10 randomly collected leaves in each replication. Although we did not calculate error rate for the studied samples, we performed Hickory test for ISSR data analysis that is Bayesian approach based method and run the program for 3 times for consistency. Moreover, we used no DNA sample in ISSR-PCR for each primer.

3.3Data analyses
3.3.1 Genetic diversity and population structure
ISSR bands obtained were coded as binary characters (presence = 1, absence = 0). The genetic diversity parameters like, Nei’s gene diversity (H), Shannon information index (I), number of effective alleles, and percentage of polymorphism (Freeland et al., 2011; Weising et al., 2005), were determined. Nei’s genetic distance was used for clustering (Weising et al., 2005; Freeland et al., 2011). Neighbor Joining (NJ) and Ward clustering as well as PCoA (Principal coordinate analysis) were used for population grouping, after 1000 times bootstrapping/ and or permutations (Podani, 2000; Freeland et al., 2011).The Mantel test was performed to check correlation between geographical and the genetic distances of the studied populations (Podani, 2000). PAST ver. 2.17 (Hamer et al., 2012) and, DARwin ver. 5 (2012) programs were used for these analyses.

AMOVA (Analysis of molecular variance) test (with 1000 permutations) as implemented in GenAlex 6.4 (Peakall and Smouse, 2006), and Nei, Gst analysis of GenoDive ver.2 (2013) (Meirmans and Van Tienderen 2004), were used to reveal significant genetic difference among the studied species (Sheidai et al., 2014).

The population genetic differentiation was studied by G'st_est = standardized measure of genetic differentiation (Hedrick, 2005), and D_est = Jost measure of differentiation (Jost, 2008) and CVA (Canonical variate analysis) method (Podani, 2000).

In order to overcome potential problems caused by the dominance of ISSR markers, a Bayesian program, Hickory (ver. 1.0) (Holsinger and Lewis, 2003), was used to estimate parameters related to genetic structure (theta B value) (Tero et al., 2003).

Bayesian based model STRUCTURE analysis (Pritchard et al., 2000), with 105 permutations and admixture method was used to study the genetic structure of populations (Sheidai et al., 2014). For STRUCTURE analysis, data were scored as dominant markers (Falush et al., 2007).

The optimum number of genetic groups (k) was determined by 1- Evanno test (Evanno et al., 2005) performed on STRUCTURE result and 2- K-Means clustering method (Sheidai et al., 2014).

3.3.2 Gene flow
Gene flow was determined by different approaches. 1- Calculating Nm an estimate of gene flow from Gst by PopGene version 1.32 (1997) as: Nm = 0.5(1 - Gst)/Gst. This approach considers equal amount of gene flow among all populations. 2- STRUCTURE analysis based on admixture model and Bayesian approach (Pritchard et al., 2000), and 3- population assignment test based on maximum likelihood as performed in Genodive ver. in GenoDive ver. 2. (2013).

Frichot et al. (2013) latent factor mixed models(LFMM) was used to check if ISSR markers show correlation with environmental features of the studied populations. The analysis was done by LFMM program Version: 1.2 (2013).

References
Brach A.R., and Song H., 2006, eFloras: New directions for online floras exemplified by the Flora of China Project, Taxon, 55(1): 188-192
http://dx.doi.org/10.2307/25065540

Bureš P., Wang Y., Horova L., and Suda J., 2004, Genome size variation in central European species of Cirsium (Compositae) and their natural hybrids, Annals of Botany, 94: 353-363
http://dx.doi.org/10.1093/aob/mch151

Evanno G., Regnaut S., Goudet J., 2005, Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study, Molecular Ecology, 14, 2611-2620

Falush D., Stephens M., and Pritchard J.K., 2007, Inference of population structure using multilocus genotype data: dominant markers and null alleles. Molecular Ecology Notes, 7: 574-578
http://dx.doi.org/10.1111/j.1471-8286.2007.01758.x

Frankham R., Ballou J.D., and DBriscoe D.A., 2002, Introduction to conservation genetics. Cambridge University Press, Cambridge, UK.

Freeland J.R., Kirk H., and Peterson S.D., 2011, Molecular Ecology (2^nded), UK: Wiley-Blackwell, pp. 449
http://dx.doi.org/10.1002/9780470979365

Frichot E., Schoville S.D., Bouchard G., and Francois O., 2013, Testingfor Associations between Loci and Environmental Gradients Using Latent Factor Mixed Models. Molecular Biology and Evolution, 30: 1687–1699
http://dx.doi.org/10.1093/molbev/mst063

Hamer Øm., Harper D.A.T.,and Ryan P.D., 2012, PAST: Paleontological Statistics software package for education and data analysis, Palaeontol Elect 4: 9

Hedrick P., 2005, Genetics of Populations, 3rd edn, Jones and Bartlett Publishers, Sudbury, MA.

Holsinger K.E., and Lewis P.O., 2003, Hickory: a package for analysis of population genetic data V1.0
http://www.eeb.uconn.edu

Hou Y., and Lou A., 2011, Population Genetic Diversity and Structure of a Naturally Isolated Plant Species, Rhodiola dumulosa (Crassulaceae), PLoS ONE, 6, e24497. doi:10.1371/journal, pone.0024497
http://dx.doi.org/10.1155/2014/476346

Jost L., 2008, GST and its relatives do not measure differentiation, Molecular Ecology, 17:4015–4026

Jump A.S., Woodward F.I., and Burke T., 2003, Cirsium species show disparity in patterns of genetic variation at their range edge, despite similar patterns of reproduction and isolation, New Phytologist, 160: 359-370
http://dx.doi.org/10.1046/j.1469-8137.2003.00874.x

Meirmans P.G., and Van Tienderen P.H., 2004, GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms,  Molecular Ecology Notes 4: 792- 794
http://dx.doi.org/10.1111/j.1471-8286.2004.00770.x

Peakall R., and Smouse P.E., 2006, GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research, Molecular Ecology Notes, 6: 288-295
http://dx.doi.org/10.1111/j.1471-8286.2005.01155.x

Petit R.J., Mousadik A., and Pons O., 1998, Identifying populations for conservation on the basis of genetic markers, Conservation Biology 12: 844–855
http://dx.doi.org/10.1046/j.1523-1739.1998.96489.x

Podani J., 2000, Introduction to the Exploration of Multivariate Data [English translation], Leide, Netherlands,Backhuyes

Pompanon F., Bonin A., Bellemain E., and Taberlet P., 2005, Genotyping errors: causes, consequences and solutions, Nature Reviews Genetics, 6: 847-846
http://dx.doi.org/10.1038/nrg1707

Pritchard J.K., Stephens M., and Donnelly P., 2000, Inference of population structure using multilocus genotype Data, Genetics, 155: 945-959

Rechinger K.H., 1982, Dracocephalum, In Flora Iranica, Rechinger K.H. (Ed.), 150. Akademische Druck-U. Verlagsanstalt, Graz, Austria, pp. 218-230

Richards C.M., 2000, Inbreeding depression and genetic rescue in a plantmetapopulation, American Naturalist, 155: 383-394
http://dx.doi.org/10.1086/303324

Setsuko S., Ishida K., Ueno S., Tsumura Y., Tomaru N., 2007, Population differentiation AND gene flow within a metapopulation of a threatened tree, Magnolia stellata (Magnoliaceae), American Journal of Botany 94: 128-136
http://dx.doi.org/10.3732/ajb.94.1.128

Sheidai M., Seif E., Nouroozi M., and Noormohammadi Z., 2012, Cytogenetic and molecular diversity of Cirsium arvense (Asteraceae) populations in Iran, Journal of Japanese Botany , 87: 193-205

Sheidai M., Zanganeh S., Haji-Ramezanali R., Nouroozi M., Noormohammadi Z., and Ghsemzadeh-Baraki S., 2013, Genetic diversity and population structure in four Cirsium (Asteraceae) species, Biologia, 68: 384-397
http://dx.doi.org/10.2478/s11756-013-0162-x

Sheidai M., Afsharb F., Keshavarzib M., Talebic S.M, Noormohammadid Z., Shafafa T., 2014, Genetic diversity and genome size variability in Linum austriacum (Lineaceae) populations, Biochemical Systematics and Ecology, 57: 20-26

Tero N., Aspi J., Siikamaki P., Jakalaniemi A., and Tuomi J., 2003, Genetic structure and gene flow in a metapopulation of an endangered plant species, Silene tatarica, Molecular Ecology, 12: 2073-2085
http://dx.doi.org/10.1046/j.1365-294X.2003.01898.x

Weising K., Nybom H., Wolff K., and Kahl G., 2005, DNA  Fingerprinting  in Plants. Principles, Methods, and Applications, (2^nd ed.), Boca Rayton, Fl., USA: CRC Press, pp. 472
http://dx.doi.org/10.1201/9781420040043

Zeng Q., Jin H.Z., Qin J.J., Fu J.J., Hu H.J., Li J.H., Yan L., Chen M., and Zhang D.W., 2010, Chemical constituents of plants from the genus Dracocephalum, Chemistry & Biodiversity, 7: 1911-1929
http://dx.doi.org/10.1002/cbdv.200900188