Study species

Salvia officinalis is a thermophilous evergreen perennial well adapted to the Mediterranean climate. Due to its shade-intolerance and low competitive ability, it thrives on sunny slopes with poor, rocky soil that are unsuitable for other more competitive plant taxa (Devetak 1963). Its natural range, as described by Greuter & al. (1986), Karousou & al. (2000) and Kintzios (2000), spans continuously along the E Adriatic coast, from Italy to Albania, rising from see level to more than 1000 m a.s.l. (Fig. 1 & 2A). Four disjunct inland populations of S. officinalis are found in gorges in the central Balkans. These sites belong to different biogeographic regions and vegetational zones as compared to those along the E Adriatic coast (Table 1) and are characterized by moderate continental climate suitable for S. officinalis and other sub-Mediterranean elements.

Salvia officinalis reproduces sexually and has poor ability for vegetative propagation (Devetak 1963). It is an outcrossing species with zygomorphic flowers well adapted to pollination by insects (Hymenoptera sp.; Fig. 2C). Although its pollinating mechanism has long been of interest to naturalists, the mode of seed dispersal is still unclear. Fructification occurs in plants more than three years old and seeds are c. 3 mm in diameter and rather heavy (the mass of 1000 seeds is c. 10 g, Devetak 1963). Upon ripening in August and September, seeds from the mother plant usually disperse via barochory to distances of several meters only, suggesting a rather limited dispersal ability.

Sampling

In spring 2010, Salvia officinalis samples were collected from a total of eight populations. Details of sampling localities, their biogeographic regions and vegetation zones are given in Table 1; their locations are shown in Fig. 1. Four sampled populations (LUS, MOR, ORJ and RUM) are located within the SE Adriatic part of the species range. The other four (JBL, MIL, PIV and SIC) are disjunct populations found in gorges in the central Balkans. JBL and PIV are geographically close to the sampled Adriatic populations, while MIL and SIC represent the most distant populations found far inland in the continental Balkans.

Young leaves were collected from nine to 11 individuals from each population; the total number of samples was 83. Sampled individuals, which were randomly selected from all populations, were evenly distributed within populations and distant at least 20 m from each other. In the field, leaves were labelled and packed in sterile coffee-filter bags which were deposited in larger bags with silica gel for drying. Fresh leaves were dried in silica gel for seven days and kept in a dark, dry place at 10 °C prior to DNA isolation.

DNA isolation, PCR amplification and sequencing

