The influence of environmental heterogeneity on the morphological and genetic diversity of Circaea lutetiana (Onagraceae) in Hyrcanian forests

Environmental gradients are important factors that can potentially influence the genetic diversity and differentiation of populations. The present study examines the effect of environmental heterogeneity of the Hyrcanian forests on populations of Circaea lutetiana L. (Onagraceae). Using morphometrics, scanning electron microscopy (SEM) of leaf epidermis, and molecular markers, we analyzed genetic diversity and differentiation among nine populations from environmentally divergent habitats. Three different gene pools were observed. Our results indicate that the genetic structure is significantly correlated to environmental factors, but not to the geographical distance. Genetic clustering in C. lutetiana is affected by temperature, humidity, elevation, and average annual rainfall. Overall, our data indicate that gene flow does not contribute to explaining spatial patterns of genetic structure and the adaptation to the environment is the main factor shaping the genetic structure of the C. lutetiana populations. The populations belonging to each of the three gene pools have similarities in microclimate parameters, despite their geographical proximity, and populations from the same genetic pool can be up to 470 km apart. This evidence, as well as morphological and genetic similarities of the populations with greater geographical distance, suggest the possibility of cryptic speciation in this species.

controlled by gene flow, genetic drift and natural selection (Eckert & al. 2008). Gene flow can increase the genetic variability via introducing new genes into local populations (Hou & Lou 2011) and is often constrained by the distance between neighboring populations and the presence of landscape barriers. More closely situated populations tend to be more genetically similar (Slatkin 1993;Hutchison & Templeton 1999;Medrano & Herrera 2008). On the other hand, landscape barriers and environmental isolation of populations can limit gene flow (Duarte & al. 2015) and lead to low genetic diversity (Spehn & Körner 2005).
Plant species facing high levels of environmental heterogeneity within their distribution range need to develop adaptive strategies suited to their particular habitat including a variety of morphological, phenological and molecular changes (Ahmad & Prasad 2011;Wheeler & al. 2015;Little & al. 2016). Such adaptive genetic changes are often fixed within populations through natural selection resulting in genetic differentiation and diversity among populations from different environmental ranges. Temperature and moisture are among the main environmental factors that can affect genetic differentiation and diversity among populations (Still & al. 2005).
Based on recent phylogenetic studies, the genus Circaea L. (Onagraceae) consists of eight species and six subspecies distributed in Eurasia and North America (Xie & al. 2009). All Circaea species are tender, broad-leaved, shadeand moisture-loving herbs. The overall area of the genus closely coincides with the area of temperate mesophilic forests of the northern hemisphere (Wagner & Hoch 2009). Three characters, indehiscent fruit, a longer floral tube, and a pubescent stem, distinguish C. lutetiana from other species in the genus (Boufford 2005).
Two reproduction modes including rhizome and seed production are expected for Circaea lutetiana (Stearns 1989;Verburg & During 1998). Its life cycle shows the features of pseudo-annuals plants which survive the winter as seeds and hibernacles produced by the rhizome apices (Verburg & al. 1996;Verburg & During 1998). In the pseudoannuals plants, new rhizomes produced by shoot form a new ramet in the next growing season, not in the same year (Verburg & al. 1996). Moreover, the flowering of C. lutetiana occur from June to August (van der Meijden 1990) could be through self-compatibly, flowers diurnally, outcrossing, and pollinating by syrphid flies and small bees, or sometimes autogamously (Wagner & al. 2007).
Circaea lutetiana is a pluri-regional species distributed as a subcosmopolitan/cosmopolitan element of Hyrcanian forests of Iran, as well as, Europe, north of Africa and southwest of Asia (Akhani & al. 2010). The Hyrcanian forests ecological region is an area of lush lowland and montane forests extending from the south of Azerbaijan to about 900 km to the east into the Iranian provinces of Gilan, Mazandaran and Golestan. The occurrence of Hyrcanian forests between the northern slope of the Alborz Mountain and the southern edge of the Caspian Sea results in a high habitat heterogeneity (Akhani & al. 2010). Circa-ea lutetiana is distributed across wide range of microhabitats with remarkable heterogeneity regarding the density of the vegetation, altitude, temperature, sun exposure degree and humidity. Some studies over such environmentally heterogeneous microhabitats have suggested the influence of climatic conditions on plant genetic structure (Nevo 2001;Hamasha & al. 2013).
For the present study, we hypothesized that the effect of climatic and altitudinal heterogeneity of the environment would have a role in shaping patterns of macro/micro-morphological and genetic diversity among populations of Circaea lutetiana. Using morphometric, leaf epidermal scanning and two genetic markers, start codon targeted (SCoT) and direct amplification of minisatellite DNA (DAMD), we aim to answer the following questions: (i) Is there a significant relationship between morphological and genetic diversity of populations; (ii) how is genetic diversity distributed among and within populations of C. lutetiana; (iii) is genetic diversity related to geographical distance of population or to differences in climate conditions. The influence of gene flow and natural selection in shaping genetic diversity will be discussed.

