Forty-five plant accessions were collected from nine geographical populations belonging to three provinces in the north of Iran. Circaea lutetiana is a pseudo annual clonal plant (Verburg & al. 2000). We tried to collect individuals from different clones in order to represent different genets of the population. Details of localities are provided in Table 1 and Fig. 3. Voucher specimens are deposited in the Herbarium of Shahid Beheshti University (HSBU). Fresh leaves were collected and used for DNA extraction and molecular study. The study habitats show remarkable environmental heterogeneity regarding vegetation density, altitude, temperature, humidity and annual rainfall (Table 1). According to the vegetation density, three types of habitats could be distinguished including deep dark and dense forest, low-density forest, and very low-density forest along the forest road. The study area spans between 49°00’ to 54°47’ longitude and 35°10’ to 41°36’ latitude, and at elevations from 210 to 1700 m. The mean year temperature (AMT) ranges between 11.83-18.6 °C. The mean yearly rainfall (AP) ranges between 27.12-139.09 mm. In contrast with other populations that exhibit wet/sub-wet conditions, the Masoleh and Zirab populations grow on the edge of the river and have a lower mean year temperature (11-12°C). Information about environmental factors of the study area was provided from the State Meteorological Organization of Iran.

A total of 10 quantitative and 5 qualitative morphological characters were investigated. For morphological analyses, the mean of 10 individuals was used as representative of each population. The studied features included leaf-base shape, petal color, seed trichome type, leaf margin shape, stigma notch type, petal/sepal length, stamen/pistil length, petal notch/petal length, leaf/petiole length, length/width of leaf, stamen/petal length, length/width of seed, length of seed tip, and distance between teeth in the same leaf.

For the scanning electron microscopy (SEM) study, a small segment from the central part of the leaf abaxial surface was placed on aluminum stubs with double-sided cellophane tape and coated with gold. The specimens were examined with a Phillips × L20 SEM. UTHSCSA Image Tools Version 3.0 was used to measure quantitative and qualitative features including the shape of guard cells, peristomatal rim, ridge and striation patterns on the cuticular membrane, stomatal pole, length/width of stomata (L/W), stomatal rim location, outer surface ornamentation of guard cells were studied (Table 2). For the terminology of epidermal features, we follow Bartholtt & al. (1998) and Pole (2010).

CTAB-activated charcoal was used for total genomic DNA extraction (Ghasemzadeh Baraki & al. 2018). The quality of DNA samples was determined by the electrophoresis on 1% agarose gel. Four DAMD primers -viz., URP9F, URP2F, URP25F, URP38F- and four SCoT primers-viz., SCoT2, SCoT3, SCoT7, SCoT18- commercialized by UBC (the University of British Columbia) were used (Collard & Mackill 2009; Kang & al. 2002).

PCR reactions for SCoT and DAMD markers were carried out in a 25 μl volume containing 10 mM Tris-HCl buffer at pH 8.3, 2.5 mM MgCl2, 1 mM dNTP mix (Cinna Gen Co, Iran), 0.2 μM of a single primer, 15-40 ng of template DNA, and 1 U of Taq DNA polymerase (Cinna Gen Co, Iran). The amplification reactions for both assays were performed with a T100 thermocycler (BIORAD, USA).

For DAMD assay, all amplification were carried out as follows: 94°C for 5 min, followed by 40 cycles of denaturation at 94°C for 1 min, annealing at 55°C for 1 min, and extension at 72°C for 90 sec. A final extension cycle at 72°C for 10 min was followed. The amplification reactions of SCoT were performed with the following program: 5 min initial denaturation at 94°C followed by 36 cycles of 1 min at 94°C, annealing at 57.5°C for 1 min, and 90 s at 72°C. The reaction was completed by a final extension step of 10 min at 72°C.

The amplification products were visualised by electrophoresis in 1% (w/v) agarose gels, followed by 0.5 µg/ml powerload staining. Fragment sizes were estimated using a 100 bp size ladder (Fermentas, Germany), photographed under UV light and scored for the presence or absence of bands.

The morphological features were standardized (Mean = 0, Variance = 1) and used to establish the Euclidean distance among taxa pairs (Podani 2000). The obtained distances were used for clustering. The Ward method of hierarchical cluster analysis (Podani 2000) was used for grouping the plant specimens based on the Euclidean distance. Alongside, PCA (principal components analysis) was carried out to identify the most variable morphological features (Podani 2000). PAST v. 2.17 (Hammer & al. 2001) was used for the above analyses.

