Phytochemicals and Their Correlation with Molecular Data in Micromeria and Clinopodium (Lamiaceae) Taxa

A study of the phytochemical and molecular characteristics of ten Micromeria and six Clinopodium taxa (family Lamiaceae) distributed in the Balkan Peninsula was carried out. The phytochemicals detected in essential oils by gas chromatography, mass spectrometry, and molecular data amplified fragment length polymorphism were used to study the taxonomic relationships among the taxa and the correlations between phytochemical and molecular data. STRUCTURE analysis revealed three genetic groups, while Bayesian Analysis of Population Structure grouped the studied taxa into 11 clusters nested in the groups obtained by STRUCTURE. Principal components analysis performed with the 21 most represented compounds in the essential oils yielded results that were partly consistent with those obtained by STRUCTURE and neighbour-joining analyses. However, their geographic distributions did not support the genetic grouping of the studied taxa and populations. The Mantel test showed a significant correlation between the phytochemical and genetic data (r = 0.421, p < 0.001). Genetic distance explained 17.8% of the phytochemical distance between populations. The current taxonomic position of several of the studied taxa is yet to be satisfactorily resolved, and further studies are needed. Such future research should include nuclear and plastid DNA sequences from a larger sample of populations and individuals.


Introduction
The genus Micromeria Benth. (Lamiaceae) includes 54 [1], 70, or only 20 [2] annual and perennial herbs, sub-shrubs, and shrubs, depending on the point of view. According to Bräuchler et al. [1,3], the distribution of Micromeria species extends from the Mediterranean to South Africa and Madagascar and from China to the Macaronesian Archipelago. Chater and Guinea [4] described 21 Micromeria species for Europe, with more than half of these species occurring in the Balkan Peninsula. The genus Micromeria is part of a complex and 0.131 (population Mj6 of M. juliana). With an average value of 0.114, M. cristata and M. juliana had the highest H E levels, while M. longipedunculata and M. croatica were on the other side of the spectrum, with an average H E value of 0.084. Of the total variability, 48.68% refers to variability between and 51.32% within populations, indicating significant differences between the studied populations. AMOVA analysis showed that intrapopulation variability was considerably higher than among populations (Table 2). Variability among populations ranged from 5.04% (M. graeca ssp. graeca) to 34.45% (M. kerneri). For most species, interpopulation variability was approximately 10%. The exceptions were the populations of M. cristata (24.82%) and M. kerneri (34.45%), which had a higher and significant interpopulation distance (ΦST).
The results of neighbour-joining (NJ) analysis are shown in Figure 1. Three genetic groups determined by STRUCTURE were marked using the same colours (blue, green, red) in the NJ tree to allow comparison of the two analyses. The NJ analysis gives significant bootstrap support to confirm that Micromeria and Clinopodium are well-differentiated groups of closely related taxa. In addition, most individuals in both groups were wellsupported, except for the M. juliana-M. kerneri-M. microphylla cluster, characterised by low differentiation among the adopted taxa. Individuals of M. juliana were grouped in the same cluster with M. kerneri from the Croatian population Gradina (Mk3). The population of M. microphylla (Mm) was divided into two clusters and associated with four populations of M. kerneri. The AMOVA results showed no statistically significant differences between M. microphylla and M. kerneri. Two populations of M. croatica, known from local Balkan literature as M. pseudocroatica (McrP), were separated from the remaining ten populations of M. croatica (Mcr) (Figure 1, Table 3). Additionally, the population of M. graeca ssp. fruticulosa (MgF) was recognised as a distinct taxon, separated from the two populations of M. graeca ssp. graeca (Mg). In the NJ tree, the population of M. cristata ssp. kosaninii (McK) and two populations of M. cristata ssp. cristata (Mc3, Mc5) were also separated.
The position of the studied taxa within the genus Clinopodium is complex. Clinopodium dalmaticum showed genetic differentiation into two subgroups. The first subgroup was formed by two Montenegrin populations (Cd), while the second was formed by three Bulgarian populations (CdB). Moreover, the Bulgarian populations were closer to C. frivaldszkyanum (Cf) than to the Montenegrin populations of the same species ( Figure 1). Three other Clinopodium species studied (C. serpyllifolium, C. pulegium, and C. thymifolium) were separated from the other taxa with high bootstrap support.
In the STRUCTURE analysis, the highest K value was observed for K = 3 (K = 1394.43), indicating the presence of three genetic clusters ( Figure S1). STRUCTURE analysis ( Figure 2A) revealed three genetic groups, shown in blue, green, and red on the NJ tree ( Figure 1). The first group included populations of Micromeria cristata, M. croatica, M. graeca, and M. longipedunculata; the second included the populations of M. juliana and M. kerneri; the third included all Clinopodium populations studied. While low levels of admixture characterised most of the studied populations, this was not the case with the M. microphylla population that was positioned between the two Micromeria clusters with high admixture levels. The genetic clusters were not geographically defined, since representatives of each cluster are found throughout the studied area.   pulegium; Cs, C. serpyllifolium; Ct, C. thymifolium; three genetic groups determined by STRUCTURE were marked using the same colours (blue, green, red) in the NJ tree to allow comparison of the two analyses.
On the other hand, BAPS analysis ( Figure 2B) revealed a congruent assignment of the studied Micromeria and Clinopodium taxa to 11 clusters nested within the groups identified by the STRUCTURE analysis. The best partitions received log-likelihoods of −182,699.06 at P = 1 (without using geographic coordinates as informative priors) and −183,284.74 at P = 1 (with spatial clustering). In general, both methods produced nearly identical results. The first two groups of the BAPS analysis were formed by M. cristata ssp. cristata
PCA analysis (Figure 3) was performed on the 21 compounds isolated from the EO, in amounts exceeding 10% per sample (population). PC1 and PC2 for the EO compounds explained 29.26% of the variance. Pearson's correlation coefficients between 21 main compounds and scores of the first two PC are shown in Table 4. The phytochemical groups obtained by the PCA were partly similar to the three genetic groups determined by NJ and STRUCTURE analyses. Among the 21 compounds, PCA detected nine components that contributed most to the differences between groups (Figure 3). Clinopodium species were mainly located in the negative region of PC1 and PC2, while Micromeria longipedunculata (Ml) was positioned among the Clinopodium species.
The main compounds in this group were menthone, pulegone, and piperitenone oxide, and the highest values were found in Clinopodium frivaldszkyanum (Cf1, menthone) and C. thymifolium (Ct1, pulegone; Ct7, piperitenone oxide) (Tables S7 and S8). Only C. dalmaticum from Bulgaria (formerly Micromeria bulgarica, CdB) had an unusual position among the Clinopodium species, positioned near M. kerneri and M. juliana in the negative region of PC1 and the positive region of PC2. The specific compounds for this phytochemical group were verbenone, caryophyllene oxide, and docosane, which were highest in three populations of M. juliana (Mj7, Mj8, Mj5) (Table S4).        (Tables S2 and S3).