Plant material and study area
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.

Morphological data
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).
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.

Data analysis
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.

Macro-and micro-morphological assay
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 correla-   tion 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 leafbase 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,    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 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, 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).

Genetic diversity
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.

Population structure and gene flow
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 accor-Gene flow between populations happens either by seed or pollen dispersal or both (Cain & al. 2000;Robertson & al. 2008). Given that seed dispersal is often found to be spatially restricted (Struik 1965), the gene flow into plant populations is mainly dependent on pollen dispersal (Bacles & Ennos 2008). We assume that environmental heterogeneity of the populations might have limited the pollen dispersal and consequently gene flow between populations. Apart of the effect of gene flow on the differentiation of populations, a trade-off between vegetative propagation and sexual reproduction, or both modes, could have existed in this pseudo-annual clonal plant (Verburg & During 1998).
The strong effects of environmental variation and microhabitat differences on the branches of rhizomes, seedling recruitment, size and number of hibernacles, genet size and weight could be expected (Verburg & al. 1996;Verburg & During 1998;Verburg & al. 2000). Verburg & al. (2000), in their study on clonal diversity in differently-aged populations of Circaea lutetiana using RAPD marker, found high genet diversity due to micro-site specialization. Our results indicate that in the populations with similar genetic structure but geographically distant, the strong selection to seedling recruitment by micro-habitat differences could be evidenced.
The environmental heterogeneity of habitats can cause intense variation in flowering time among adjacent populations depending on their altitude, amount of sun exposure, humidity, temperature, and vegetation density (Cortés & al. 2013;Duarte & al. 2015). Phenological variations in flowering set and flowering period can limit the pollen exchange and therefore decrease the degree of gene flow among populations (Franks & Weis 2009). Such kind of reproductive isolation may have contributed to the limitation of gene flow dance 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.

Environmental factors and genetic diversity
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.