The obtained SCoT and DAMD bands were coded as binary characters (presence = 1, absence = 0). The genetic diversity parameters like allele diversity (Weising & al. 2005), Nei’s gene diversity (H), Shannon information index (I), the number of effective alleles, and percentage of polymorphism (Freeland & al. 2011) were determined using the GenAlex v. 6.4 software in each population for a dominant molecular marker (Peakall & Smouse 2006). Nei’s genetic distance was used for clustering (Weising & al. 2005; Freeland & al. 2011). Neighbour Joining (NJ) clustering was used for grouping (Freeland & al. 2011) after 100 times bootstrapping. The NJ tree with the Dice coefficient showed a better grouping than one done with the Jaccard coefficient. PAST and GenAlex were used for these analyses. Genetic differentiation of the studied specimens and populations was studied by AMOVA (analysis of molecular variance) test (with 1000 permutations) as performed in GenAlex v. 6.4 (Peakall & Smouse 2006).

Genetic admixture among populations was analyzed by model-based clustering as performed by STRUCTURE v. 2.3 (Pritchard & al. 2000). We used the admixture ancestry model under the correlated allele frequency model. Data were scored as dominant markers and analysis followed the method suggested by (Falush & al. 2007). Indirect evaluation of the gene flow (Whitlock & McCauley 1999) among populations was checked by reticulation analysis in the DARwin v. 5.0 (Perrier & Jacquemoud-Collet 2006) with a reticulation tree.

A Mantel test was performed to check for correlation between geographical and genetic distances of the studied populations (Podani 2000). PAST v. 2.17 (Hammer & al. 2001) and Genealex v. 6.4 (Peakall & Smouse 2006) were used for these analyses. The Pearson coefficient of correlation was determined between geographical features (longitude and latitude) and genetic diversity parameters. A canonical correspondence analysis (CCA) was done using PAST v. 2.17c (Hammer & al. 2001), to determine the relative importance of geographical factors in the spatial organization of genetic diversity among genotypes.

In total 23 macro (see material methods section) and micro-morphological characteristics of Circaea lutetiana populations were studied (see leaf SEM features in Table 2). Among them 17 features were selected based on PCA-biplot (Fig. 1). PCA analysis of morphological features revealed that the first three PCA components comprised about 70.38% of the total variability of the studied populations. In the first PCA axis with about 33% of the total variation, morphological traits like the stomatal pole, leaf-base shape, petal color possessed the highest correlation while in the second PCA axis, characters like seed trichomes, ridge pattern of the cuticular membrane, and stomatal rim location; in the third PCA axis characters like leaf margin and stigma notch possessed the highest correlation. Therefore, these morphological characters were the most variable morphological characteristics among the studied populations. PCA-biplot of morphological characters (Fig.1) separated the populations studied into distinct groups. It also showed that the HJ specimen was differentiated from the other populations due to its seed trichomes, while the SA, ZI, and MA specimens were differentiated from the others due to the “T-pieces” feature of the stomatal pole and petal color (white-pink to pink) (Fig. 3).

