Genetic Diversity and Structure of Dalmatian Pyrethrum (Tanacetum cinerariifolium Trevir. /Sch./ Bip., Asteraceae) within the Balkan Refugium

Dalmatian pyrethrum (Tanacetum cinerariifolium Trevir. /Sch./ Bip.) is an outcrossing, perennial insecticidal plant, restricted to the eastern Adriatic coast (Mediterranean). Amplified fragment-length polymorphisms (AFLP) were used to investigate the genetic diversity and structure within and among 20 natural plant populations. The highest level of gene diversity, the number of private alleles and the frequency down-weighted marker values (DW) were found in northern Adriatic populations and gradually decreased towards the southern boundary of the species range. Genetic impoverishment of these southern populations is most likely the result of human-related activities. An analysis of molecular variance (AMOVA) indicated that most of the genetic diversity was attributed to differences among individuals within populations (85.78%), which are expected due to the outcrossing nature of the species. A Bayesian analysis of the population structure identified two dominant genetic clusters. A spatial analysis of the genetic diversity indicated that 5.6% of the genetic differentiation resulted from isolation by distance (IBD), while 12.3% of the genetic differentiation among populations followed the pattern of isolation by environmental distance (IBED). Knowledge of the genetic diversity patterns of the natural populations and the mechanism behind these patterns is required for the exploitation and possible conservation management of this endemic and economically important species.


