Development of New Microsatellite Markers for Salvia officinalis L. and Its Potential Use in Conservation-Genetic Studies of Narrow Endemic Salvia brachyodon Vandas

Nine new microsatellite markers (SSR) were isolated from Salvia officinalis L. A total of 125 alleles, with 8 to 21 alleles per locus, were detected in a natural population from the east Adriatic coast. The observed heterozygosity, expected heterozygosity, and polymorphic information content ranged from 0.46 to 0.83, 0.73 to 0.93 and 0.70 to 0.92, respectively. New microsatellite markers, as well as previously published markers, were tested for cross-amplification in Salvia brachyodon Vandas, a narrow endemic species known to be present in only two localities on the Balkan Peninsula. Out of 30 microsatellite markers tested on the natural S. brachyodon population, 15 were successfully amplified. To obtain evidence of recent bottleneck events in the populations of both species, observed genetic diversity (HE) was compared to the expected genetic diversity at mutation-drift equilibrium (HEQ) and calculated from the observed number of alleles using a two-phased mutation model (TPM). Recent bottleneck events were detected only in the S. brachyodon population. This result suggests the need to reconsider the current threat category of this endemic species.

The primary goal of this research was to identify new microsatellite markers for Dalmatian sage and to establish, together with previously developed microsatellite markers, a set of SSR markers for future population genetics studies of this species. The secondary goal was to examine the possibility of using microsatellite markers developed for Dalmatian sage in population genetics studies of short-tooth sage. This second objective is highly dependent on the total number of available markers. To meet both goals and to demonstrate that the loci contain sufficient variation for individual discrimination, natural populations of each species from the Pelješac peninsula were studied.

Development of New Microsatellite Markers for Dalmatian Sage
In total, 3840 colonies were screened for dinucleotide repeats in Dalmatian sage (Salvia officinalis L.). After removing low quality reads, 235 unique sequences remained in a total length of 107,875 bp. Out of 235 unique sequences, 224 contained GA or GT microsatellite repeats. High quality PCR primer pairs were designed for 15 dinucleotide microsatellite loci. Eleven of these were polymorphic, while four were either monomorphic or did not amplify at all. Finally, nine primer pairs had suitable amplification patterns and signal intensity and were used to screen 25 individuals representing a natural Dalmatian sage population ( Table 1). The screen resulted in a total of 125 alleles, 8 to 21 alleles per locus, with an observed heterozygosity from 0.46 to 0.83, and an expected heterozygosity from 0.73 to 0.93. The DNA sequences of these microsatellite loci were deposited into GenBank under accession numbers JX440363 to JX440371. Eight out of nine microsatellites showed a high polymorphic information content (PIC) of more than 0.75. It should be noted that the PIC value of the remaining locus (SoUZ021) was rather high (0.70), indicating that all nine loci could be very useful in assessing the genetic diversity and population structure of Dalmatian sage. Three out of nine newly developed microsatellite markers (SoUZ022, SoUZ025 and SoUZ027) showed significant deviations from Hardy-Weinberg expectations (HWE) after application of the sequential Bonferroni corrections. These three loci also exhibited an overall excess of homozygotes and null allele frequencies using Brookfield's formula [20]. They varied from 0.14 (SoUZ025) to 0.21 (SoUZ027). However, bearing in mind the distribution range of the tested S. officinalis population as well as the range of altitudes of sampled individuals, this heterozygote deficiency is more likely due to population structure than to locus-specific phenomenon (e.g., scoring error or null alleles). In accordance with this opinion, we recommend a more detailed sampling in future population genetic studies of Dalmatian sage than was performed in this study. If using this sampling design results in the same loci continuing to show an overall excess of homozygotes and null allele frequency, corrections for null alleles or their exclusion from the study can be used as a last resort [21][22][23].

