Development and characterization of 16 polymorphic microsatellite markers from Taiwan cow-tail fir, Keteleeria davidiana var. formosana (Pinaceae) and cross-species amplification in other Keteleeria taxa

Background Keteleeria davidiana var. formosana (Pinaceae), Taiwan cow-tail fir, is an endangered species listed on the IUCN Red List of Threatened Species and only two populations remain, both on the Taiwan Island. Sixteen polymorphic microsatellite loci were developed in an endangered and endemic gymnosperm species, Keteleeria davidiana var. formosana, and were tested in an additional 6 taxa, K. davidiana var. calcarea, K. davidiana var. chienpeii, K. evelyniana, K. fortunei, K. fortunei var. cyclolepis, and K. pubescens, to evaluate the genetic variation available for conservation management and to reconstruct the phylogeographic patterns of this ancient lineage. Findings Polymorphic primer sets were developed from K. davidiana var. formosana using the modified AFLP and magnetic bead enrichment method. The number of alleles ranged from 3 to 16, with the observed heterozygosity ranging from 0.28 to 1.00. All of the loci were found to be interspecifically amplifiable. Conclusions These polymorphic and transferable loci will be potentially useful for future studies that will focus on identifying distinct evolutionary units within species and establishing the phylogeographic patterns and the process of speciation among closely related species.


Findings
Background For conservation management, the evaluation of genetic diversity and the identification of distinct evolutionary units in endangered taxa should allow for the development of more efficient conservation strategies [1][2][3]. Over the past century, human activities associated with the lumber industry and agriculture have destroyed and fragmented forest habitats at altitudes of 500-2500 m on Taiwan Island. Several gymnosperm species have decreased in population size and are now under threat, including Cycas taitungensis [4][5][6], Podocarpus nakaii [1], Amentotaxus formosana [2], and Keteleeria davidiana var. formosana [7]. Extreme population declines can affect the genetic diversity of a species [8]. The loss of genetic diversity carries serious evolutionary concerns regarding both the ability to adapt to changes on oceanic islands and longer-term interactions with other organisms [9]. Thus, molecular markers are used to evaluate endangered species in greater depth, with the purpose of applying this genetic diversity information to the conservation and restoration of biodiversity [10].
Keteleeria davidiana var. formosana (Pinaceae), Taiwan cow-tail fir, is an endangered species listed on the IUCN Red List of Threatened Species [11] and derived from K. davidiana, including var. davidiana and var. calcarea, occurred in China by geographical disjunction and leaf scars obscurely or obviously protruding on branchlets [12]. The taxa has only two wild populations, the Dawu Taiwan Keteleeria Forest Reserve (Dawu) and the Pinglin Taiwan Keteleeria Nature Reserve (Pinglin), which are located in Northern and Southern Taiwan, respectively [13]. A previous inventory study has found that less than a hundred individuals remain in the wild, few seedlings are found in the field, and adult trees are vulnerable to lightning damage or death [14]. Additionally, continuous logging over the past century in the low-elevation broadleaf forests has restricted K. davidiana var. formosana to narrow habitats and has decreased its population size [15]. Therefore, the conservation strategies for Keteleeria davidiana var. formosana should include in situ protection and management in conjunction with ex situ approaches such as clonal orchards. Hence, the population genetics of this species should be evaluated before in situ and ex situ conservation management begins. In this study, we developed polymorphic microsatellite loci from Taiwan cow-tail fir to estimate the population structure within the species and to identify distinct evolutionary units within the population. In addition, the transferability of the microsatellite primers developed in this study was tested in other taxa, specifically, three species and three varieties of the Keteleeria genus.

Sampling and DNA extractions
Twenty-five individuals were sampled from each of the two remaining populations of K. davidiana var. formosana from Dawu and Pinglin, Taiwan (Table 1). To test the transferability of the markers, two to three individuals of three species, K. pubescens, K. fortune, and K. evelyniana, and three varieties, K. fortunei var. cyclolepis, var. calcarea, and K. davidiana var. chienpeii, were sampled from the field. Permission for tissue collection was obtained from the Council of Agriculture, Republic of China, authorities. All taxa of the Keteleeria genus are diploid with 24 chromosomes [13]. The sample size, location, and storage herbarium for the voucher specimens are listed in Table 1. Total genomic DNA was extracted from silica-dried leaf powder using the Plant Genomic DNA Extraction Kit (RBC Bioscience, Taipei, Taiwan).