Based on macro and micro-morphological traits, the studied specimens were separated in the Ward tree of morphological features (Fig. 2). Two main branches resulted. One included four populations arranged in two sub-groups: subgroup I, including SA and ZI, subgroup II including MA and NA. The rest of the populations were placed in another branch, which also had two sub groups: populations HJ became detached from others and placed in one sub-branch. As seen in the PCA biplot (Fig. 1), each of the studied populations showed distinct morphological trait(s). The leaf-base shape was truncate to truncate-cuneate in the ZI, SA, and NA populations, while the others had a round-truncate leaf base. Petal color was white-pink to pink in HJ, ZI, SA, MA, and white in the other populations. Considering the measurement of denticulation of leaf margin, there were three groups: slightly denticulate (HJ, MA, NA), denticulate (LI, LA, SH, ZY), and heavily denticulate (ZI, SA). The notching on stigma was deeper in ZY, LA, SH, and lower in the other populations. This classification does not agree with geographical distance. Although morphological variation was observed within the population, we tried to consider the average of traits in the population. In the studied specimens, there were two types of stigma knotch including deeply (in ZY, LA, SH) and shallowly (Fig. 3h, 3f, 3e) knotched. In all populations, except for populations HJ, MA, ZI, and SA, whose petals were white-pink or pink in color, the petals were completely white (Fig. 3i, 3c, 3b, 3a). As seen in Fig. 3 leaf-margin of the specimens showed three character states including slightly denticulate (in HJ, MA, NA; Fig. 3i, 3c, 3d), denticulate (LA, LI, ZY, SH; Fig. 3f, 3g, 3h, 3e) to heavily denticulate (ZI, SA; Fig. 3b, 3a). The following two main types of epidermis were recognized. Type I (Ia: SA, ZI and Ib: SH, HJ) had a fine striation on the outer surface of guard cell and slight to normal ridged on the cuticular membrane. In type I the stomatal rim was located nearly to the epidermal surface (Fig. 3a, 3b, 3e, 3i). Type II (IIa: NA, ZY and IIb: MA, LI, LA) had normal striation on outer surface of guard cells and rarely exhibited striations on the cuticular membrane. In this type the stomatal rim was raised in relation to the surface of the epidermis (Fig. 3d, 3h, 3c, 3g, 3f).

The amplification of primers 4 SCoT and 4 DAMD resulted in 148 bands. The genetic diversity parameters related to the studied populations are shown in Table 3. The highest values for the mean number of effective alleles (Ne = 1.36, 1.31) and gene diversity (He = 0.2, 0.18) were in populations LA and SA, respectively. These populations also showed the highest values for the genetic polymorphism percentage (52.44, 47.54). The lowest values for the mean number of the effective alleles and gene diversity occurred in LI (Ne = 1.23, He = 0.12) and SH (Ne = 1.24, He = 0.13). These populations showed the lowest values for the genetic polymorphism percentage (28.38, 33.78). The mean number of different alleles over all loci (Na) ranged between 0.62 (LI)-1.1 (LA). Shannon’s information index ranged between 0.17 (LI)-0.29 (LA).

The AMOVA results indicated significant molecular difference (FPT = 0.46, P = 0.01) among populations indicating a great level of genetic differentiation. It also revealed that 46% of the total genetic diversity was related to inter-population differences, while 54% was related to intra-population differences.

The Dice distances between pairs of populations were calculated based on the 148 analysed bands. The Neighbor Joining (NJ) dendrogram (Fig. 4) consisted of two major clusters each containing two subclusters. Populations LA, HJ nested within subcluster 1 and MA and ZI occurred in subcluster 2 of first main cluster. The members of populations ZY and LI formed subcluster 1, while populations SH, SA, NA occurred in subcluster 2 of the second major cluster.

The Mantel test revealed no significant correlation (r = 0.02, p = 0.5) between genetic distance and geographical distance (Fig. 5). According to the Evanno’s method, the STRUCTURE analysis indicated K = 3 as the most likely number of gene pools. A STRUCTURE plot based on k = 3 is shown in Fig. 6. In general, populations SA, NA and SH were genetically more alike and form the first genetically similar group. The same holds true for populations LI, ZY which compose the second group. Populations LA, HJ, MA, and ZI formed the third genetic group. Some of the populations were highly genetically discrete (SA, SH, NA ZY, LI, ZI), while the others showed a high degree of genetic admixture due to common ancestral alleles (MA, LA, HJ). This Bayesian approach analysis was in accordance with NJ tree and PCA plot. The reticulation method of networking agreed with the least square method and determined similar/shared alleles between populations. Also, the reticulogram plot (Fig. 7) indicated a limited degree of gene flow among populations HJ-MA, LA-LI, and ZI-ZY. None of these three pairs are located at close geographical distance between each other.

The CCA analysis was used to investigate the influence of environmental factors on the differentiation of populations. According to the CCA plot (Fig. 8), the environmental features of the different habitats affected population’s aggregation (Axis 1 = 45.18 % of variance; Axis 2 = 31.78 %). The clustering of LA, HJ, NA, SH, and SA populations was influenced by mean annual temperature (AMT). Within this group mean annual rainfall (AP) separated SH and SA populations from the others. The clustering of ZI and MA was influenced by humidity. Elevation influenced the grouping of LI and ZY.