Cross-Amplification in Narrow Endemic Salvia brachyodon
Including previously published di-and tri-nucleotide microsatellite loci [24,25] and the primers published in this study, we were able to test 29 Dalmatian sage microsatellite markers for cross-amplification in short-tooth sage. The amplification rate was 52%. The 15 successfully amplified microsatellite markers were used for in-depth analysis in natural populations of both species ( Table 2). The development of new microsatellite markers, which were described earlier in this paper, proved to be entirely justified because as many as four of these markers were among the 15 polymorphic markers that were successfully amplified in short-tooth sage. A total of 87 alleles were observed across 15 loci. The number of alleles per locus ranged from 3 to 9, the observed heterozygosity ranged from 0.33 to 0.92, and the expected heterozygosity ranged from 0.33 to 0.84. Only two microsatellite loci (SoUZ026 and SoUZ002) had low polymorphic information contents (PIC) of 0.31 and 0.47. One of the 15 microsatellite loci exhibited significant deviations from HWE (SoUZ020) and the presence of null alleles with a frequency of 0.18. Five out of the 105 tests for linkage disequilibrium were significant (p < 0.01) after applying sequential Bonferroni corrections (SoUZ006/SoUZ007, SoUZ009/SoUZ011; SoUZ014/SoUZ009; SoUZ006/SoUZ009 and SoUZ014/SoUZ020). Deviations from HWE and linkage disequilibrium are the result of primer-site mismatch (null alleles), which are common in cross-amplified species, or as a consequence of specific population structures (e.g., clonality). It is noteworthy that during fieldwork, dense patches of S. brachyodon individuals and stolons at the soil surface or below ground were observed (Figure 1c). Therefore, it is very likely that clonal individuals exist in S. brachyodon populations and that a significant result for linkage disequilibrium is a consequence. However, in our sample set of 25 individual short-tooth sage plants, none were genetically identical. More intensive population genetics research and tests for clonal propagation [26,27] will allow us to elaborate on short-tooth sage reproduction and offer a better explanation for deviations from HWE and linkage disequilibrium observed in this study.

Population Genetics Parameters and Structures of Natural Populations of Two Closely Related Species
Based on our descriptive statistics (Table 2), we can conclude that the number of alleles as well as the expected heterozygosity or genetic diversity were prominently higher in Dalmatian than in short-tooth sage. Only 29 of 234 detected alleles were common to both species studied ( Table 2). Because microsatellites are known to have very rapid evolutionary rates, prominent differences between Dalmatian and short-tooth sage in the length of alleles at each microsatellite locus were expected. Therefore, the alleles that are common to both species may be the same due to coincidence (identity-in-state) rather than because of a common origin (identity-by-descent).  This study is among those that have shown that rare species have less genetic variation than widespread species [28,29]. However, by virtue of selecting the most polymorphic microsatellites, number of alleles tends to be higher in the species from which they were originally developed [30][31][32] and it is impossible to estimate to which extent cross-amplification procedure itself also contributed to the reduced genetic variation observed in the population of short-tooth sage. If the occurrence of null alleles is uncommon, the microsatellites successfully amplified in a related species proved to be a very useful tool in population genetic studies. The greatest benefit of this study was the opportunity to explore population genetics phenomenon such as genetic bottleneck in rare and widespread congeners [33,34]. In this case, the differences in microsatellite variability between species do not influence the outcome of the analysis. The observed genetic diversity (H E ) is compared to the expected genetic diversity at mutation-drift equilibrium (H EQ ) in the same way regardless of the actual number of alleles at a locus and there is no advantage in using highly variable markers, as shown by the power analysis [35]. Tests for evidence of a recent bottleneck based on microsatellite loci and assuming the two-phased mutation model (TPM) were applied to both Dalmatian and short-tooth sage populations. A statistically significant departure from mutation-drift equilibrium was detected in short-tooth sage (Wilcoxon test; p = 0.02), suggesting that this population underwent a recent bottleneck that reduced its genetic diversity. The results of the same test, when applied to Dalmatian sage and based on 29 (all loci) or 15 (only those loci that were amplified in short-tooth sage) microsatellite loci, were not significant in either case (p = 0.99 and 0.81, respectively). This result was expected for Dalmatian sage population because this is a large population with an extended distribution range.
According to recent data, short-tooth sage from the Pelješac peninsula is considered a near threatened (NT) species of Croatian flora [16]. According to the IUCN Standards and Petitions Subcommittee [36], near threatened species are close to being qualified as vulnerable. Evidence of a genetic bottleneck in the population on the Pelješac peninsula should prompt new testing of IUCN quantitative criteria used in determining whether a species is truly endangered. If these criteria prove that short-tooth sage is truly threatened, then we hope that this study will contribute to better protection.