Mantel Test
The correlations between AFLP and the phytochemical matrices of dissimilarity were calculated using the Mantel test. A significant correlation was observed between the phytochemical and molecular data (r = 0.421, p Mantel < 0.001). According to the same test, 17.8% (R = 0.178) of the phytochemical distance between populations could be explained by genetic distance ( Figure S2).

Discussion
Several phylogenetic conclusions can be drawn from the genetic diversity and STRUC-TURE results. The NJ analysis separated Micromeria from Clinopodium taxa and reinforced the recent transfer of species from the section Pseudomelissa (genus Micromeria) to the genus Clinopodium by Bräuchler et al. [15]. Although the distinction was confirmed between the Micromeria and Clinopodium groups, it is questionable whether it is substantial enough to label these groups as separate genera. STRUCTURE analysis indicated the presence of three genetic clusters: two within Micromeria and a third of Clinopodium species. If Clinopodium is treated as a separate genus, then the other two Micromeria groups should also be considered as such. Since the previous Clinopodium-Micromeria segregation was based on a small sample size and analysis of a single cpDNA region [3], the results presented here are even more relevant.
The NJ analysis showed genetic differentiation of Clinopodium dalmaticum in the Montenegrin and Bulgarian populations. Regarding the EO composition, PCA also separated Montenegrin from Bulgarian populations of C. dalmaticum (Figure 3). These results suggest variability within C. dalmaticum, with the note that the Bulgarian populations were previously considered to be Micromeria bulgarica [19,20]. Variability within C. dalmaticum was previously described by Vandas [45], who identified M. bulgarica in the area of Usunža and Krivska River (North Macedonia). Chater and Guinea [4] and Ančev [19] identified two subspecies of M. dalmatica (now C. dalmaticum): M. dalmatica ssp. dalmatica and M. dalmatica ssp. bulgarica (Velen.) Guinea. On the other hand, Bräuchler et al. [1,15] concluded that M. bulgarica is a synonym of C. dalmaticum. Since C. dalmaticum showed both molecular and phytochemical separation among its populations, research on this species should be continued. The Bulgarian population of C. dalmaticum was closer to C. frivaldszkyanum than to the Montenegrin populations of C. dalmaticum (Figure 1). This suggests that the taxonomic relationships within these two or three taxa of the genus Clinopodium require further clarification. Given the considerable geographic distance between the Montenegrin and Bulgarian populations of C. dalmaticum, the existence of two geographically distinct groups is not unusual. Such a refugia-within-refugia model developed by Gómes and Lunt [46] for the Iberian Peninsula has also been applied to some species from the Balkan Peninsula [47][48][49].
In studies of individual species or groups of closely related species in the Balkan Peninsula, the question of the presence of (micro)refugial areas that protected local populations during unfavourable climatic conditions of glaciation cycles cannot be avoided. In the AFLP analyses, DW markers (Table 1)  Similar to the STRUCTURE results, no spatial structuring was observed of populations characterised by high frequencies of DW markers, suggesting that no single refugial area can be identified within the Balkan Peninsula. Instead, there appear to have been numerous microrefugia scattered over large areas. The northern Adriatic coastal area does not harbour any of this microrefugia, as low levels of these markers characterised these populations.
The NJ analysis separated Micromeria species into several genetic groups and raised questions about the systematic position of certain taxa. Populations of M. cristata were separated into two statistically significant different subgroups ( Figure 1, Table 3 (Figure 3). Additionally, the difference between M. croatica and M. pseudocroatica was greater than the difference among the ten populations of M. croatica. Although Bräuchler et al. [1] concluded that M. pseudocroatica is only a synonym of M. croatica, future genetic research should verify whether the differences presented herein might be due to geographical isolation. Namely, populations of the disputable M. pseudocroatica are located on the Pelješac Peninsula and are partly isolated from continental populations of M. croatica.
The population of Micromeria microphylla (Mm) was divided into two clusters associated with four populations of M. kerneri, although the differences between the species were not statistically significant (Table 3). On the other hand, the PCA analysis of EO compounds showed that M. microphylla (Mm) is quite different from populations of M. kerneri (Figure 3). Different habitat conditions can also explain the differences in EO composition between these species. The complexity in this group is further increased by one population of M. kerneri (Mk3), which is divided into two clusters closely related to M. juliana (Figure 1). The AMOVA (Table 3)  STRUCTURE analysis (Figure 2A) detected a general distribution in three genetic groups. Their geographic distributions did not support the genetic grouping of the studied taxa and populations. Surprisingly, there was virtually no spatial structuring of the recognised genetic clusters from the STRUCTURE analysis, as representatives from all groups are found across large regions, mixed in a seemingly random fashion. Such a result is hard to explain, but it confirms that the distribution ranges of the studied taxa do not follow the levels of their relatedness. It should be noted that a similar situation was also detected, where numerous populations of genetically well-supported taxa were scattered over large areas without any signs of spatial groupings. Such a result suggests the presence of strong gene flow barriers among closely related taxa, eliminating any possibility for interspecies hybridisation and consequent fusion of these taxa into spatially structured clusters. The exception is mentioned in the Micromeria juliana-M. kerneri-M. microphylla complex, where these barriers are weak at best. As such, there is currently no clear explanation for the obtained results. This is possibly a consequence of contrasting evolutionary histories and environmental conditions experienced by these taxa that have resulted in the development of strong reproductive isolation mechanisms.
Another set of results enabled a more straightforward conclusion. Not only were the majority of analysed taxa well-supported (except for the Micromeria juliana-M. kerneri-M. microphylla complex), but the recognition of a few additional taxa is now possible, thus opening the possibility for systematic repositioning within the studied groups. Within the C. dalmaticum group, two well-supported taxa were identified, one formerly known as M. EO content is known to depend on the developmental stage of the plant and the collection site [50]. To exclude the influence of the plant's developmental stage, the aboveground plant parts of all investigated taxa were collected for isolation of EO during flowering time. The composition of EO of Micromeria and Clinopodium taxa were investigated in all populations of the studied taxa. In general, the results presented in this study are consistent with patterns reported in the literature. In the composition of EO of M. cristata collected in Serbia, Bulgaria, Greece, and North Macedonia, the most abundant compound was borneol (14.11-26.28%) (Table S1). Kostadinova et al. [51] also identified borneol (6.1%) in a sample of M. cristata from Bulgaria, while its isomer isoborneol (11.3%) was the most abundant in the sample collected in Serbia [52]. The extent to which subspecies can differ in EO composition was shown by Çarikçi [53], who studied three subspecies of M. cristata, namely M. cristata ssp. cristata, M. cristata ssp. phyrigia P. H. Davis, and M. cristata ssp. orientalis P. H. Davis. In all three subspecies, the main constituents of EO were borneol and camphor [53]. Thus, the compounds α-muurolol and pulegone that were predominant in the taxon M. cristata ssp. kosaninii (Table S1) were not identified in the subspecies from Turkey. These differences are not surprising, considering the geographical distance and habitat conditions. Borneol was also one of the main compounds in the studied samples of Micromeria croatica (Tables S2 and S3) from Croatia, Bosnia and Herzegovina, and Montenegro, followed by the compounds E-caryophyllene and caryophyllene oxide. Caryophyllene oxide was the main compound in most of the studied populations of M. croatica, according to Slavkovska et al. [34], Kremer et al. [54], and Vuko et al. [55]. The EO of M. graeca ssp. graeca analysed here was rich in α-bisabolol (Table S3), while in the same taxon from Greece, the main component was epi-α-bisabolol [56]. In the composition of the ten EO samples of M. juliana, the most abundant volatile components were E-caryophyllene and caryophyllene oxide (Table S4). Similarly, these two compounds also dominated the EO composition of M. juliana from Anatolia, Turkey [53]. Caryophyllene oxide (12.81-23.46%) was the most abundant compound in the five M. kerneri oils studied, followed by α-pinene (12.3-16.13%) (Table S5). A previous study also showed that the EO composition of M. kerneri and M. juliana was characterised by a high concentration of oxygenated sesquiterpenes, with caryophyllene oxide as the most abundant compound [42]. In this study, the composition of EO of M. microphylla was reported for the first time, dominated by eudesem-7-(11)-en-4-ol (22.91%) ( Table S6). The peculiarity of the oil composition is not surprising, considering the isolation of this population in the central Adriatic (Table 1, Figure 4).
subspecies of M. cristata, namely M. cristata ssp. cristata, M. cristata ssp. phyrigia P. H. Davis, and M. cristata ssp. orientalis P. H. Davis. In all three subspecies, the main constituents of EO were borneol and camphor [53]. Thus, the compounds α-muurolol and pulegone that were predominant in the taxon M. cristata ssp. kosaninii (Table S1) were not identified in the subspecies from Turkey. These differences are not surprising, considering the geographical distance and habitat conditions.
Borneol was also one of the main compounds in the studied samples of Micromeria croatica (Tables S2 and S3) from Croatia, Bosnia and Herzegovina, and Montenegro, followed by the compounds E-caryophyllene and caryophyllene oxide. Caryophyllene oxide was the main compound in most of the studied populations of M. croatica, according to Slavkovska et al. [34], Kremer et al. [54], and Vuko et al. [55]. The EO of M. graeca ssp. graeca analysed here was rich in α-bisabolol (Table S3), while in the same taxon from Greece, the main component was epi-α-bisabolol [56]. In the composition of the ten EO samples of M. juliana, the most abundant volatile components were E-caryophyllene and caryophyllene oxide (Table S4). Similarly, these two compounds also dominated the EO composition of M. juliana from Anatolia, Turkey [53]. Caryophyllene oxide (12.81-23.46%) was the most abundant compound in the five M. kerneri oils studied, followed by α-pinene (12.3-16.13%) (Table S5). A previous study also showed that the EO composition of M. kerneri and M. juliana was characterised by a high concentration of oxygenated sesquiterpenes, with caryophyllene oxide as the most abundant compound [42]. In this study, the composition of EO of M. microphylla was reported for the first time, dominated by eudesem-7-(11)-en-4-ol (22.91%) ( Table S6). The peculiarity of the oil composition is not surprising, considering the isolation of this population in the central Adriatic (Table 1, Figure 4).  Clinopodium dalmaticum is endemic to the Balkan Peninsula, and is widespread in Bulgaria, Montenegro, and Greece, including Crete [57]. In this study, the volatile compounds of samples collected in Montenegro and Bulgaria were identified. In the composition of the isolates from Montenegro, the predominant compound was piperitone, while the Bulgarian samples were rich in E-caryophyllene, α-pinene, and thymol (Table S7). The most abundant compounds in the EO of C. frivaldszkyanum, C. pulegium, C. serpyllifolium, and C. thymifolium were pulegone and piperitenone oxide, making oils of these species extremely rich in oxygenated monoterpenes (53-85.31%) (Tables S7 and S8). According to Zheljazkov [58], pulegone was one of the main constituents in the EO of C. frivaldszkyanum from the Bulgarian populations of Shipka and Uzana.
The Mantel test showed a significant correlation between the phytochemical and AFLP data (r = 0.421, p < 0.001). The literature on this topic is diverse. According to Slavkovska et al. [34], the composition and quantity of EO of Micromeria species distinguished section Pseudomelissa from Eumicromeria. The EO of species from the section Pseudomelissa was dominated by oxygenated monoterpenes of the menthane type, while various terpene compounds were dominant in species from the section Eumicromeria [34]. Multivariate analysis (PCA and UPGMA) of compositions determined in the EO of M. kerneri and M. juliana separated the populations of these two species [42]. Feulner et al. [36] determined strong and significant correlations between AFLP data and floral scent volatiles at the population level (r = 0.791, p = 0.004) and individual level (r = 0.823, p < 0.001) in Sorbus taxa (family Rosaceae). Xavier et al. [59] found a significant correlation between volatile chemical classes and genetic traits of Aniba Aubl. species in the Amazon region in Pará State (Brazil). Additionally, AFLP profiles of 11 Hypericum species and cultivars were correlated with their levels of phytochemical markers (chlorogenic acid, hyperforin, hypericin, pseudohypericin, and rutin) determined in their methanolic extracts enabling true-to-type identification and marker-assisted breeding programmes [60]. Investigations of 20 populations of four Thymus L. species native to Hungary found only partial similarities between dendrograms generated by hierarchical cluster analysis based on DNA patterns and EO samples [61]. On the other hand, in a study of Ophrys L. taxa (Orchidaceae), Stökl et al. [62] did not find any correlation between scent data and DNA-molecular data. Trindade et al. [63,64] concluded that there was no correlation between the chemical analysis of EO of Thymus caespititius Brot. from the Azores and molecular data sets. Volatile and molecular analysis of Juniperus brevifolia (Seub.) Antoine from the same archipelago also showed no correlation between chemical and molecular data sets [65]. Finally, Emami-Tabatabaei et al. [66] studied the possible correlation between AFLP data and the EO profile obtained by GC-MS of Lutea elbursensis Mozaff from northern Iran, concluding that the chemical composition of EO cannot be used as a reliable taxonomic tool.

Plant Material
Randomly selected samples of wild-growing plants of Micromeria and closely related Clinopodium species were collected during the blooming period from June to August 2018. Voucher specimens of herbal material were deposited in the Fran Kušan Herbarium, Faculty of Pharmacy and Biochemistry, University of Zagreb, Croatia (Table 1, Figure 4). For molecular analysis, several young leaves from 3 to 11 plants per population were collected on a dry day. Immediately after collection, leaves were dried in plastic bags containing silica gel and stored for further use in DNA analysis. Additionally, above-ground shoots with leaves and flowers were harvested and mixed to obtain a randomly selected sample. The collected plant parts were air-dried and protected from direct sunlight for 15 days at 22 • C and 60% relative humidity. From each locality, 50 g of air-dried plant material was hydro-distilled for 3 h in a Clevenger-type apparatus. The EO obtained was dried over anhydrous sodium sulphate.

AFLP Data Analysis Within-Population Diversity
To construct a binary matrix, the obtained AFLP fragments were scored as present (1) or absent (0). Diversity within populations was assessed by calculating the proportion of polymorphic markers (%P), the number of private markers (N pr ), and the frequency downweighted marker values (DW) [69] using the AFLPdat package in R [70]. The Shannon information index of each population was calculated as I = −Σ (p i log 2 p i ), where p i is the phenotypic frequency [71,72]. In addition, genetic diversity (H E ) was calculated using a Bayesian approach [73], assuming the Hardy-Weinberg equilibrium due to outcrossing (F IS = 0) as implemented in AFLP-Surv v. 1.085 (Vekemans, X., Laboratoire de Génétique et Ecologie Végétale, Université Libre de Bruxelles, Bruxelles, Belgium) [74]. The overall mismatch error rate for all AFLP primer combinations was 2.5%.
Analysis of molecular variance (AMOVA) [78] was used to partition the total genetic variance among and within populations of each taxon and between closely related taxa, among populations within taxa and within populations. The variance components were tested with 10,000 permutations in Arlequin ver. 3.5.2.2 (Excoffier, L., Lischer, H., Institute of Ecology and Evolution, University of Berne, Bern, Switzerland) [79].
Population structure was assessed using two Bayesian clustering approaches implemented in STRUCTURE v2.3.4 (Pritchard Lab., Stanford University, Stanford, CA, USA) [80] and BAPS v6.0 (Corander, J., Cheng, L., Marttinen, P.; Sirén, J.; Tang, J., Department of Mathematics and Statistics, University of Helsinki, Helsinki, Finland) [81]. In STRUC-TURE, 30 runs were performed for each K by setting the number of clusters (K) from 1 to 21. Each run consisted of a burn-in period of 200,000 steps followed by 1,000,000 Monte Carlo Markov Chain (MCMC) replicates assuming an admixture model and correlated allele frequencies. The calculations were performed on the Isabella computer cluster at the University of Zagreb, University Computing Centre (SRCE). The most probable number of K was selected by calculating ∆K [82] in StructureSelector [83]. StructureSelector was also used to cluster and average the results of independent runs using the approach described by Kopelman et al. [84]. BAPS was applied for population mixture analysis without the geographic origin of samples as an informative prior ('Clustering of Individuals') and with this prior ('Spatial Clustering of Individuals') [85]. The maximum number of clusters (K) was set to 20, and each run was replicated 10 times. Population admixture analysis [86] was performed with the default settings.

Gas Chromatography and Mass Spectrometry (GC-MS) Analyses
The EO of each Micromeria and Clinopodium sample obtained by hydro-distillation were collected for each sample in a pentane/diethyl ether mixture and analysed by gas chromatography and mass spectrometry (GC-MS). GC was performed using a gas chromatograph (model 3900; Varian Inc., Lake Forest, CA, USA) and a mass spectrometer (model 2100T; Varian Inc.). The MS conditions were ion source temperature 200 • C, ionisation voltage 70 eV; mass scan range: 40-350 mass units. The carrier gas was helium. Two columns were used: nonpolar VF-5 ms and polar capillary columns CP Wax 52. The conditions for the VF-5 ms column were temperature 60 • C isothermal for 3 min, then increased to 246 • C at a rate of 3 • C min −1 , and held isothermal for 25 min. The CP Wax 52 column conditions were: temperature 70 • C isothermal for 5 min, then increased to 240 • C at a rate of 3 • C-min −1 , and held isothermal for 25 min. The injection volume was 2 µL and the split ratio was 1:20. The triplicate analyses of individual peaks were identified by comparing the retention indices of the n-alkanes with literature data and authentic standards [87,88].

Principal Component Analysis
Principal component analysis (PCA) was based on 21 significant constituents of the EO. PCA was performed using the PRINCOMP procedure in SAS v9.3 (SAS Institute Inc., Cary, NC, USA) [89], and the biplot showing the populations and oil constituents (as vectors) was constructed using the first two principal components.

Mantel Test
The Mantel test [90] was used to test the correlation between genetic and biochemical data matrices. Pairwise genetic distances between populations were calculated using Nei's standard genetic distance (DNei) in AFLP-Surv v1.085 [69]. Biochemical differences were calculated as Euclidean distances between populations for the first two principal components of the PCA of EO constituents. The significance level was assessed after 10,000 permutations in NTSYS-pc v2.21L [91].

Conclusions
STRUCTURE analysis based on AFLP genetic data grouped the studied ten Micromeria and six closely related Clinopodium taxa distributed in the Balkan Peninsula into three genetic groups. The first cluster included all Clinopodium taxa, while Micromeria species were divided into two clusters. In general, their geographic distributions did not support the genetic grouping of the studied taxa and populations. Numerous populations of genetically well-supported taxa were also found scattered over large areas with no evidence of spatial groupings. Such a result suggests that substantial gene flow barriers exist between closely related taxa, precluding any possibility for inter-species hybridisation and consequent fusion of these taxa into spatially structured clusters. An exception is the M. Juliana-M. kerneri-M. microphylla complex, where these barriers are weakest. Generally, groups of taxa were much less supported than individual taxa, indicating their concurrent dispersal and approximately the same time of origin.
The results also showed that the current taxonomic position of certain species requires stronger resolution. Within the C. dalmaticum group, two well-supported taxa were identified, one formerly known as M. dalmatica and the second as M. bulgarica. The species M. graeca ssp. fruticulosa (formerly M. fruticulosa) emerged as a well-supported taxon and not a representative of M. graeca. Two populations of M. croatica, formerly recognised as M. pseudocroatica, were also clearly differentiated from any other taxon. Within the M. cristata group, the taxon M. cristata ssp. kosaninii (formerly M. kosaninii) lacks the needed support for its recognition as either a species or subspecies. Although further studies are needed on some species within the genera Clinopodium and Micromeria, the AFLP data obtained in this research provide a good starting point for future studies. Such a study should include nuclear and plastid DNA sequences on a larger sample of populations and individuals. Finally, the Mantel test showed a significant correlation between phytochemical and AFLP data.