Total genomic DNA was isolated following Protocol 2 of Aleksić & al. (2012). Two chloroplast intergenic spacers, 3′rps16-5′trnK and rpl32—trnL, located in large and small single copy regions, respectively, were amplified using primers described in Shaw & al. (2007): a forward primer rpS16x2F2: 5′ - AAAGTGGGTTTTTATGATCC - 3′ and a reverse primer trnK^(UUU)xl: 5′ - TTAAAAGCCGAGTACTCTACC - 3′ for amplification of the former region, and a forward primer trnL^(UAG): 5′ - CTGCTTCCTAAGAGCAGCGT - 3′ and a reverse primer rpL32-F: 5′ - CAGTTCCAAAAAAACGTACTT C - 3′ for amplification of the latter region. PCR amplifications of both regions were performed in 25-µl volumes containing: 25 ng template DNA; 1 × Taq Buffer with (NH₄)₂SO₄ (Fermentas, Vilnius, Lithuania); 2.5 mM MgCL₂; 0.2 mM dNTPs; 0.1 µM of each forward and reverse primer; 0.80% BSA (Bovine Serum Albumin, Promega, St. Louis, U.S.A.); and 0.025 U of Taq DNA polymerase (Fermentas, Vilnius, Lithuania). For both loci, we used identical PCR amplification profiles: an initial denaturation at 94 °C for 10 min; followed by 35 cycles of denaturation at 94 °C for 45 s; annealing at 53 °C for 1 min; extension at 72 °C for 1 min; and a final extension of 10 min at 72 °C. The presence of a specific PCR product was confirmed by gel electrophoresis on 2 % agarose gels. PCR products were sequenced with the forward primer by Macrogen Europe, Amsterdam, the Netherlands ( http://dna.macrogen.com/eng/) via Sanger sequencing using a 96-capillary 3730xl DNA Analyzer automated sequencer (Applied Biosystems Inc., U.S.A.). Electropherograms were edited manually using Chromas Lite ( http://www.technelysium.com.au/chromas_lite.html). Three and six different sequences detected in overall sample at rpl32—trnL and 3′rps16-5′trnK, respectively, were deposited in GenBank (accessions JQ771324 to JQ771332).

Genetic diversity estimates

Sequences generated from each cpDNA intergenic spacer were aligned separately by Muscle (Edgar 2004) in MEGA5 (Tamura & al. 2011), the alignments were refined manually, and used for assembling a concatenated matrix comprising both cpDNA regions.

Levels of inter- and intra-population genetic diversity were quantified by calculating the number of haplotypes (h), Nei's nucleotide diversity (π) and haplotype diversity (Hd) indices using DnaSP v.5.10 (Librado & Rozas 2009). We also calculated within-population diversity (H_s) and total gene diversity (H_T) using PERMUTE (Pons & Petit 1996).

Genetic differentiation, spatial genetic structuring and phylogeographic analyses

Pairwise population differentiation statistics (F_ST) were calculated in ARLEQUIN 3.5 (Excoffier & al. 2005), and statistical significance was evaluated using 1000 permutations. Those values were used to calculate the effective number of female migrants per generation as an indirect estimate of average levels of seed flow relative to genetic drift under the infinite island model. For haploid data, Nm = ½(1/F_ST — 1), where N is the female effective population size and m is the female migration rate (Slatkin 1993). Nm was calculated in the overall sample and between all pairs of populations.

Fig. 1.

Distribution of Salvia officinalis within the Balkans, studied populations and distribution and abundance of cpDNA haplotypes. The map was constructed using the Flora Croatica Database, BEOU herbarium database and authors' fieldwork observations. Data on distribution of S. officinalis in Albania are incomplete. Studied populations are marked with black circles.

We tested for a correlation between linear pairwise population genetic distances, as measured by F_ST, and the logarithm of geographic distances between pairs of populations using a Mantel test (Mantel 1967) in IBDWS (Isolation by Distance Web Service) v. 3.21 (Jensen & al. 2005) with 30 000 permutations to determine the significance of genetic differentiation by distance (isolation-by-distance, IBD pattern). Undefined values that occurred when identical, monomorphic populations were involved in the pairwise calculations were converted to 0.0001 for correlations.

SAMOVA analysis (Dupanloup & al. 2002) was used to define groups of populations that were geographically homogeneous and maximally differentiated from each other. SAMOVA does not require an a priori definition of populations, but instead searches for emergent group structures based solely on the genetic data. The program repeatedly seeks the composition of a user-defined K of groups of geographically adjacent populations that maximizes F_CT and minimizes F_SC. We ran 100 simulated annealing processes for K ranging from K = 2 to K = 6, and recorded the configuration with the largest F_CT values and lowest F_SC as the best grouping of populations. In addition, SAMOVA was used to identify barriers between populations.

A genetic distance matrix of pairwise population F_ST values was used to perform a hierarchical analysis of molecular variance (AMOVA, Excoffier & al. 1992) using ARLEQUIN 3.5. AMOVA was employed to estimate and partition the total molecular variance into three components: among lineages (F_CT), among population within each lineage (F_SC), and within populations (F_ST), whose statistical significance was tested by 10 000 permutations. To perform AMOVA, we assumed different groupings of populations as indicated in previous analyses. AMOVA may assist in finding the optimal grouping of populations recognized as the one with the highest F_CT and the lowest F_SC values, similarly as in SAMOVA, and can reveal homogeneity of detected groupings.

We used CONTRIB (Petit & al. 1998) to calculate genetic differentiation for unordered (G_ST) and ordered alleles (N_ST) and to test the phylogeographic signal across the species range by comparing these estimates. A phylogeographic structure is implied if N_ST is higher than G_ST because closely related haplotypes are more often found in the same geographic area than would be expected by chance. A statistical significance between G_ST and N_ST was obtained by 1000 permutations.

A statistical parsimony network revealing genealogical relations among cpDNA haplotypes was constructed using TCS 1.21 (Clement & al. 2000). Sequence gaps were treated as a fifth character state after recoding indels longer than 1 bp as single base-pair indels. TCS calculates the number of mutational steps among all pairs of haplotypes and then joins the most similar haplotypes into a network in which their combined probability is greater than 95 % (Templeton & al. 1992).

Phylogenetic analysis and divergence time estimates

For phylogenetic analysis we used a Bayesian approach in BEAST v.1.6.2 (Drummond & Rambaut 2007) because it accounts for the genealogical uncertainty due to the stochastic nature of the coalescence process and does not require potentially misleading outgroups for rooting. We constructed a majority rule consensus tree, assessed support for the clades with Bayesian posterior probabilities (PP) and estimated divergence time between detected lineages. Indels were excluded from these analyses. The best-fitting model of sequence evolution for our data set was determined using MEGA5. The smallest Bayesian Information Criterion (BIC) score was obtained for T92 and it was low for the HKY model as well, while the Akaike Information Criterion (AIC) supported the TN93 and HKY substitution models. Therefore, we used the HKY model in the BEAST analysis. To avoid over parametrization, we selected a demographic model of constant population size as a tree prior for modelling changes in population size through time. Tajima's relative rate test in MEGA5 was used to test the equality of evolutionary rates among three chosen haplotypes (A, B and F), which led to the acceptance of molecular clock hypothesis. In order to keep the number of estimated parameters to a minimum, we used a strict clock which is the most appropriate model for intraspecific data. The estimate of divergence time between detected lineages was based on published substitution rates of 1.01 × 10^-9 substitutions per site per year for synonymous sites of cpDNAs in seed plants (Graur & Li 2000). These values approximate the evolutionary rates of introns and noncoding spacers of organelle DNA (Chiang & al. 2009) and have been applied to another member of the mint family, Thymus herbabarona Loisel. (Molins & al. 2011). The input file for the BEAST analysis was constructed using the BEAUti interface of the BEAST package and the file with parameter settings was executed in BEAST. Following a burn-in of 1 000 000 steps, all parameters were sampled once every 100 steps from 5 000 000 MCMC steps. We used TRACER v.1.6.2 to test for stationarity of model parameters and to verify adequate sample sizes, which all exceeded 200, suggesting acceptable mixing and sufficient sampling. The posterior distributions of estimated divergence times for the nodes were summarized in TreeAnnotator 1.4.2, and trees were visualized using FigTree 1.0 (Rambaut 200).

Table 1.

Sampling localities of Salvia officinalis populations and genetic variation in two cpDNA regions.

Table 2.

Salvia officinalis haplotypes, haplotype frequencies and variant sites in a 1216-base-pair (bp) matrix comprising the cpDNA regions rpl32—trnL and 3′rps16—5′trnK.

Demographic analyses

The neutrality Fu's F_s test (Fu 1997) and Tajima's D test (Tajima 1989) implemented in ARLEQUIN 3.5 were performed for the Salvia officinalis lineages to detect departures from population equilibrium. The significance of both test statistics was tested with 1000 bootstrap replicates. The same program was used to assess the demographic history of lineages by a mismatch distribution analysis. Rapid expansion in the recent past is expected to generate a unimodal distribution of the pairwise differences between sequences, while multimodal distribution is more typical of a population in mutationdrift equilibrium (Rogers & Harpending 1992). The sum of squared deviations (SSD) between the observed and the expected mismatch distributions was computed and P values were calculated as the proportion of simulations producing larger SSD than the observed SSD. A significant P value rejects the fit of the data to the sudden expansion model.

CpDNA polymorphism and genetic diversity estimates

Both chloroplast intergenic spacers were successfully amplified in all Salvia officinalis individuals. The aligned sequence length of 3′rps16-5′trnK was 658 bp and harboured 14 informative sites (six transitions, four transversions and four length mutations comprising a four-bp indel, and three duplications of three- or four-bp motifs). Within the 558 bp of the aligned rpl32—trnL sequence, eight informative sites were found (six transitions, one duplication of a six-bp motif, and one poly-T stretch with a variable number of nucleotides, i.e. a microsatellite) (Table 2). The length of the concatenated matrix was 1216 bp harbouring 22 informative sites (base substitutions and length mutations).

In total, seven haplotypes were identified (A to G, Table 2), and their distribution among Salvia officinalis populations was non-random (Fig. 1; Tables 1 and 2). The most abundant haplotype B (frequency of 0.531, Table 2) was fixed in four inland populations (JBL, MIL, PIV and SIC, the latter population comprised an additional rare haplotype A closely related to B) and common in one Adriatic population (LUS). Haplotype C (frequency of 0.213) was exclusive in Adriatic populations. It was prevalent in populations LUS and MOR and rare in ORJ and RUM. Haplotype G (frequency of 0.193) was another exclusively Adriatic type found only in ORJ and RUM. Private rare haplotypes were found in one inland population (haplotype A, SIC) and in three Adriatic populations (haplotype D, MOR; haplotype E, ORJ; and haplotype F, RUM). The number of haplotypes per population ranged from one in the majority of inland populations to three in Adriatic populations ORJ and RUM (Tables 1 and 2).

Haplotype diversity (Hd) and nucleotide diversity (π) per population ranged from 0.000 (JBL, MIL and PIV) to 0.509 (LUS), and from 0.0000 (JBL, MIL and PIV) to 0.0036 (ORJ), respectively, with a total of 0.641 and 0.0035, respectively in the overall sample (Table 1). The average population genetic diversity, H_s, was 0.221 (SE = 0.074), and total gene diversity, H_T, of 0.695 (SE = 0.120).

Genetic differentiation, spatial genetic structuring and phylogeographic analyses

Pairwise population F_ST values ranged from 0.000 (between invariable inland populations fixed for an identical haplotype) to 0.934 (between Adriatic population MOR and two inland populations, MIL and PIV) (Table 3). F_ST values were non-significant between all pairs of inland populations and between two pairs of Adriatic populations, LUS-MOR and ORJ-RUM (Table 3). Accordingly, historical Nm estimated using F_ST values was > 1.00 among populations which were not significantly differentiated and < 1.00 between highly and significantly differentiated pairs of populations (Table 3).

Table 3.

Genetic differentiation and effective seed flow among Salvia officinalis populations.

The IBD pattern was not detected. Although the correlation coefficient between linear pairwise population F_ST values and the logarithm of the geographic distances between populations was positive, it was low and insignificant (r² = 0.003, P = 0.351).

SAMOVA analysis revealed the existence of at least two ancient lineages in Salvia officinalis. Assuming two lineages, at K = 2, F_CT value was high and significant (0.799, P = 0.03, Table 4) and F_SC value was relatively high and significant (0.336, P = 0, Table 4). All four inland populations (JBL, MIL, PIV and SIC) and two Adriatic populations (LUS and MOR) were distinguished from the remaining Adriatic populations (ORJ and RUM; the graphical representations of SAMOVA results are shown in Fig. 3). We delineated the former group of populations as inland-Adriatic lineage (IAL) and the latter group as purely Adriatic lineage (PAL). However, at K = 3, F_CT value slightly increased and remained significant (0.813, P = 0.003, Table 4) while F_SC value rapidly decreased and remained significant (0.012, P = 0, Table 4). Assuming three lineages in S. officinalis, populations comprising IAL would be further divided into two lineages, one comprising inland populations (JBL, MIL, PIV and SIC) and the other comprising Adriatic populations (LUS and MOR). Such a grouping of populations coincided with the biogeographic regions and vegetational zones in which all populations were found (see Table 1). The analysis, which takes into account geographic locations of populations in addition to genetic data, revealed that the percentage of molecular variation assigned to variation among populations within lineages was low but significant at K = 3.

Table 4.

Summary of AMOVA and SAMOVA statistics for alternative groupings of populations.

Fig. 2.

Salvia officinalis in nature. — A: part of the Adriatic natural population from Mt Orjen (ORJ), Montenegro; B: an individual plant from population ORJ; C: a zygomorphic flower of S. officinalis from population ORJ. — All photographs taken on 13 June 2012 by D. Stojanović.

Fig. 3.

Graphical representation of SAMOVA analyses and obtained groupings of populations at K = 2-5. Red lines indicate barriers between populations detected at each K. — A: K = 2; first group: ORJ & RUM; second group: JBL, LUS, MIL, MOR, OST & PIV. — B: K = 3: first group: ORJ & RUM; second group: LUS & MOR; third group: JBL, MIL, OST & PIV. — C: K = 4; first group: ORJ; second group: RUM; third group: LUS & MOR; fourth group: JBL, MIL, OST & PIV. — D: K = 5; first group: ORJ; second group: RUM; third group: LUS; fourth group: MOR; fifth group: JBL, MIL, OST & PIV.

Further increase of K in SAMOVA analyses resulted in separation of ORJ from RUM (K = 4), MOR from LUS (K = 5), and finally JBL from the other inland populations (K = 6). However, at K = 4-6, the increase of F_CT values was negligible (0.815, 0.818 and 0.800, respectively) as well as the decrease of F_SC values (-0.040, -0.095 and -0.090, respectively). In addition, at K = 2, SAMOVA suggested barriers between inland population PIV and geographically close Adriatic populations LUS and ORJ, and between RUM and its neighbouring populations (Adriatic populations LUS, MOR and ORJ, and an inland population JBL). At K = 3, barriers were found between PIV and all geographically close Adriatic populations, and between MOR and RUM but not between MOR, LUS and ORJ.

The partitioning of molecular variance into three components was performed assuming two and three lineages in Salvia officinalis (Table 4). In both analyses, high and similar amounts of molecular variance were attributed to variation among lineages (79.87 % and 81.30 %, respectively), and both F_CT values were high but at K = 2, F_CT value was marginally significant (0.799, P = 0.05) while it was significant at K = 3 (0.813, P = 0.03). However, under the assumption of three lineages, the percentage of variation assigned to variation among populations within each lineage was non-significant and lower (0.22%, F_SC = 0.012, P > 0.05) than under the assumption of two lineages (IAL and PAL) (6.77%, F_SC = 0.336, P = 0.00). Also, lower and non-significant F_ST values were obtained within each lineage under the assumption of three lineages than under the assumption of two lineages. Thus, AMOVA analysis supported the existence of three lineages in S. officinalis. Nonetheless, at K = 2, AMOVA revealed that inland-Adriatic lineage of S. officinalis was more heterogeneous than PAL (IAL: F_ST = 0.697, P = 0; PAL: F_ST = 0.004, P > 0.05).

The comparison between N_ST and G_ST showed a significant difference (N_ST = 0.764, SE = 0.046; G_ST = 0.682, SE = 0.066; P = 0.000), suggesting a phylogeographic structure.

Statistical parsimony analysis revealed that haplotypes B, C, D and E are at one- or two-mutational steps from the hypothetical haplotype 1 (hhl), which occupies the internal position at the haplotype tree and can be considered ancestral, while haplotype G was separated by 13 mutational steps from hhl (Fig. 4). Haplotypes A (differing from haplotype B in three mutations), B. C and D were exclusively found in populations comprising IAL with the exception of haplotype C, which was also found in two individuals from PAL. These haplotypes were all more closely related to hhl than was haplotype G. Haplotype G was exclusively found in SE Adriatic populations comprising PAL. PAL shared only haplotype C with IAL, and harboured two rare, private haplotypes, E (ORJ) and F (RUM). Based on statistical parsimony analysis, we were unable to infer whether two or three independent lineages were present in Salvia officinalis.

Phylogenetic analysis and divergence time estimates

The majority rule consensus tree from the strict clock Bayesian analysis with posterior probabilities of branches is presented at Fig. 5. In this tree, the root was placed between haplotypes F and G, which are both found in populations comprising PAL, and haplotypes A, B, C, D and E, which are all found in populations comprising IAL with the exception of haplotype E. Posterior probabilities (PP) of branches corresponding to these two lineages were very strong (1.00 and 0.99, respectively) unambiguously supporting the existence of two ancient lineages in Salvia officinalis, IAL and PAL. The grouping of haplotypes A and B found exclusively in inland populations (JBL, MIL, PIV and SIC) was strongly supported (PP of 1.00), whereas the separation of haplotypes C and D found in Adriatic populations LUS and MOR was poorly supported (PP of 0.58).

The BEAST dated the Most Recent Common Ancestor of the present sample to 3.267 Mya (lower 95 % Highest Posterior Density, HPD, = 1.584 Mya, upper 95 % HPD = 5.247 Mya), suggesting that the initial differentiation of the Salvia officinalis gene pool occurred sometime in the Pliocene. The diversification of IAL to inland and Adriatic sub-lineages was estimated to occur in the mid-Pleistocene (1.188 Mya (0.495 Mya, 2.079 Mya)), while the differentiation of PAL occurred 0.396 Mya (0.099 Mya, 0.990 Mya). Within IAL, differentiation of the Adriatic and inland sub-lineages was dated to 0.693 Mya (0.198 Mya, 1.287 Mya) and 0.495 Mya (0.099 Mya, 0.891 Mya), respectively.

Demographic analysis

Neither Tajima's D values, which were negative and non-significant for both lineages (-0.15338, P = 0.70383 in IAL, -0.55528, P = 0.40450 in PAL), nor Fu's F_s values, which were positive and non-significant for both lineages (0.89142, P = n.a. in IAL, 4.80511, P = 0.96100 in PAL), deviated from the expectation of neutrality (Table 5). The significant SSD values (0.06540, P = 0.01833 for IAL and 0.20714, P = 0.04559 for PAL) rejected a sudden expansion model of Salvia officinalis lineages (Table 5), as also revealed by mismatch distributions of both lineages, which were multimodal and inconsistent with the bell-shaped curve expected for an expanding population (data not shown).

Fig. 4.

Statistical parsimony analysis of two cpDNA regions in Salvia officinalis. Small black circles indicate hypothetical haplotypes and each solid line represents one mutational step; A–G = haplotypes found in studied populations; hhl = hypothetical haplotype that is ancestral to the overall sample.

Pliocene evolutionary history of Salvia officinalis

Walker & al. (2004) demonstrated recently that the genus Salvia, with over 900 species, has a polyphyletic origin, and that their Clade I, comprising the type species of the genus, Salvia officinalis, is monophyletic — a view supported also by further studies (Walker & Sytsma 2007; Will & Claßen-Bockhoff 2014). The centre of origin of S. officinalis, as inferred from the geographic distribution of Walker's Clade I, must be close to the current natural range of the species in the Balkans, a region characterized by a distinct and rather dynamic palaeogeographic history as compared to the rest of Europe (Pantić 1984; Utescher & al. 2007; Savić 2008) that has led to its current exceptional richness in diverse physiographic and climatic conditions (Stanković 1960).

The Balkan Peninsula, currently occupying an area of c. 520 000 km², arose from an archipelago within the Tethyan ocean and became connected to the S part of Palaeo-Europe since the late Oligocene (25 Mya, Pantić 1984). According to fossil records (Pantić 1984; Utescher & al. 2007), “Mediterranean” and “sub-tropical” taxa were rather abundant in the Balkans during the middle Miocene (15–11 Mya), while the presence of vast water bodies in the Para-Tethys realm in the late Miocene (11–6 Mya) enabled retention of a specific palaeoclimate and persistence of warmth-loving palaeoflora in this region as compared to W and C Europe (Utescher & al. 2007). It is worth mentioning that the aridification of the climate in S Europe during the upper Pliocene (around 3.6 Mya) has led to the diversification of Thymus herba-barona, endemic to the W Mediterranean islands, into three major clades (Molins & al. 2011). At that time, i.e. during the upper Pliocene, the withdrawal of the Pannonian Sea from the central parts of the Balkan Peninsula was in progress, and that event may have triggered the radiation and diversification of numerous taxa within this region. Assuming that the ancestors of Salvia officinalis were already present within the Balkans, the estimated time of divergence of an ancestral S. officinalis gene pool into two distinct lineages of 3.267 Mya (1.584 Mya, 5.247 Mya) would suggest that the withdrawal of the Pannonian Sea may have been an essential milestone during the Pliocene evolutionary history of this species that led to the occurrence of two lineages, currently represented by populations of the sub-structured inland-Adriatic lineage (IAL) and a more homogeneous purely Adriatic lineage (PAL).

Table 5.

Results of the mismatch distribution analyses and neutrality tests in the populations overall, and separately in two Salvia officinalis lineages.

Fig. 5.

Majority rule consensus tree from a strict clock Bayesian analysis. Node heights correspond to median ages; numbers along branches are Bayesian posterior probabilities; identical haplotypes are collapsed as triangles, and their vertical extensions are proportional to the number of individuals.

Pleistocene evolutionary history of Salvia officinalis and indigeneity of its inland populations

The evolutionary history of two Salvia officinalis lineages during the turbulent Pleistocene was distinct and complex. A star-like radiation of closely related haplotypes comprising IAL (A, B, C and D) from the hhl (Fig. 4) would imply past adaptive radiation and expansion of this lineage (see Posada & Crandall 2001; Schaal & al. 2003). Given the estimated divergence time of IAL into inland and Adriatic sub-lineages (1.188 Mya (0.495 Mya, 2.079 Mya)), the diversification of this lineage was most likely triggered by the onset of the Pleistocene, as initially hypothesized by Adamović (1908, 1909). That author assumed that S. officinalis was widely and more continuously distributed throughout the Balkans in the past, and that the repeated southward migrations of continental temperate forests during the Pleistocene led not only to its retreat to inland canyons and gorges but also to its displacement from more continental habitats to sites closer to the Adriatic coast. The described scenario is supported by the current spatial distribution of IAL haplotypes, i.e. by the prevalence of alternative derived but closely related haplotypes in geographically distant IAL sub-lineages (B in inland, C in Adriatic sub-lineage). The ancestral hhl haplotype may have been lost over time during apparently rather dynamic historical changes of IAL populations associated with range alternations (e.g. fragmentation, genetic drift, founder and bottleneck effects). At present, IAL displays gradual size reductions (as inferred from the observed relation of a relatively high haplotype diversity, Hd = 0.426, and a relatively low nucleotide diversity, π = 0.00106) rather than a recent demographic expansion (as demonstrated by both mismatch distribution analysis and neutrality tests).

However, due to the conflicting views with regard to the origin of Salvia officinalis inland disjunctions, it is possible that the inland sub-lineage of IAL is, actually, an artefact. Nonetheless, an important outcome of our work, which undermines both hypotheses on an anthropogenic origin of inland S. officinalis populations, is that all these populations harbour haplotype B, which is rather rare or absent in Adriatic populations. In case of a human-mediated origin of disjunct inland populations of S. officinalis, more abundant Adriatic types would be expected in these populations established either from seeds collected from natural populations by travellers and/or healers and dropped later on during their journeys (a seed transfer scenario) or from seeds/plants collected also from natural populations and transferred to the gardens of medieval monasteries in this region from which they may have expanded throughout the surrounding areas (naturalization of cultivated plants scenario). The latter scenario, however, is less likely because of a low competing ability of S. officinalis (Devetak 1963), which precludes colonization of sites already occupied by more competitive species. In addition, we found exceptionally low seed flow in this species (overall Nm = 0.15), further supporting previous assumptions on rather poor dispersal ability of this species (Devetak 1963). On the other hand, assuming that haplotype B was the prevalent chloroplast type in a previously widely and continuously distributed inland sub-lineage of S. officinalis, this type would be expected in remnant populations upon decreases in range and population sizes, because it is well known that the most abundant types are less sensitive to size reductions (Luikart & al. 1998). However, at this point, we cannot exclude the possibility that haplotype B may be present in other S. officinalis populations not used in the present study (but see later).

Nonetheless, the indigeneity of disjunct inland populations of Salvia officinalis is supported by other lines of evidence as well. Based on biogeographic data, numerous inland gorges throughout the Balkans, including those where this species is currently present, have been delineated as putative glacial refugia (Stevanović & Vasić 1995). Along with S. officinalis, the gorges in question currently host other sub-Mediterranean taxa such as Artemisia alba Turra, Arum italicum Mill., Micromeria cristata (Hampe) Griseb., Ostrya carpinifolia Scop., Paliurus spina-christi Mill., Ruta graveolens L., Ziziphora capitata L. and several endemic taxa as well. Therefore, even if we assume that S. officinalis has been introduced to these sites by humans, it is highly unlikely that all these sub-Mediterranean taxa, few of which have been used by humans, have been established at same sites by human mediation. Furthermore, all disjunct inland S. officinalis populations display characteristics of relict populations that are rather common throughout the entire Mediterranean basin (Stevanović & Vasić 1995; Petit & al. 2005a), such as isolation, small sizes, and low levels of genetic variability.

The Adriatic populations of IAL — LUS and MOR, whose ancestors most likely reached the Adriatic coast prior to the mid-Pleistocene diversification of IAL (1.188 Mya (0.495 Mya, 2.079 Mya)), currently thrive on the SE Adriatic coast where more continuous favourable habitats and suitable climate are available today as they were in the past (Utescher & al. 2007; see also Surina & al. 2011). These currently large, viable and genetically diverse populations (Table 1), however, do not display signs of a sudden demographic expansion, as inferred from the mismatch distribution analysis and neutrality tests. However, it is worth mentioning that the inferences on treatment of these populations as a sub-linage of IAL and not as an independent lineage were not straightforward in SAMOVA, AMOVA (Table 4) and statistical parsimony analysis (Fig. 4) as well as based on population pairwise F_ST values (Table 3). This would suggest that the long-term persistence of a species in a refugial region, such as the Balkans, may result in rather complex contemporary genetic and/or phylogeographic patterns that have to be interpreted with caution.

The other two studied Adriatic populations — ORJ and RUM — represent the second Salvia officinalis lineage — PAL, which is characterized by a relatively high genetic diversity (Hd = 0.403, π = 0.00257). These currently large and viable populations, found within different biogeographic regions and vegetational zones as compared to Adriatic populations of IAL (Table 1), have not experienced a sudden expansion in the recent past. They harbour an exclusive haplotype G, separated by 13 mutational steps from hhl. The occurrence of two rare haplotypes in PAL — E and F, which were both more closely related to hhl and haplotypes from IAL than was haplotype G, supports the common origin of both S. officinalis lineages, and implies that the history of PAL populations was most likely characterized by strong and potentially multiple founder and/or genetic drift effects along with subsequent expansions. This is because such historical population dynamics would enable not only the current prevalence of a derived haplotype G in PAL but also the persistence of haplotypes E and F in that lineage. The finding of haplotype C in one individual from both ORJ and RUM is more puzzling, and may be explained either by immigration or retention of shared polymorphisms. Although an independent estimation from other markers is required to support our assertion about the retention of shared polymorphisms, this assumption is more likely because of the exceptionally low seed flow in S. officinalis (overall Nm = 0.15) and observed barriers between RUM and MOR (but not between LUS and ORJ.

Putative glacial refugia of Salvia officinalis

The apparent phylogeographic structure in Salvia officinalis was detected also by the comparison of N_ST and G_ST values (N_ST = 0.764, SE = 0.046; G_ST = 0.682, SE = 0.066; P = 0.000). The IBD pattern, however, was lacking (r² = 0.003, P = 0.351), at least partly due to the current overlap of IAL and PAL on the SE Adriatic coast, where geographically close populations may not necessarily be genetically close as well. The observed phylogeographic break is indicative of a localized glacial refugium, as demonstrated in other taxa as well (e.g. Gómez & Lunt 2007 and references therein; Horn & al. 2009; see also Médail & Diadema 2009; Surina & al. 2011). It is worth mentioning that the same region, the SE Adriatic coast, was recognized as a glacial refugium of another thermophilous plant species, Edraianthus tenuifolius (Waldst. & Kit.) A. DC. (Campanulaceae), and that the northward post-glacial expansion of this species along the E Adriatic coast advanced from that region (Surina & al. 2011). It is tempting to hypothesize that a similar past population dynamic was present in S. officinalis as well, and that would imply that haplotypes common in the main glacial refugium of this species are to be expected along the E Adriatic coast.

In addition, our molecular and previous biogeographic data (Stevanović & Vasić 1995) support the view that all inland sites with relict populations of Salvia officinalis should be regarded as refugial as well, as already discussed.

Conservation implications

Salvia officinalis is characterized by a rather high level of genetic diversity in the plastid genome, H_T = 0.695, similarly as observed in other native aromatic/medicinal plants (e.g. Tarayre & al. 1997; Yuan & al. 2010; Molins & al. 2011). Given the high genetic diversity at the nuclear DNA level as well (Židovec 2004; Radosavljević & al. 2011, 2012; Greguraš 2013), S. officinalis per se cannot be viewed as an endangered species, and thus, it is expected to keep its reproductive fitness and to evolve and adapt in response to environmental changes over time (Frankham & al. 2002) if further anthropogenic deterioration of habitats and over-exploitation may be prevented or at least limited.

However, given the exceptionally high genetic differentiation at the plastid genome in Salvia officinalis , G_ST = 0.682, which is in the range of values found in Thymus herba-barona fragmentarily distributed in the W Mediterranean islands (Molins & al. 2011), as well as variable levels of genetic diversity at the population level (i.e. genetically variable vs genetically depleted populations), it is necessary to implement conservation measures that would, on the one hand, take into account the historic distinctiveness of populations and, on the other hand, focus on currently most endangered populations, such as those comprising the inland IAL sub-lineage. This is important because in aromatic and/or medicinal plants, genetic distinctiveness of populations may coincide with their differentiation in relation to essential oil content and composition (e.g. Skoula & al. 1999; Trindade & al. 2009). Previous studies revealed that S. officinalis plants from the Sićevo gorge (SIC), recognized as subsp. multiflora (Gajić 1973; Greuter & al. 1986), appear to have a distinct essential oil composition (Couladis & al. 2002), and their loss may have profound effects on future breeding efforts and commercial utilization of this species.

As already mentioned, currently the most threatened Salvia officinalis populations are relict and refugial inland populations, which are small in size, genetically depleted (except SIC), and isolated. Although they have been protected in situ (data from the Institute for Nature Conservation of Serbia), they lack natural regeneration (Stojanović 2002) and are still extensively used as freely available natural resources of these plants. The most threatened is population MIL currently comprising c. 20 individuals. Therefore, we recommend an urgent action for safeguarding relict and refugial populations of S. officinalis via cost-effective long-term storage of seeds (Klumpp 2005), while further actions would require the establishment of ex situ collections. However, the mixing of plants originating from genetically well-differentiated populations should be avoided until it is known whether an outbreeding depression may pose a problem (Fenster & Dudash 1994; Hufford & Mazer 2003). On the other hand, Adriatic populations of S. officinalis, which are genetically diverse, viable and less threatened by harvesting due to their large sizes, should be preserved in situ.

Our results may have implications for systematic and taxonomic considerations in Salvia officinalis, because our molecular data do not support further subdivision of the inland sub-lineage of IAL and the delineation of individuals from Serbia (SIC and MIL) as subsp. multiflora, as previously inferred from morphological data (Gajić 1973; Greuter & al. 1986). Nonetheless, it has been demonstrated recently that DNA barcodes have failed to accurately identify morphologically well-differentiated species because, in cases of low effective mutation rates, even reciprocally monophyletic species may share identical haplotypes (see Velzen & al. 2012 and references therein). On the other hand, a profound distinctiveness of individuals comprising IAL and PAL at the molecular level may reflect divisions at the intraspecific level that have not been evident from morphological data. Those findings, along with views on the presence of typical S. officinalis in Italy (Pignatti 1982; Reales & al. 2004; Conti & al. 2005), which would suggest an amphi-Adriatic distribution of this species, trigger additional questions and necessitate further studies.