Experimental Section
This study was carried out on 25 Dalmatian sage and 25 short-tooth sage plants from the natural populations originating on the Pelješac peninsula (Croatia). New microsatellites for S. officinalis L. were developed from genomic DNA libraries enriched for GA and GT repeats as described earlier [37], but with several modifications.
Genomic DNA samples were extracted from dried leaves using the GenElute Plant Genomic DNA Miniprep Kit (Sigma-Aldrich, St. Louis, MO, USA). Nine restriction enzymes were used for genomic DNA digestion (HaeIII, MseI, Sau3AI, RsaI, AluI, HinfI, EcoRV, BglII and EcoRI) (New Englands Biolabs, Ipswich, MA, USA). Single-stranded overhangs of restriction fragments were removed using mung bean nuclease and dephosphorylated using calf intestinal alkaline phosphatase (CIP) (New England Biolabs, Ipswich, MA, USA). Phosphorylated linkers were prepared from SNXfor and SNXrev primers using T4 polynucleotide kinase (Thermo Fisher Scientific Inc., Waltham, MA, USA) [38]. Ligation of linkers to DNA fragments was performed by combining double-stranded SNX linkers, DNA fragments, XmnI restriction enzyme and T4 DNA ligase (New England Biolabs, Ipswich, MA, USA). Long (GA) n and (GT) n probes were constructed in a PCR extension reaction, followed by their attachment to small (5 × 5 mm) pieces of nylon membrane (Nytran ® Super Charge, Schleicher & Schuell BioScience GmbH, Dassel, Germany) and overnight hybridization of DNA fragments containing microsatellite regions. After the nylon membranes were rinsed, the microsatellite fragments were ligated into the pGEM-T Easy Vector (Promega Corporation, Madison, WI, USA), and heat-shock transformation into XL10-Gold Ultracomponent Cells (Agilent Technologies, Stratagene, La Jolla, CA, USA) was performed. The resulting culture was spread on LB-agar plates containing ampicillin, IPTG and X-gal. After overnight incubation, white bacterial colonies were transferred by toothpick into 384-well plates containing Luria-Bertani (LB) freezing media (LB broth + 13 mM KH 2 PO 4 , 36 mM K 2 HPO 4 , 1.7 mM sodium citrate, 6.8 mM (NH 4 ) 2 SO 4 , 4.4% v/v glycerol) for long-term storage. Libraries were transferred onto nylon membranes and screened by Southern analysis using Cy5-and Cy3-labeled 30 bp oligonucleotides with microsatellite core repeats. Positives were detected after stringent washing by scanning the blots using an Ettan DIGE Imager (GE Healthcare Biosciences, Pittsburgh, PA, USA).
Positive clones were randomly selected from the libraries and used for plasmid isolation (Wizard Plus SV Minipreps, Promega Corporation, Madison, WI, USA). Sequencing of plasmid isolates was performed by means of T7 and SP6 universal primers using Big Dye chemistry and an ABI 3730XL analyzer (Applied Biosystems, Foster City, CA, USA). Sequences were edited and assembled using CodonCode Aligner software version 2.0.6 (CodonCode Corporation, Dedham, MA, USA). Microsatellite repeats in sequences were localized by MISA PERL SCRIPT [39]. PCR primer pairs flanking microsatellite repeats were designed using the PRIMER 3 program [40]. Because some of the SSR markers could have been monomorphic or might not have amplified well, new microsatellite PCR primers were first tested on five randomly chosen Dalmatian sage DNAs. Only polymorphic SSR markers with good amplification were tested on the complete set of DNA samples from both natural populations studied using a tailed primer protocol [41]. PCR amplification was performed on the GenAmp ® PCR System 9700 (Applied Biosystems, Foster City, CA, USA) using a two-step protocol with an initial touchdown cycle. The cycling condition were as follows: For each microsatellite locus, the average number of alleles per locus (N a ), the observed heterozygosity (H O ), the expected heterozygosity or genetic diversity (H E ), and the polymorphism information content (PIC) [42] were calculated using PowerMarker V3.23 [43]. GENEPOP version 3.4 [44] was used to test genotypic frequencies for conformance to Hardy-Weinberg expectations (HWE) and to test the loci for gametic disequilibrium (LD). Sequential Bonferroni adjustments [45] were applied to correct for the effect of multiple tests using SAS release 8.02 (SAS Institute Inc., Cary, NC, USA). Each locus was evaluated for the presence of null alleles, scoring errors, and allelic dropout using Micro-Checker version 2.2.3 [46]. The program BOTTLENECK version 1.2.02 [35] was used to test for evidence of recent bottleneck events in the populations of both species on the basis of this theoretical expectation [47]. The observed genetic diversity (H E ) was compared to the expected genetic diversity at mutation-drift equilibrium (H EQ ) and calculated from the observed number of alleles under the intermediate Two-Phase Model (TPM), assuming 30% multistep changes. The Two-Phase Model was applied because it has been shown to be the most appropriate for microsatellite DNA data [48]. Based on the number of loci in our dataset, the Wilcoxon signed-rank test [49] was chosen for the statistical analysis of heterozygote excess or deficiency, as recommended by Piry et al. [47].

Conclusions
Nine new microsatellite markers (SSR) were isolated from Dalmatian sage (Salvia officinalis L.). The observed parameters indicate that all loci may be very useful for assessing genetic diversity and population structures of Dalmatian sage. New microsatellite markers, as well as previously published markers, were tested for cross-amplification in short-tooth sage (Salvia brachyodon Vandas). The amplification rate was 52%. The greatest benefit of this study was the opportunity to compare genetic variation in rare and widespread congeners. This study is among those that have shown that rare species have less genetic variation than widespread species. Recent bottleneck events detected in the short-tooth sage population using a two-phased mutation model (TPM) highlight the need to reconsider the current threat category of this endemic species.