Salvia officinalis (Lamiaceae), common or Dalmatian sage, is a Mediterranean aromatic and medicinal plant used in medicine since ancient times. Knowledge on current genetic patterns and genealogical history of its natural populations is required for both breeding efforts and species conservation. We used sequences of two chloroplast intergenic spacers, 3′rps16—5′trnK and rpl32—trnL, from 83 individuals from eight natural populations to distinguish between anthropogenic vs natural origin of four disjunct inland populations found outside of the main Adriatic range of the species. We found seven haplotypes, high total gene diversity (HT = 0.695) and genetic differentiation (GST = 0.682), as well as a phylogeographic structure with two lineages, a sub-structured inland-Adriatic lineage (IAL, comprising inland and Adriatic sub-lineages) and a purely Adriatic lineage (PAL). All four inland and disjunct populations, which comprised the inland sub-lineage of IAL, were almost fixed for a distinct haplotype genealogically closely related to the ancestral haplotype and displayed other features of relict populations. Along with previous biogeographic data and other lines of evidence, assumptions on their anthropogenic origin were rejected. At present, a less diverse IAL (Hd = 0.426, π = 0.00106) and a more diverse PAL (Hd = 0.403, π = 0.00257), whose divergence was dated to the Pliocene (3.267 Mya), do not exhibit signs of recent demographic expansions and overlap on the SE Adriatic coast, a region delineated as the main glacial refugium of S. officinalis. Conservation measures accounting for the historical distinctiveness of populations and focusing on currently the most threatened populations are recommended.
Salvia officinalis L. (Lamiaceae), common or Dalmatian sage, is a Mediterranean aromatic and medicinal herb used in medicine for more than two millennia (Osbaldeston & Wood 2000). Due to the invaluable importance of its bioactive components to humans, it is not surprising that the vast majority of previous studies on this species were focused on the discovery of various bioactive compounds and the characterization of their biological activities (e.g. Karousou & al. 2000; Kintzios 2000; Baričević & al. 2002; Böszörményi & al. 2009 and references therein; Lipman 2009), and that substantial efforts have been made to facilitate selection of plants for commercial production by assessing correlations between the morphological features and/or environmental conditions and essential oil content (e.g. Corsi & Botega 1999; Maksimović & al. 2007; Jug-Dujaković & al. 2012), and to produce cultivars with increased content of relevant bioactive components through breeding and other biotechnology programs (e.g. Kintzios 2000; Nell & al. 2009).
The demands of global markets for Salvia officinalis and its various products, however, are increasing, and they apparently cannot be satisfied by commercial production only. This poses a severe threat to the species because the natural populations have been viewed as a freely available natural resource of these valuable plants and have been exposed to over-exploitation especially in developing countries (Baričevič & al. 2002; Lipman 2009). If such a trend continues, a rapid erosion of natural genetic resources of S. officinalis may be expected in the near future. In order to preserve extant gene pool of this species and to facilitate future breeding efforts, knowledge on levels of its genetic diversity, genetic differentiation and genealogical history is required (Moritz 1994; Frankham & al. 2002).
To date, the levels of genetic diversity and genetic differentiation have been surveyed in cultivated forms and gene bank accessions of Salvia officinalis (e.g. Bazina & al. 2002; Khalil & al. 2005; Bruna & al. 2006; Echeverrigaray & Agostini 2006; Farkas & al. 2008; Böszörményi & al. 2009; Mader & al. 2010). Studies of this kind focusing on natural populations are scarce and involve more or less localized sets of natural populations of this species (Židovec 2004; Radosavljević & al. 2011, 2012; Greguraš 2013) while phylogeographic surveys, which may shed more light on species evolutionary history, have not been undertaken yet. Genetic and/or phylogeographic studies in S. officinalis, however, may be rather challenging, not only in this species but also in other medicinal and/or aromatic plants extensively utilized by humans because of the potential anthropogenic impact on current species distributions (cf. Olsen & Schaal 1999) that can also affect conservation practice (Moritz 1994; Frankham & al. 2002).
According to the Med-Checklist (Greuter & al. 1986), the natural range of Salvia officinalis is limited to the Balkans, a view that we follow in the context of the present study (but see the Discussion for additional remarks). Within the Balkans, S. officinalis occurs continuously along the E Adriatic coast (Devetak 1963; Hedge 1972; Greuter & al. 1986; Karousou & al. 2000; Kintzios 2000). However, it is also found within the central Balkans, in several disjunct and isolated gorges characterized by a specific microclimate and vegetation types and with a high level of specific endemism (Stevanović & Vasić 1995). The origin of these geographically distant inland populations of S. officinalis has been puzzling scientists for more than a century. Adamović (1908, 1909) argued that they represent the remains of a formerly wider continuous range of this species and therefore should be regarded as relict. If this is the case, then more ancient, distinct and scarce genetic types would be expected in these populations as compared to types currently present in Adriatic populations. Furthermore, genetic differentiation in these putative inland glacial refugia may be expected. It is worth noting that S. officinalis individuals from the Sićevo gorge, delineated as S. officinalis subsp. multiflora Gajić based on morphological data (Gajić 1973; Greuter & al. 1986), appear to have a distinct essential oil composition (Couladis & al. 2002).
However, it has also been hypothesized that the disjunct inland populations of Salvia officinalis are not indigenous to the Balkans. Two scenarios have been put forward to explain the occurrence of this Mediterranean herb at several disjunct sites outside of the main range of the species, and both involve human mediation. Given the documented cultivation of S. officinalis in medieval gardens in Europe (Dweck 2000) as well as in numerous monasteries in the Balkans, it has been suggested that these cultivated plants may have escaped from cultivation and expanded throughout the Balkans in the past (cf. Devetak 1963; Dweck 2000). On the other hand, due to the archaeological findings of the major Roman Empire roads Via militaris and Via Egnatia in the Sićevo gorge and the Crni Drim gorge, Mt Jablanica, respectively (cf. Syme & Birley 1999), in which S. officinalis is currently present, it has been postulated that this species may have established itself at those sites from seeds unintentionally dropped by traders and/or healers in transit. The status of S. officinalis as indigenous to the central Balkans would be undermined if these inland populations were to harbour derived genetic variants found to be common in Adriatic populations.
In the context of our phylogeographic study, we characterize the conflicting inferences made with regard to the origin of inland disjunctions of Salvia officinalis, localize glacial refugia of this species and provide guidelines for species conservation. Given the probable maternal plastid inheritance in all Lamiales sampled by Corriveau & Coleman (1988) and the preferential utility of such genomes in phylogeographic surveys in plants (Schaal & al. 1998), we tested two chloroplast (cpDNA) intergenic spacers in 83 S. officinalis individuals originating from four natural populations from a more continuous, SE Adriatic part of the species range, as well as from four disjunct inland sites in the central Balkans. To accomplish these objectives, we assessed: (1) levels of cpDNA genetic diversity in S. officinalis populations and species per se; (2) genetic structuring and potential phylogeographic breaks; (3) levels of seed flow among populations; and (4) divergence patterns and their time frames, which have led to the establishment of the detected contemporary genetic/phylogeographic structuring in this species.
Material and methods
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.
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 (NH4)2SO4 (Fermentas, Vilnius, Lithuania); 2.5 mM MgCL2; 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 (Hs) and total gene diversity (HT) using PERMUTE (Pons & Petit 1996).
Genetic differentiation, spatial genetic structuring and phylogeographic analyses
Pairwise population differentiation statistics (FST) 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/FST — 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.
We tested for a correlation between linear pairwise population genetic distances, as measured by FST, 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 FCT and minimizes FSC. We ran 100 simulated annealing processes for K ranging from K = 2 to K = 6, and recorded the configuration with the largest FCT values and lowest FSC as the best grouping of populations. In addition, SAMOVA was used to identify barriers between populations.
A genetic distance matrix of pairwise population FST 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 (FCT), among population within each lineage (FSC), and within populations (FST), 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 FCT and the lowest FSC values, similarly as in SAMOVA, and can reveal homogeneity of detected groupings.
We used CONTRIB (Petit & al. 1998) to calculate genetic differentiation for unordered (GST) and ordered alleles (NST) and to test the phylogeographic signal across the species range by comparing these estimates. A phylogeographic structure is implied if NST is higher than GST because closely related haplotypes are more often found in the same geographic area than would be expected by chance. A statistical significance between GST and NST 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).
Sampling localities of Salvia officinalis populations and genetic variation in two cpDNA regions.
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.
The neutrality Fu's Fs 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, Hs, was 0.221 (SE = 0.074), and total gene diversity, HT, of 0.695 (SE = 0.120).
Genetic differentiation, spatial genetic structuring and phylogeographic analyses
Pairwise population FST 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). FST 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 FST values was > 1.00 among populations which were not significantly differentiated and < 1.00 between highly and significantly differentiated pairs of populations (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 FST values and the logarithm of the geographic distances between populations was positive, it was low and insignificant (r2 = 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, FCT value was high and significant (0.799, P = 0.03, Table 4) and FSC 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, FCT value slightly increased and remained significant (0.813, P = 0.003, Table 4) while FSC 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.
Summary of AMOVA and SAMOVA statistics for alternative groupings of populations.
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 FCT values was negligible (0.815, 0.818 and 0.800, respectively) as well as the decrease of FSC 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 FCT values were high but at K = 2, FCT 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%, FSC = 0.012, P > 0.05) than under the assumption of two lineages (IAL and PAL) (6.77%, FSC = 0.336, P = 0.00). Also, lower and non-significant FST 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: FST = 0.697, P = 0; PAL: FST = 0.004, P > 0.05).
The comparison between NST and GST showed a significant difference (NST = 0.764, SE = 0.046; GST = 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.
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 Fs 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).
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 km2, 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).
Results of the mismatch distribution analyses and neutrality tests in the populations overall, and separately in two Salvia officinalis lineages.
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 FST 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 NST and GST values (NST = 0.764, SE = 0.046; GST = 0.682, SE = 0.066; P = 0.000). The IBD pattern, however, was lacking (r2 = 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.
Salvia officinalis is characterized by a rather high level of genetic diversity in the plastid genome, HT = 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 , GST = 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.
This work was financially supported by the Ministry of Education and Science of the Republic of Serbia, Research Grants 173021 and 173030. The authors would like to thank D. Lakušić for constructive comments that helped to improve an earlier version of the manuscript, J. Šinžar-Sekulić for help in assembling the distribution map based on literature and herbarium data, and J. D. Lockwood for language editing. Peter B. Phillipson (Muséum National d'Histoire Naturelle, Paris) and Yasaman Salmaki (University of Tehran) are also thanked for their reviews of an earlier version of the manuscript.
- L. Adamović 1908: Die Bedeutung des Vorkommens der Salbei in Serbien. — Bot. Jahrb. Syst. 41: 175–179. Google Scholar
- L. Adamović 1909: Die Vegetationsverhältnisse der Balkanländer (Mösische Länder): umfassend Serbien, Altserbien, Bulgarien, Ostrumelien, Nordthrakien und Nordmazedonien. — Leipzig: W. Engelmann. Google Scholar
- D. Baričevič , J. Bernath , L. Maggioni & E. Lipman 2002: Report of a working group on medicinal and aromatic plants. — Rome: International Plant Genetic Resources Institute. Google Scholar
- S. Bruna , A. Giovannini , L. De Benedetti , M. C. Principato & B. Ruffoni 2006: Molecular analysis of Salvia spp. through RAPD markers. — Acta Hort. 723: 157–160. Google Scholar
- F. Conti , G. Abbate , A. Alessandrini & C. Blasi 2005: An annotated check-list of the italian vascular flora. — Roma: Edizioni Palombi. Google Scholar
- Z. Devetak 1963: Prilog poznavanju kadulje u Jugoslaviji. — Radovi Poljoprivrednog Fakulteta u Beogradu XII: 241–261. Google Scholar
- A. C. Dweck 2000: The folklore and cosmetic use of various Salvia species. — Pp. 1–25 in: S. E. Kintzios (ed.), Sage. The genus Salvia. [Medicinal and Aromatic Plants — Industrial Profiles 14]. — Amsterdam: Harwood Academic. Google Scholar
- S. Echeverrigaray & G. Agostini 2006: Genetic relationships between commercial cultivars and Brazilian accessions of Salvia officinalis L. based on RAPD markers. — Rev. Bras. Pl. Med. 8: 13–17. Google Scholar
- L. Excoffier , G. Laval & S. Schneider 2005: Arlequin ver. 3.0: an integrated software package for population genetics data analysis. — Evol. Bioinform. Online 1: 47–50. Google Scholar
- L. Excoffier , P. E. Smouse & J. M. Quattro 1992: Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. — Genetics 131: 479–491. Google Scholar
- A. Farkas , N. Papp , G. Horvath , T. S. Nemeth , I. Szabo & T. Nemeth 2008: RAPD based genetic diversity among Salvia officinalis L. populations. — Farmacia 56: 339–343. Google Scholar
- Y.-X. Fu 1997: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. — Genetics 147: 915–925. Google Scholar
- M. R. Gajić 1973: Varijabilnost žalfije (Salvia officinalis L.). Nova podvrsta (S. officinalis L. subsp. multiflora Gajić subsp. nova). — Glasn. Prir. Muz. Beogradu, C, 7: 29–51. Google Scholar
- D. Graur & W. H. Li 2000: Fundamentals of molecular evolution. — Sunderland: Sinauer Associates. Google Scholar
- D. Greguraš 2013: Genetic diversity and population structure of Dalmatian sage (Salvia officinalis L.). — Dissertation, Faculty of Agriculture, University of Zagreb. Google Scholar
- W. Greuter , H. M. Burdet & D. Long 1986: Med-Checklist: a critical inventory of vascular plants of the circum-mediterranean countries. 3. Dicotyledones (Convolvulaceae-Labiatae). — Genève: Conservatoire et Jardin botaniques de la Ville de Genève; Berlin: Secretariat Med-Checklist, Botanischer Garten & Botanisches Museum Berlin-Dahlem. Google Scholar
- I. C. Hedge 1972: Salvia L. — Pp. 188–192 in: T. G. Tutin , V. H. Heywood , N. A. Burges , D. M. Moore , D. H. Valentine , S. M. Walters & D. A. Webb (ed.), Flora europaea 3. — Cambridge: University Press. Google Scholar
- R. Karousou , E. Hanlidou & S. Kokkini 2000: The sage plants of Greece: distribution and infraspecific variation. — Pp. 27–53 in: S. E. Kintzios (ed.), Sage. The genus Salvia. [Medicinal and Aromatic Plants — Industrial Profiles 14]. — Amsterdam: Harwood Academic. Google Scholar
- S. E. Kintzios 2000: Sage. The genus Salvia. [Medicinal and Aromatic Plants — Industrial Profiles 14]. — Amsterdam: Harwood Academic. Google Scholar
- R. T. Klumpp 2005: Seed and pollen storage: European focus. — Pp. 601–622 in: T. Geburek & J. Turok (ed.), Conservation and management of forest genetic resources in Europe. — Zvolen: Arbora Publishers. Google Scholar
- E. Lipman 2009: Report of a working group on medicinal and aromatic plants. — Rome: Bioversity International. Google Scholar
- N. Mantel 1967: The detection of disease clustering and generalize regression approach. — Cancer Res. 27: 209–220. Google Scholar
- A. Molins , G. Bacchetta , M. Rosato , J. A. Rossello & M. Mayol 2011: Molecular phylogeography of Thymus herba-barona (Lamiaceae): Insight into the evolutionary history of the flora of the western Mediterranean islands. — Taxon 60: 1295–1305. Google Scholar
- T. A. Osbaldeston & R. P. A. Wood 2000: Dioscorides de materia medica: being an herbal with many other medicinal materials written in Greek in the first century of the common era; a new indexed version in modern English. — Johannesburg: Ibidis Press. Google Scholar
- N. K. Pantić 1984: On the evolution of land flora based on plant fossils from the territory of Serbia. — Pp. 191–246 in: M Jankovic (ed.), Végétation de la République Socialiste de Serbie. Tome 1. Partie générale. — Beograd: Académie Serbe des Sciences et des Arts. Google Scholar
- S. Pignatti 1982: Flora d'Italia 1-3. Bologna: Edagricole. Google Scholar
- O. Pons & R. J. Petit 1996: Measuring and testing genetic differentiation with ordered versus unordered alleles. — Genetics 144: 1237–1245. Google Scholar
- I. Radosavljević , Z. Šatović , J. Jakše , B. Javornik , D. Greguraš , M. Jug-Dujaković & Z. Liber 2012: Development of new microsatellite markers for Salvia officinalis L. and its potential use in conservation-genetic studies of narrow endemic Salvia brachyodon Vandas. — Int. J. Molec. Sci. 13: 12082–12093 Google Scholar
- A. R. Rogers & H. Harpending 1992: Population growth makes waves in the distribution of pairwise genetic differences. — Molec. Biol. Evol. 9: 552–569. Google Scholar
- I. R. Savić 2008: Diversification of the Balkan fauna: its origin, historical development and present status. — Pp. 57–78 in: S. E. Makarov & R. N. Dimitrijević (ed.), Advances in arachnology and developmental biology: papers dedicated to Professor Božidar P. M. Ćurčić, Monographs 12. — Belgrade: Institute of Zoology, University of Belgrade; Sofia: BAS; Vienna: Faculty of Life Sciences, University of Vienna; Belgrade: SASA & UNESCO MAB Serbia. Google Scholar
- S. Stanković 1960: The Balkan lake Ohrid and its living world. — Den Haag: Dr. W. Junk. Google Scholar
- V. Stevanović & V. Vasić 1995: Biodiverzitet Jugoslavije sa pregledom vrsta od medunarodnog značaja. — Beograd: Biološki fakultet i Ecolibri. Google Scholar
- D. Stojanović 2002: Infraspecijska varijabilnost Salvia officinalis L. (Lamiaceae) u Jugoslaviji. — Master's thesis, University of Belgrade. Google Scholar
- R. Syme & A. R. Birley 1999: The provincial at Rome: and, Rome and the Balkans 80BC-AD14. — Exeter: University Press. Google Scholar
- F. Tajima 1989: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. — Genetics 123: 585–595. Google Scholar
- M. Tarayre , P. Saumitou-Laprade , J. Cuguen , D. Couvet & J. D. Thompson 1997: The spatial genetic structure of cytoplasmic (cpDNA) and nuclear (allozyme) markers within and among populations of the gynodioecious Thymus vulgaris (Labiatae) in southern France. — Amer. J. Bot. 84: 1675–1684 Google Scholar
- A. R. Templeton , K. A. Crandall & C. F. Sing 1992: A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. — Genetics 132: 619–6. Google Scholar
- H. Trindade , M. M. Costa , S. B. Lima , L. G. Pedro , A. C. Figueiredo & J. G. Barroso 2009: A combined approach using RAPD, ISSR and volatile analysis for the characterization of Thymus caespitititus from Flores, Corvo and Graciosa islands (Azores, Portugal). — Biochem. Syst. Ecol. 37: 670–677. Google Scholar
- V. Židovec 2004: Varijabilnost prirodnih populacija mirisave kadulje (Salvia officinalis L.). — Dissertation. Faculty of Agriculture, University of Zagreb. Google Scholar