Isolation of microsatellite DNA loci and identification
To isolate the microsatellites efficiently, we used the modified AFLP [16] and magnetic bead enrichment method [1,17,18] to select microsatellite loci. Genomic DNA from one individual of K. davidiana var. formosana was treated with the restriction enzyme MseI (Promega, Madison, Wisconsin, USA) and separated using 1% Nusieve® 3:1 agarose gel (FMC Bio Products, Rockland, ME, USA) electrophoresis. The cut DNA fragments ranged from 400 to 1000 bps, and the purified fragments were extracted from the agarose gel using the HiYield™ Gel PCR DNA Fragments Extraction Kit (RBC Bioscience). Purified DNA fragments were ligated to a double-stranded MseI-adaptor (complementary oligo A: 5′-TACTCAG GACTCAT-3′; 5′phosphorylated oligo B: 5′-GACGAT GAGTCCTGAG-3′) using the Quick Ligation™ Kit (New England Biolabs, MA, USA) at 25°C for 5 minutes. The ligation product was used as the template DNA for the enrichment of the partial genomic library, and the adaptorspecific primers, named Mse I-N (5′-GATGAGTCCT GAGTAAN-3′), were used to perform 20 cycles of prehybridization PCR amplification. The PCR mixture contained 20 ng template DNA, 10 pmol adapter-specific primer, 2 μL 10 × reaction buffer, 2 mM dNTP mix, 2 mM MgCl 2 , 0.5 U Taq DNA polymerase (Promega), and sterile water was added to reach a total volume of 20 μL. The PCR protocol was set at 94°C for 5 min, followed by 18 cycles of 94°C for 30 s, 53°C for 1 min, and 72°C for 1 min, using a Labnet MultiGene 96-well Gradient Thermal Cycler (Labnet, Edison, New Jersey, USA). The amplicons were denatured and hybridized to four separated 5′-biotinylated oligonucleotide probes, (AG) 15 , (AC) 15 , (TCC) 10 , and (TTG) 10 , in 250 μL of hybridization solution at 68°C for 1 h, and 1 mg of Streptavidin MagneSphere Paramagnetic Particles (Promega) was added to the mixture solution to capture the hybridizations at 42°C for 2 h. The four enriched microsatellite DNA fragments were washed with high-salt and low-salt solutions and then used as templates for 25 cycles of PCR amplification using the adaptor-specific primers; the amplification protocol was same as that used for the prehybridization PCR. The PCR products were purified using the HiYield™ Gel PCR DNA Fragments Extraction Kit (RBC Bioscience) and then cloned using the pGEM®-T Easy Vector System (Promega). The plasmid DNAs that were purified from white colonies were digested by the restriction enzyme EcoRI (Promega) and screened using 1% agarose electrophoresis to determine the size of the inserted DNA fragments.  [19], and a pair of specific primers for each microsatellite locus was designed using FastPCR software version 6.4.18 [20].

DNA amplification and genotyping
The optimal annealing temperature was evaluated by gradient PCR using a Labnet MultiGene 96-well Gradient Thermal Cycler (Labnet) for a temperature range from 50°C to 65°C. For each Keteleeria taxon, two to three individuals were evaluated for the optimal annealing temperature. The polymorphism of the two remaining populations of K. davidiana var. formosana was evaluated using 25 individuals from each population (Table 1). To test the optimal annealing temperature, a 20 μL reaction mixture containing 20 ng template DNA, 0.2 μM each of the reverse and forward primers, 2 μL 10 × reaction buffer, 2 mM dNTP mix, 2 mM MgCl 2 , 0.5 U Taq DNA polymerase (Promega), and sterile water was amplified using a Labnet MultiGene 96well Gradient Thermal Cycler (Labnet). The PCR program was set at 94°C for 5 min, followed by 30 cycles of 94°C for 40 s, a temperature gradient ranging from 50 to 65°C for 60 s, and 72°C for 60 s, and a final extension of 72°C for 10 minutes [21]. Amplicons were checked by 1% agarose electrophoresis to isolate the target DNA bands, which were confirmed by sequencing [18]. To examine the genetic polymorphisms, the PCR cocktail was amplified using the previously described program with the temperature gradient replaced by the optimal annealing temperature (Ta). The amplicons were checked by 1% agarose electrophoresis to isolate the target DNA bands, which were confirmed by sequencing [18]. The amplicons were separated by electrophoresis on a 10% polyacrylamide gel (acrylamide: bisacrylamide 29: 1, 80 V for 14-16 hours) using a 25 or 50 bp DNA Step Ladder (Promega) to determine the allele size. The bands were then imaged under UV light using the Flo Gel FGIS-3 fluorescent gel image system (Top BIO Co., Taipei, Taiwan), and the sizes of the PCR products were identified using Quantity One software version 4.62 (Bio-Rad Laboratories, Hercules, California, USA).