DISCUSSION
According to the results of the Mantel and STRUC-TURE analyses, genetic differentiation among studied populations of Circaea lutetiana was not related to geographic distance. We could see a clear genetic split between neighboring populations, for example MA and LI, ZI and SA, HJ and SH, ZY and NA. Moreover, there are populations with similar genetic structure but geographically distant (LI and ZY; MA and ZI; HJ and LA). The absence of isolation-by-distance strongly indicates that gene flow between populations is infrequent (Hamasha & al. 2013). The NJ tree and Bayesian STRUCTURE analyses revealed that populations of Circaea lutetiana in the Hyrcanian region of Iran grouped in three different clusters. According to the CCA results, the genetic clustering of the populations is related to the similarity in microclimate conditions. The first cluster comprises the LI and ZY populations that are located far away towards the west and east of the Hyrcanian region. In the CCA results, both are separated from the rest of the populations due to their high elevation.
The second cluster comprises the SH, SA and NA populations, with a similar gene pool. According to CCA analysis, these populations are distinguished by their average annular rainfall that is the lowest among all populations. The four remained populations including LA, HJ, ZI and MA fall into the third cluster. Within this cluster, the ZI and MA populations are characterized by the relatively cool temperature and high level of humidity (along river bank), while the LA and HJ populations belong to the habitats with the highest average annual temperature. The two last populations show the most admixed genetic pool among the studied populations. Regarding the environmental features of the clusters, in disturbed areas the populations may be younger and seedling recruitment can increase the genetic diversity (i.e. LA, HJ, MA, but not ZY and ZI). In accordance with Verburg & al. (2000) study, the high genetic diversity in the LA, HJ, and MA populations may be due to habitat disturbance (caused by forest road) and recolonization of young ramets in empty patches. Based on Verburg & al. (2000) the empty patches due to disturbance could be increase genet diversity and genet size distributions.
All populations in deep, dark and dense forests (LI, NA, SA, SH) show low genetic diversity and maybe clonal reproduction by hibernacles of older genets effectively compete avoiding seedling recruitment. The phenotypic responses to various environments may also consist of highly special physiological, developmental as well as reproductive adaptations that improve plant function in those environments (Bradshaw 1965;Schmitt & al. 1999).
Based on Losos & Glor (2003) opinion, the formation of new taxonomic ranking depends on the morphological variation and geographical separation among population. Both morphological and genetic similarities of the two pair of populations including MA, ZI and LI, ZY may indicate that cryptic speciation under the influence of microclimate factors are occurring in different geographical locations.
Thus, the populations belonging to each of the three clusters have similarities in terms of microclimate parameters, despite their geographical proximity, and the population from the same genetic cluster can be up to 470 km apart (e.g. MA and ZI or LI and ZY). We therefore suggest that environmental heterogeneity has shaped the morphological and structural diversity of Circaea lutetiana in Hyrcanian forests. A similar study on Camisonia benitensis P.H.Raven (Onagraceae) with microsatellite markers, indicated evidence of cryptic genetic subdivision that did not correlate with habitat type, watershed, or physical distance between populations and populations from the same genetic cluster can be up to 29 km distant from one another (Dick & al. 2014).
It is well-known that the populations of a given species facing different environmental conditions may undergo genetic changes to adapt to their local conditions (Hufford & Mazer 2003) (Khan & al. 2015). Such adaptive genetic changes are often fixed within population through natural selection resulting in multiple, genetically distinct populations within a single species (Antonovics & Bradshaw 1970;Slatkin 1985;Mitton & al. 1997;Khan & al. 2015). There are a multitude of plant adaptations for tolerating or capitalizing on environmental heterogeneity. When the environmental heterogeneity is related to the availability of plant resource like water, light, and nutrients, it can directly affects rates of resource uptake, and thus growth, reproduction, and survivorship. On the other hand, heterogeneity in physical parameters like temperature, humidity, and wind speed may affects the plant function by modifying metabolic rates, stomatal control, or morphology, genet size and weight, phenology, flowering time, seedling recruitment, and thus resource processing (Fox & al. 2012;Verburg & During 1998;Verburg & al. 2000).

CONCLUSION
SCoT and DAMD analyses revealed significant genetic differentiation among the studied populations. However, the Mantel test showed no correlation between genetic distance and geographical distance. Moreover, with regard to macro/micro-morphological features, similar results were observed. Our study indicates that local adaptation to microhabitats played the most important role in the diversification of Circaea lutetiana populations in the Hyrcanian region of Iran. This may indicate the presence of cryptic species or subspecies in C. lutetiana, but a more thorough molecular study based on sequence markers is needed to confirm that conclusion. From a conservation perspective, given the positive and strong correlation between the genetic differentiation and environmental heterogeneity, the maintenance of populations in each region is crucial to ensure the preservation of genetic diversity in this species.