Introduction
Conservation of genetic diversity is essential for any species long-term survival [1], providing the ability to respond to selective pressures caused by both anthropogenic and natural environmental changes [2]. Patterns of genetic diversity and population structure are specifically shaped in response to many factors and are result of both past history and present processes [3]. Primarily, climatic oscillations during the Quaternary causing repeated restriction-recolonization events have left distinguishing signatures in modern species genetics [4]. These oscillations have had a considerable impact on the current distribution and genetic diversity in numerous animal and plant species [5][6][7][8] in which altitudinal and latitudinal range shifts were a likely response to intense climatic changes [9][10]. Detailed phylogeographical studies have revealed that three European peninsulas of the Mediterranean area (Apennines, Balkans and Iberia) have served as refugia for many temperate species and a postglacial recolonization source for central and northern Europe [3,7,[9][10][11]. However, this traditional interpretation of southern European peninsulas as being single glacial refugia has been modified, and the model ''refugia within refugia'' was developed [10,[12][13][14][15]. Genetic divergence expressed as a rarity index (DW; frequency down-weighted marker values) as an indicator of the historical processes of a species [16] enables the distinction between old vicariance and recent dispersal. As a consequence of long-term isolation, populations in glacial refugia are characterized by high levels of DW values and genetic diversity [17], whereas populations that were established after the most recent glacial cycle show little genetic variation [18] as an outcome of successive founder events during postglacial colonization [6,9]. However, zones of secondary contact of divergent lineages can also exhibit high levels of genetic diversity [11,19].
The Balkan Peninsula itself represents a unique ecological and biogeographical phenomenon in Europe possessing the tremendous number of plant and animal species and pronounced richness in endemic taxa as well. This remarkable diversity is the result of the variety of its regions, complex geological history, and interactions between populations, species and ecosystems [20]. In the last decade various phylogeographic and molecular studies have been conducted in this area with the aim of elucidating the genetic relationships of its flora and fauna [10,[21][22][23][24][25].
It is widely acknowledged that geographical and environmental factors can also play significant roles in the modeling of speciesspecific inter-and intra-population genetic structure. Spatially separated populations may experience isolation by distance (IBD) phenomena in which geographical barriers and physical distances cause restricted pollen and seed dispersal, i.e., gene flow [26][27][28][29][30]. Greater genetic divergence between populations inhabiting different environments suggests that various present and/or past ecological conditions may significantly impact the genetic differentiation of local populations [31][32], forming an isolation by environmental distance (IBED) pattern [33]. A significant positive partial correlation, after eliminating the effect of geographic distance, indicates the contribution of climatic gradients to patterns of genetic divergence [33][34].
Nevertheless, species are greatly affected by direct or indirect human activities including population fragmentation and overharvesting pressures, all resulting in loss of genetic diversity and the alteration of their population genetic structure [35]. Moreover, founder effects and conscious or unconscious human selections for high yield and quality through cultivation of wild plants likewise impose profound impacts on the genetic diversity patterns by narrowing genetic base of the plant material [36][37][38].
For our investigation we have chosen Dalmatian pyrethrum (Tanacetum cinerariifolium Trevir. /Sch./ Bip.), an insecticidal perennial plant. Dalmatian pyrethrum is an outcrossing diploid (2n = 18) [39], self-incompatible and thermophytic plant species of the family Asteraceae [40]. This species is restricted to the eastern coast of the Adriatic Sea, and its distribution range expands from the coastal region to 200 m.a.s.l., but in some cases, it is also present in higher mountainous Mediterranean zones at elevation above 500 m.a.s.l. Continuously distributed populations of Dalmatian pyrethrum can be found in extremely degraded habitats with shallow rocky soils. In Croatia, abundant populations can be found in the southern parts of the Istrian Peninsula (Premantura), the Kvarner islands (Krk, Cres and Lošinj), the Velebit and Biokovo mountains and along the Dalmatian coastal region and its islands (Brač, Hvar, Biševo, Vis, Korčula, Lastovo and Mljet) [41]. Apart from Croatia distribution range expands to the southern parts of Bosnia and Herzegovina, and coastal regions of Montenegro and Albania [42]. Plants produce the natural insecticide pyrethrin, which attacks the nervous systems of the insects causing a knock-down effect and death [43].
The cultivation and use of Dalmatian pyrethrum and its products have a documented history in Croatia. The dust from grinding dried flowers has traditionally been used in Croatian households and agriculture [44]. Originally, flowers were excessively gathered from the wild resources; however, due to the high demand for plant material, cultivation began in 1850s near Dubrovnik and rapidly expanded along the Dalmatian coastal region and the islands [45]. The first fields were established with seeds gathered from the wild. Dalmatian pyrethrum was primarily reproduced through vegetative splits and later by directly sowing the seeds [46]. It was grown in the middle and southern parts of Dalmatia, including Zadar, Biograd, Š ibenik, Split, Hvar, Brač, Makarska, Korčula and Dubrovnik, as well as the northern Adriatic islands Krk and Cres [44]. In 1914, Croatia became the leading producer with high quality plant material. The highest yield was achieved in 1926 with approximately 1,359 tons of dried pyrethrum flowers [46].
The goal of this study was to analyze the genetic diversity and structure within and among 20 populations of Dalmatian pyrethrum across their geographical range in Croatia and to determine impact of (i) Pleistocene climatic oscillations, (ii) geographic (IBD) and (iii) environmental isolation (IBED), and (iiii) human activities through overexploitation and cultivation on the obtained pattern of genetic diversity. AFLP markers have found the widest application in assessing the variability and genetic structure of populations in a variety of plants [47][48][49][50][51][52][53] and once again have proven to be invaluable tools in describing the former. Despite its historical significance and value, the genetic diversity and structure of Dalmatian pyrethrum has never been investigated; therefore, this study using AFLP markers represents the first such attempt. Dalmatian pyrethrum is nationally strictly protected by The Croatian Ordinance on the proclamation of protected and strictly protected wild taxa (Official Gazette 99/09; http:// narodne-novine.nn.hr/clanci/sluzbeni/dodatni/403046.pdf). Today, as an endemic species, Dalmatian pyrethrum is threatened by anthropogenic habitat loss and degradation related to factors such as urbanization and habitat conversion; therefore, the need and importance of these species' genetic resources protection must be recognized even more. The gathered information could be of utmost importance for future activities regarding the conservation and breeding programs of this endemic and economically important species.

Ethics statement
Collecting of all samples from natural habitats in Croatia including Nature Park Biokovo was permitted by the authority of Directorate for Nature Protection of the Ministry of Culture of the Republic of Croatia (532-08-01/6-06-02; 07.07.2006).

Sampling of plant populations
A total of 411 specimens of T. cinerariifolium were sampled from 20 locations along the Adriatic coast in Croatia ( Figure 1). Between 13 and 25 samples per location were collected as fresh leaf tissue ( Table 1). The vast majority of the species populations is distributed in Croatia, therefore we focused our study exclusively on a wide sampling of Croatian populations.

DNA isolation and AFLP analysis
The total genomic DNA was isolated from 90 to 110 mg of fresh leaves tissue with DNA-GenElute Plant Genomic DNA Miniprep Kit (Sigma Aldrich, Steinheim, Germany), according to manufacturer's instructions. Prior to isolation leaves tissue was quickly frozen in liquid nitrogen (2196uC) and ground to a fine powder with TissueLyzer homogenizer set (Qiagene, Hilden, Germany) at 30 Hz for 1 minute. The DNA concentrations were measured using a Qubit Fluorometer (Invitrogen, Carlsbad, California). The AFLP analysis was carried out following the original protocol that was described by Vos et al. [54] with some modifications [55]. Restriction of the genomic DNA was performed using the restriction enzymes EcoRI and MseI. Specific adapters were ligated using T4 DNA Ligase (Fermentas, Thermo Scientific). The restriction digestion and adapter ligation, preamplification and selective amplifications were performed in a GeneAmp PCR System 9600 (Applied Biosystems, Carlsbad, SAD). Six primer combinations were selected for amplification; VIC-EcoRI-ACGC+MseI-AGA, NED-EcoRI-AGA+MseI-CAG, VIC-EcoRI-ACG+MseI-CGA, FAM-EcoRI-ACC+MseI-CGA, NED-EcoRI-AGA+MseI-CGA, PET-EcoRI-ATG+MseI-CGA. The amplified fragments were separated by capillary electrophoresis in an ABI3130xl Genetic Analyzer (Applied Biosystems, Carlsbad, USA). The amplified DNA fragments in the size range of 50 to 125 bp were scored with GeneMapper 4.0 software (Applied Biosystems, Carlsbad, USA) as being present (1) or absent (0) in order to create a binary matrix.

Data analysis
The selection of markers for further statistical analysis and the estimation of the error rate per primer combination were performed as suggested by Herrmann et al. [56]. Thirteen DNA samples were used as duplicated samples, two of which were chosen for multiple controls, while 14 samples represented negative controls. The output from GeneMapper (fragment size and peak heights) was imported into the marker selection algorithm scanAFLP 1.3r (www-leca.ujf-grenoble.fr/logiciels.htm).
The genetic diversity within populations was estimated by determining the percentage of polymorphic bands (%P), the number of private alleles (N pr ) and Shannon's information index (I) [57][58]. Shannon's information index was used to measure the total diversity (H T ) as well as the average intra-population diversity (H P ) and the diversity within (H P /H T ) and among populations [(H T -H P )/H T ]. As a measure of divergence, the frequency downweighted marker values (DW) [23] were calculated using the software AFLPdat [59].
The marker allelic frequencies were estimated from the observed frequencies of the fragments using a Bayesian method with non-uniform prior distribution of allele frequencies proposed by Zhivotovsky [60] as implemented in AFLP-Surv ver. 1.0 [61].
This approach assumes Hardy-Weinberg equilibrium (F IS = 0) justified by the outcrossing nature of Dalmatian pyrethrum. The obtained allelic frequencies were used in the analysis of genetic diversity within and between populations following the treatment of Lynch and Milligan [62]. The population genetic structure was described in terms of the total gene diversity (H T ), the average gene diversity within populations (H W ), the average gene diversity among populations in excess of that which was observed within populations (H B ), and Wright's F ST statistics. The expected heterozygosity or Nei's gene diversity (H E ) [63] was calculated for each population as well as pairwise Nei's genetic distances (D NEI72 ) [62,64]. The unrooted tree was constructed using the Fitch-Margoliash method in the PHYLIP ver. 3.6b software package [65]. A total of 1,000 distance matrices were computed by bootstrapping using the program AFLP-Surv, and the bootstrap values were calculated from FITCH and CONSENSE within PHYLIP.
An analysis of molecular variance (AMOVA) [66] was used to partition the total genetic variance within and among populations using Arlequin ver. 2.000 [67]. The variance components were tested statistically by non-parametric randomization tests using 10,000 permutations. Pairwise population comparisons that were examined with AMOVA resulted in Q ST values that were equivalent to the proportion of the total variance that was partitioned between the two populations and could be interpreted  as the inter-population distance average between any two populations [68]. The genetic structure was assessed using the model-based clustering approach implemented by the software STRUCTURE 2.3.3 [69], which simultaneously identifies clusters and assigns individuals to clusters using a Bayesian approach. Ten runs of STRUCTURE were performed by setting the number of clusters from 1 to 21 (K = 1-21). Each run consisted of a burn-in period of 200,000 steps followed by 10 6 MCMC (Monte Carlo Markov Chain) replicates assuming an admixture model and correlated allele frequencies. No prior information was used to define the clusters. The calculations were carried out at the Bioportal of the University of Oslo (http://www.bioportal.uio.no/). The most likely number of clusters (K) was selected by calculating an ad hoc statistic DK based on the rate of change in the log probability of the data between successive K values as described by Evanno et al. [70] and implemented in Structure-sum-2011 [71]. The runs with the maximum likelihood were selected, and by averaging the estimated membership coefficients of the individuals, the proportion of ancestry of each population in each of the clusters was calculated.
A further population mixture analysis was conducted using BAPS ver. 5.3 [72] (a) without using the geographic origin of the samples as an informative prior ('clustering of individuals') and (b) with this prior ('spatial clustering of individuals') [73]. BAPS was run with the maximal number of clusters (K) set to 20, and each run was replicated 10 times. The results of the mixture analysis were used as input for population admixture analysis [74] with the default settings to detect admixture between the clusters.
Isolation by distance (IBD) among the populations was tested using the method of Rousset [75]. Mantel's test [76] was used to compute and test the linear correlation between the matrix of the natural logarithm of geographic distances (in km) between pairs of populations and the matrix of pairwise F ST /(1-F ST ) ratios.
The plausibility of patterns of isolation by environmental distance (IBED) was evaluated, where differences in the bioclimatic variables among the sampling sites would better explain the genetic distance patterns [33,77]. The climate data for 25 sampling sites were obtained from the WorldClim database at a spatial resolution of approximately 1 km 2 [78]. The ecological characteristics of the sampling sites were described using 19 bioclimatic variables (11 temperature-and 8 precipitation-related) representing the annual trends, seasonal variations and extremes in temperature and precipitation (Table 2). Bioclimatic variables were used for the Principal Component Analysis (PCA) and subsequently to calculate the environmental distances among populations based on principal component scores rather than original values to avoid intercorrelation among variables. PCA was performed using the PROC PRINCOMP procedure in SAS. The number of principal components was determined by checking the eigenvalues of the principal components and the cumulative proportion of the explained variance. The biplot was constructed using two principal components showing populations and bioclimatic variables (as vectors). The standardized scores of the first three principal components (eigenvalue .1) were multiplied by the root of their eigenvalues, and the pairwise Euclidean distance between populations were calculated to obtain an environmental distance matrix.
A simple Mantel's test was used to test the linear correlation between the matrix of environmental and genetic [F ST /(1-F ST )] distances, and a three-way Mantel's test was conducted between the same distances matrices while considering the geographical distances among populations. Thus, the correlations between the residual environmental and genetic distances among the analyzed populations represent isolation by environmental distance (IBED) pattern. The significance levels were obtained after 10,000 permutations as implemented in NTSYS-pc ver. 2.10s [79].

ScanAFLP
The initial number of polymorphic bands that were obtained with six selective primer combinations equaled 850, while after the applying scanAFLP, the remaining dataset consisted of 593 polymorphic markers (43.35% of the initial number of polymorphic markers). The estimated error rate per primer combination ranged from 0.68% to 3.50% with a mean of 1.83%, which is within the typical range found in AFLP datasets [80].  (Table 1, Figure 2B).

Genetic structure among populations
A moderate level of population differentiation was observed (F ST = 0.078) ranging from 0.018 between pops. 10 and 11 to 0.161 between 6 and 18.
The AMOVA analysis showed that although the most of the genetic diversity was attributable to differences among individuals within populations (85.78%), the Q ST values among populations were highly significant (Q ST = 0.141 p,0.0001) confirming the existence of population differentiation. The Q ST values ranged from 0.049 between pops. 4 and 11 to 0.278 between pops. 5 and 16.
Nei's distance (D NEI72 ) varied from 0.002 between pops. 10 and 11 to 0.022 between pops. 6 and 18 with an average value of 0.010. The pattern of grouping in the Fitch-Margoliash tree suggests a direct relationship between the populations and their geographical origins. This pattern revealed the existence of two clades (northern vs. middle and southern Adriatics) supported by a bootstrap value of 88%. Inside the northern clade, the inland pops. 5 and 6 along with the population 2 separated from the northern island populations with a bootstrap value of 72%. The pops. 9 and 11 separated from pops. 8 and 4 with a bootstrap value of 88%. The populations from Mt. Biokovo (pops. 16, 17 and 18) grouped together and showed a clear separation from the remaining southern populations with a high bootstrap support of 99% ( Figure 3).
The results of the STRUCTURE analysis are reported in Figure 2C and 2D. The average estimates of the likelihood of the data, conditional on a given number of clusters, [ln P(X|K)], were obtained for each of the 10 independent runs for each K (from K = 1 to 21). The estimates of [ln P(X|K] increased progressively as K increased. Using these estimates, the highest DK value was observed for K = 2 (302.59). The second highest DK value was Table 2. Correlations between 19 environmental variables (BIO01-BIO19) and the first three principal components. Thus, the hypothesis of having more than three clusters (gene pools) was far from optimal. The proportions of membership (ancestry) of each individual in each of the gene pools were calculated for K = 2 to 3 based on the run with the highest [ln P(X|K)], and populations were assigned to a particular gene pool (A, B or C). The results were in accordance with the results that were obtained using the distancebased method. Representative populations of each gene pool that were identified in the STRUCTURE analysis were located in proximity of each other. At K = 2, the representative populations of gene pool A were mainly those from the northern Adriatic, while the populations from the middle and southern Adriatic were representative of a separate gene pool B. The most admixed populations were 9 and 11 with almost equal proportions of ancestry from the two gene pools ( Figure 2C). At K = 3, a newly formed gene pool C was identified and was predominant in the pops. 16, 17 and 18 from Mt. Biokovo ( Figure 2D).
The BAPS mixture with or without spatially informative priors resulted in the congruent assignment of individuals to two clusters. The best partitions received log marginal likelihoods of 2 51646.2626 at P = 1 (without using geographical coordinates as informative priors) and 251793.7943 at P = 1 (spatial clustering). The results were consistent with the results of the STRUCTURE analysis ( Figure S1). Isolation by distance (IBD) and isolation by environmental distance (IBED) A significant correlation was found between the genetic and geographical distances (r = 0.236; P Mantel = 0.012), suggesting that 5.6% of the genetic differentiation between populations could be explained by geographical distances ( Figure 4A).
A Principal Component Analysis based on the correlation matrix revealed that first three principal components had an eigenvalue .1 and explained more than 90% of the total variation (Table 2). Eight temperature-related variables (BIO01, BIO03, BIO05, BIO06, BIO08, BIO09, BIO10 and BIO11) and one precipitation-related variable (BIO15) were highly positively correlated (r.0.70; P,0.001) with the first principal component (PC1), while four precipitation-related variables (BIO12, BIO14, BIO17 and BIO18) were highly negatively correlated (r,20.70; P,0.001). Strong positive correlations (r.0.70; P,0.001) were observed between the second principal component (PC2) and two precipitation-related variables (BIO13 and BIO16). The biplot that was constructed by the two principal components showing populations and 19 environmental factors (as vectors) is presented in Figure 5. PC1, explaining 60.79% of the total variation, separated the middle (pops. 7-11) and southern (pops. 13-15, 19 and 20) Adriatic populations occupying warmer habitats from pops. 12, 16, 17 and 18, which were sampled at higher elevation (516 to 1,335 m.a.s.l.). Along PC2, explaining 21.04% of the total variation, the northern Adriatic populations (pops. 1-6), occupying habitats that were characterized by greater rainfall, were separated from the remaining populations. The average Euclidean distance between the pairs of populations that was calculated based on the standardized scores of the first three principal components was 5.34, ranging from 0.08 (between pops. 16 and 17) to 10.12 (between pops. 8 and 16, respectively).
A greater correlation was found between the genetic and environmental distances (r = 0.376; P Mantel ,0.001) than between the genetic and geographical distances ( Figure 4B). After considering the effect of the geographical distances in a three-way Mantel test, the correlation between the residual genetic and residual environmental distances remained significant (r = 0.351; P Mantel = 0.005) suggesting that isolation by environmental distance (IBED) explains 12.3% of the observed genetic pattern ( Figure 4C). The removal of the effect of environmental variation resulted in a still significant but much lower correlation between the genetic and geographic distances (r = 0.189; P Mantel = 0.039).

Discussion
Identification of genetic diversity patterns is of fundamental importance for effective management and species conservation   [7], and as we expand our knowledge about the genetic diversity and structure of natural populations, it subsequently becomes important to understand the mechanisms behind the determined patterns. Therefore, the genetic diversity and population genetic structure of Dalmatian pyrethrum must be analyzed both in the context of demographic history and biological traits of the species [3]. Additionally, understanding the impact of geographical distances and consequently restricted gene flow among the analyzed populations as well as climatic characteristics of the sampling sites can contribute to a better comprehension of the relationships between the populations. Furthermore, the fact that human historical activities might have influenced the observed pattern must not be excluded.
The results of our study show genetic differentiation of T. cinerariifolium into two ( Figure 2C) to three ( Figure 2D) distinct groups. At K = 2 it is plausible that the north-south split in the study area is due to widely spaced canyon of karstic river Zrmanja in northern Dalmatia as karst rivers are often major barriers to the otherwise continuous distribution of species. The obtained results are somewhat similar to the phylogeographical patterns detected in Edraianthus tenuifolius [10], Campanula pyramidalis [25] and Cardamine maritima [81] showing either phylogeographical or taxonomical split in the area of river Neretva valley (central Dalmatia). However, the fact that north-south split might be the artefact of genetic impoverishment (discussed further in the text) of central and southern population must not be neglected. With additional splitting of southern gene pool at K = 3, a new pattern that is congruent with the survival of the species in multiple refugial areas, followed by a period (s) of isolation occurs. Three genetically and geographically distinct and well-defined groups of Dalmatian pyrethrum populations therefore support the ''refugia within refugia'' model that was already confirmed in numerous studies in the Balkan Peninsula [10,21,23,24]. The restricted gene flow among inter-refugial populations has led to genetic differentiation, resulting in well-defined gene pools [82]. Among the 20 populations that were included in this study, the highest levels of gene diversity, private alleles and DW were found in group of northern Adriatic populations from Kvarner bay. The genetic structure of these populations, which are grouped together within the same cluster as revealed by the STRUCTURE analysis, suggests their ancient origin and significant isolation through a recent period of evolution [17] Such northerly located microrefugia for a thermophytic species like Dalmatian pyrethrum is quite unexpected but is strongly supported by the obtained results. Previous research in this region on Dalmatian wall lizard (Podarcis melisellensis) [21], nose-horned viper (Vipera ammodytes) [83] and Edraianthus tenuifolius [10] found no evidence of a refugium that far north. Our results could possibly be explained by numerous glaciations and interglaciation periods that are characteristic of the Quaternary and that were accompanied by vast species migrations. During the last interglacial period (approximately 130,000 yr. B. P. to 116,000 yr. B. P.), the global climate was warmer than at present [84][85], and it is possible that the populations that were initially located in the more southern and thus warmer parts of Adriatic basin expanded far to the north-west. By tracking suitable microhabitats some populations survived during the cold period, and isolation from their pool of origin eventually has led to strong genetic differentiation that can now be detected through northern Adriatic gene pool. These microhabitats in Kvarner bay served as microrefugia, and at the end of the last glacial period, these protected and long-lived populations recolonized the northern Adriatic regions. However, a major flaw of this hypothesis is that in this scenario, there should be at least some evidence of ancient populations in the southern parts of the Adriatic basin. Generally, in more southern populations, the levels of gene diversity, number of private alleles and DW decrease. Overall, the detected pattern is not in concordance with the general predictions that the highest levels of diversity can be detected among southern regions of species distribution range [5,82], thus suggesting either a more recent origin of these populations or some past event which resulted in abrupt changes of their genetic structure. At this point, human activity and the excessive gathering of this species during the last century, primarily in the southern and, to a lesser extent, in the northern parts of Adriatic, must be considered. In fact, until the 1850s, the only reliable source of pyrethrum material was from wild plants [45]; hence, it is very likely that overexploitation alone is responsible for such low gene diversity and DW. An exception is the population 19 located in the southern parts of the species' range. Because it exhibits high levels of gene diversity accompanied by low levels of private alleles and DW, this population is likely the zone of secondary contact between the divergent lineages and is anthropogenic in origin. Additionally, some populations show a quite admixed pattern comprising large fractions of different gene pools (pops. 8, 9, 11, etc.). The specific structure of these populations could also easily be a consequence of the admixture of divergent lineages through secondary contact after postglacial colonization, however, more likely explanation of the obtained results could lie in the cultivation practices of Dalmatian pyrethrum, which began in 1850s near Dubrovnik (Southern Dalmatia) and rapidly expanded along the Dalmatian coastal region and the islands [45]. These populations could have been founded through human activity by the exchange of seed material between farmers from different geographical regions. It is well documented that anthropogenic activities foster genetic homogenization and decrease genetic differentiation among populations [36,[86][87]. Therefore, we hypnotize that for southern populations overexploitation and process of cultivation are responsible for substantial loss of private alleles, levels of gene diversity and DW. Moreover, extensive seed exchange between different populations caused admixture as observed in pops. 8, 9 and 11. Thus, this phylogeographical study could represent typical example of humans' shaping the genetic structure of investigated plant species. The genetic impoverishment as a result of either overexploitation of wild resources or cultivation practices has also been demonstrated for American ginseng (Panax quinquefolius L.) [35], Rheum tanguticum [88], Scutellaria baicalensis [37], Eucommia ulmoides [89], Moringa peregrina [90], etc.
The lowest values of gene diversity were detected in the populations from Mt. Biokovo (pops. 16, 17 and 18) (H E = 0.092-0.096), which were also identified as a separate gene pool by the Bayesian cluster analysis. The Mt. Biokovo populations (pops. 16, 17 and 18) were sampled in geographically and ecologically marginal areas of the species distribution, and it is generally expected that in such peripheral populations, the amounts of diversity may be reduced [91]. This reduction can be the result of, small population size, increased isolation and consequently reduced gene flow [92][93]. At such high elevation (greater than 1200 m.a.s.l.) with atypical climatic conditions for the species, it is highly unlikely that the spontaneous formation and development of populations of such a pronouncedly thermophytic species could be expected in the present day. However, during some past periods that were warmer than present (Holocene Climatic Optimum) [94], such altitudinal shifts in species distribution are quite expected and have already been documented for the western Balkans [10]. Therefore, these populations can be treated as relict populations, which are ecologically isolated from other populations. Evidently, a sufficient period of time has passed for this group of populations to diverge and form a separate, well-defined gene pool.
Evidence that both geography and environmental heterogeneity contributed significantly to the obtained patterns of genetic divergence was found. The IBD analysis showed positive and significant correlations, which could in principal be the result of successive postglacial colonization from the refugia [28]. However, a large part of unexplained genetic variations among the populations still remains after accounting for the geographical distances. Furthermore, a significant correlation between environmental and genetic variation was found, and 12.3% of the genetic variation was explained by environmental heterogeneity. The populations from Mt. Biokovo (pops. 16, 17 and 18), located at significantly higher elevation (from 1295 to 1335 m.a.s.l.), inhabiting environmentally distinct habitat in comparison to the remaining populations, are typical cases of environmental isolation. It is well acknowledged that microclimate strongly influences plant phenology [95][96] whereas differences in temperature, photoperiod and moisture are major factors in shaping patterns of flowering and fruiting [97]. Due to cooler temperatures and higher precipitation the populations from Mt. Biokovo tend to flower later (cca. 20 days) than do the lowelevation populations (authors observations); therefore, there is no phenological overlap between them, resulting in genetic isolation and differentiation of Mt. Biokovo populations (pops. 16,17,18). In addition, great influence of sea surface on coastal area ecological conditions and consequently on vegetation cover [25] could have also served as a stimulus for such strong differentiation between high altitude and lowland populations. Bathymetric configuration of Adriatic Sea floor is characterized by very deep valley with steep slopes in southern region (South Adriatic Basin) and relatively shallow and plain central and northern region [98][99]. For that reason northern Adriatic part has experienced more pronounced substantial shifts in shoreline during the sea-level fluctuations, while the southern regions remained almost unchanged [100]. Therefore, in the past when the sea level was significantly lower than today [99,101], populations from middle and especially southern parts of Dalmatia were located in proximity of the sea as they are today. However, this is not the case for high altitude areas as they were never located near sea surface, which might have resulted in divergent evolution particularly observed in Mt. Biokovo populations (pops. 16, 17 and 18). The strong influence of environmental factors on genetic divergence of populations has already been documented for several plant [77,102] and animal species [31,33]. For example, Temunović et al. [77] determined genetic structuring of Croatian Continental and Mediterranean ash trees populations (Fraxinus angustifolia Vahl) due to contrasting ecological conditions of the sampling sites. Similarly, strong fine-scale spatial genetic structure associated with habitat factors was detected in continuous population of Biscutella laevigata [102]. For the Franciscana dolphins genetically isolated populations in environmentally distinct areas of Argentina were observed [33] and for the grey wolf (Canis lupus) genetic differentiation among local populations was correlated to contrasting climatic conditions [31].
In our study the detected IBED pattern explained more than twice the genetic divergence than did the IBD pattern (12.3 vs. 5.6%, respectively). Therefore it is clear that the revealed genetic differentiation in Dalmatian pyrethrum is more strongly influenced by the bioclimatic conditions of the sampling sites rather than by geographic distances per se.

Conclusion
The aim of this study was to analyze the genetic diversity and structure of natural Dalmatian pyrethrum (Tanacetum cinerariifolium Trevir. /Sch./ Bip.) populations from Croatia. Overall, our results show that, as a thermophytic plant, Dalmatian pyrethrum could have survived the Quaternary climatic oscillation along the eastern Adriatic coast by tracking suitable habitats, i.e., microrefugia. Northern Adriatic populations, which were identified as a separate cluster by the STRUCTURE analysis, are characterized by high levels of gene diversity and DW and should be treated as conservation priority sites. Human influence through overexploitation and cultivation has resulted in the low genetic variability of southern Adriatic populations. The results also provide evidence that both IBD and especially IBED contributed to the genetic differentiation among some populations. Figure S1 Genetic structure of T. cinerariifolium populations. (A) using software Structure and assuming K = 2, (B) using Structure and assuming K = 3; (C) using BAPS without spatially informative prior, and (D) using BAPS with spatially informative prior. Each individual plant is represented by a single vertical line divided into colors. Each color represents one cluster, and the length of the colored segment shows the individual's estimated proportion of membership in that cluster. White lines separate populations that are labelled below the figure (1-20). (TIF)