Data analysis
The genetic variation indices, including the number of alleles (Na), the number of effective alleles (Ne), the observed heterozygosity (Ho) and expected heterozygosity (He), Shannon's information index (H), and fixation index (F IS ) were estimated using GenAlEx version 6.5 [22]. The Hardy-Weinberg equilibrium (H WE ) was tested using Arlequin software version 3.5.1.2 [23].

Development of polymorphic microsatellite markers
A total of 392 microsatellite loci containing repeat motifs were discovered, with maximum and minimum lengths of 992 bps and 119 bps, respectively, and an average sequence length of 477 bps. We designed 102 primer pairs between the up-and down-flanking regions of the motifs based on the primer design parameters computed using FastPCR software version 6.4.18 [20]. To test the optimal annealing temperatures, which were obtained using gradient temperature PCRs, template DNA was derived from three individuals of each of the Keteleeria taxa, which are listed in Table 1. We tested the polymorphisms for K. davidiana var. formosana using 25 individuals from each population. Finally, we selected 16 polymorphic loci from the 102 microsatellites based on the detection of unambiguous polymorphic amplicons using fixed annealing temperature PCR and a polyacrylamide gel genotyping protocol. The characteristics of the 16 polymorphic microsatellite loci are listed in Table 2. Of the 16 loci, 11 are complete microsatellite loci, consisting of 5 with a dinucleotide motif, 5 with a trinucleotide motif, and 1 with a tetranucleotide motif. Of the 5 remaining loci, 3 carried a compound motif and 2 carried an interrupted motif. The sequences of 16 loci reported in this paper are available from GenBank (accession numbers: HG518488-HG518503) ( Table 2).
Pinglin populations, respectively, revealing population differentiation between the two remaining populations of the species. The genetic variability, including the means of the observed and expected heterozygosity (Table 3), of K. davidiana var. formosana was high compared with that of other Pinaceae species such as Pinus koraiensis (0.38 and 0.57) [24] and P. massoniana (0.27 and 0.65) [25,26], but was similar to that of P. pinaster (0.65 and 0.83) [27]. Unfortunately, no data for other Keteleeria taxa are available for the comparison of genetic variability. However, the high observed and expected heterozygosity values implied that the species was historically larger in population size and was more common in low-elevation forests than it is at present [7].

Test of transferability
To test the transferability of these microsatellite loci, we tested these primers in six other cow-tail fir taxa: K. fortunei var. cyclolepis, K. pubescens, K. davidiana var. calcarea, K. davidiana var. chienpeii, K. fortune, and K. evelyniana (Table 1). Two to three samples of each taxon were used in the evaluation of cross-amplification. All loci were transferable to the six taxa of Keteleeria (Table 4), and the annealing temperatures are listed on Table 2. However, few alleles were recovered in the six other Keteleeria taxa sampled here were caused by a small sampling effect. Nonetheless, the systematic and phylogeographic patterns among the Keteleeria species remain unclear. The transferability of these microsatellite loci among different Keteleeria taxa indicates the usefulness of these molecular genetic markers for interspecific and intraspecific research on such topics as phylogeography, speciation and introgression.

Conclusions
For conservation purposes, 16 new polymorphic microsatellite loci were developed from K. davidiana var. formosana. The genetic variation indices evaluated using these 16 polymorphic microsatellite loci for the two remaining populations of this endangered species are potentially useful for future studies that will focus on identifying distinct evolutionary units within the populations for conservation management. The interspecies transferability of these microsatellite loci may also be useful for future research aiming to reconstruct the phylogeographic patterns and the process of speciation among closely related species.

Availability of supporting data
The microsatellite sequences are available through the National Centre for Biotechnology Information (see http://www.ncbi.nlm.nih.gov/). The accession numbers on the repository are the following: GenBank accession number HG518488 through HG518503.