Isolation and characterization of microsatellite loci for Smilax sieboldii (Smilacaceae).

PREMISE OF THE STUDY: Polymorphic microsatellite markers were developed for Smilax sieboldii (Smilacaceae), a member of the S. hispida group with a biogeographic disjunction between eastern Asia and North America, to study the phylogeography and incipient speciation of this species and its close relatives. METHODS AND
RESULTS: Transcriptome sequencing produced 47,628 unigenes. Seventeen loci were developed from 122 randomly selected primer pairs. Polymorphism and genetic variation were evaluated for 68 accessions representing five populations of S. sieboldii. The number of alleles per locus ranged from four to 18; the expected heterozygosity ranged from 0.59 to 0.92. Twelve loci were successfully amplified in five related species: S. scobinicaulis, S. californica, S. hispida, S. moranensis, and S. jalapensis.
CONCLUSIONS: These novel expressed sequence tag-derived microsatellite markers will facilitate further population genetic research of S. sieboldii and its close allies of the S. hispida group.

The Smilax hispida group is a well-supported clade including six species in Smilacaceae (Qi et al., 2013) with a disjunct distribution including eastern Asia (S. sieboldii Miq. and S. scobinicaulis C. H. Wright), western North America (S. californica (A. DC.) A. Gray), eastern North America (S. hispida Raf.), and Mexico (S. moranensis M. Martens & Galeotti and S. jalapensis Schltdl.). Smilax sieboldii is a typical element of temperate broad-leaved forests that occurs widely in mainland China, Taiwan, Japan, and Korea. Previous studies based on two cpDNA intergenic regions indicated that at least four biogeographic lineages exist, with each lineage containing at least one private haplotype. This phylogeographic structure is considered to be related to the historical fluctuation of climate and sea level (Zhao et al., 2013). However, this study was limited by the lack of nuclear markers. Therefore, polymorphic microsatellite markers will enhance our understanding of population genetic diversity and historical demography (e.g., gene flow, genetic bottlenecks) and will allow for connecting these patterns to geological and environmental changes.

Existing microsatellite markers for Smilax species (Xu et al., 2011; Martins et al., 2013) showed limited transferability and polymorphism for the S. hispida group due to phylogenetic distance. Therefore, in the current study we aimed to develop more polymorphic and transferable expressed sequence tag–simple sequence repeat (EST-SSR) markers from the transcriptome, which contains abundant ESTs, based on a high-throughput sequencing approach.

METHODS AND RESULTS

Transcriptome sequencing

Fresh young leaves of one wild accession of S. sieboldii were collected at Tianmu Mountain, Zhejiang Province, China (Appendix 1), and frozen in liquid nitrogen. RNA was extracted using TRIzol Reagent (Invitrogen Life Technologies, Carlsbad, California, USA) and treated with DNase (TaKaRa Bio, Shuzo, Kyoto, Japan) following the manufacturer’s instructions. A 2 × 150-bp paired-end RNA-Seq library was prepared following the normalized eukaryote transcriptome library preparation protocol of the Beijing Genomics Institute (Shenzhen, China) and sequenced on the Illumina HiSeq 2500 platform (Illumina, San Diego, California, USA). A total of 65,863,062 raw reads were generated and uploaded to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (accession SRP095761). The raw data were filtered using FASTX-TOOLKIT version 0.0.14 (Gordon and Hannon, 2010) by removing adapter sequences and low-quality reads with >5% unknown bases and/or >15% low-quality bases (quality value <20). Remaining reads were assembled into 66,482 transcripts using TRINITY version 2.3.2 (Grabherr et al., 2011), which were then clustered into 47,628 unigenes with TGICL version 2.1 (Pertea et al., 2003).

Microsatellite development

Using the MIcroSAtellite identification tool (MISA) (Thiel et al., 2003), microsatellite regions in the unigenes were screened according to the following criteria for repeat numbers: dinucleotide repeats ≥6, trinucleotide repeats ≥5, and tetranucleotide, pentanucleotide, and hexanucleotide repeats ≥4. Primers were designed for the screened microsatellite loci using Primer3 (Untergasser et al., 2012) with the default parameter settings. A total of 9263 microsatellite sequences were obtained, from which 2252 primer pairs were designed. Of these, 122 primer pairs were randomly selected and their forward primers were synthesized with one of three different universal primers (5′-CACGACGTTGTAAAACGAC-3′, 5′-TGTGGAATTGTGAGCGG-3′, or 5′-CTATAGGGCACGCGTGGT-3′) (Boutin-Ganache et al., 2001; Sakaguchi and Ito, 2014). To prevent primer dimers, hairpin structures, and mismatches, the best matches of forward primers and universal primers were selected using OLIGO version 6.67 (Molecular Biology Insights, Cascade, Colorado, USA).

We selected 12 accessions from various populations (Appendix 1) to test the effectiveness of primer amplification and to preliminarily assess genetic variation. Total genomic DNAs were extracted from silica-dried leaves using Plant DNAzol (Invitrogen Life Technologies). PCR amplifications were performed following the standard protocol of the Tsingke PCR kit (Tsingke Biotech Company, Beijing, China) in a final volume of 10 μL, which contained approximately 5 ng of DNA, 5 μL of 2× PCR Master Mix, 0.1 μM of forward primer, 0.4 μM of reverse primer, and 0.3 μM of fluorescently labeled universal primer (FAM, ROX, HEX, TAMRA; Table 1). The PCR thermal profile involved an initial denaturation at 95°C for 5 min; followed by 35 cycles of 94°C for 40 s, 58°C for 30 min, 72°C for 30 s; and a final 10-min extension step at 72°C. Fragment lengths of PCR products were analyzed on a 3730xl DNA Analyzer (Applied Biosystems, Foster City, California, USA) with GeneScan 500 LIZ as an internal reference (Applied Biosystems). Electrophoresis peaks were scored using GeneMarker version 2.2.0 (SoftGenetics, State College, Pennsylvania, USA). A total of 17 primer pairs with stable repeatability and high variation were selected for further analysis. All primer sequences obtained from this study were submitted to GenBank (Table 1).

Table 1.

Characteristics of 17 newly developed microsatellite loci in Smilax sieboldii.

Locus	Primer sequences (5′–3′)	Repeat motif	Allele size range (bp)b	Fluorescent dyec	GenBank accession no.	Functiond	Organism	E-value
SS2	F: ACTGTAGGAGTTGAGCACAGAGG	(GA)17	60–100	FAM	KY404961	Auxin response factor 15	Oryza sativa subsp. japonica	0
	R: AGATTCGGGAAAACAGAGGAAT
SS5	F: CAACCCAAAACAAAACAAGAGAG	(AG)12	96–132	TAMRA	KY404962	Hydrolase protein 30	Arabidopsis thaliana	5E-24
	R: GATACACGGGTAACCACCACC
SS19	F: ACTTTGCCTATTAAGCATCCGTT	(CT)10	116–154	ROX	KY404963	Polygalacturonase inhibitor	Pyrus communis	5E-98
	R: AGTACTGCTTCCTCCACAACAAG
SS20	F: AACACACGATCTCAAAGAAGAGC	(GAA)15	89–122	FAM	KY404964	Protein FAF-like, chloroplastic	Arabidopsis thaliana	6E-18
	R: CGTCGTCATCTTCTTCTCTGTTT
SS21	F: GAATCCTTTCGCTTAGGGAAGT	(CT)12	107–137	TAMRA	KY404965	Probable ADP-ribosylation factor GTPase-activating protein AGD14	Arabidopsis thaliana	2E-15
	R: CACAAAGAATAAAAGAACGCTCG
SS33	F: AGTAGGATCCCAGCTTTCTTGAG	(AG)11	141–179	HEX	KY404966	Uncharacterized protein At4g08330, chloroplastic	Arabidopsis thaliana	2E-32
	R: CTCTCTCATCCCCAAATGTTTCT
SS43	F: CAAGTATCCACAACGAAAACCAT	(GA)11	154–180	HEX	KY404967	Oxygen-evolving enhancer protein 2, chloroplastic	Fritillaria agrestis	7E-119
	R: GTGGAGGAAACATGCAGTTGAT
SS74	F: GACGGCACCAAGAGAAGAAT	(CTG)8	181–241	FAM	KY404968	—	—	—
	R: GTGGATATCATCACCTCGGG
SS95	F: GTAGAGGCGCTGGGTTCC	(TGG)8	135–180	ROX	KY404969	Sulfated surface glycoprotein 185	Volvox carteri PE	3E-06
	R: GCCAAGCTCTGGAAGAACAC
SS100	F: GATTAGTGAGAGCTTGGCGG	(GAG)9	137–170	TAMRA	KY404970	Threonine-protein kinase-like protein At5g23170	Arabidopsis thaliana	2E-64
	R: ATGCACCAACTCCTTCCAAC
SS103	F: ACCATCTGTCCCAGTTGCAT	(TGG)10	263–281	ROX	KY404971	E3 ubiquitin-protein ligase At1g12760	Arabidopsis thaliana	8E-24
	R: CTCCCGAGGTTGTCAAAGAG
SS108	F: AAAGGCCCCCAATTATCATC	(TGC)13	106–124	FAM	KY404972	Formin-like protein 5	Oryza sativa subsp. japonica	5E-29
	R: CGGCTGGAGAAGATGAACTC
SS109	F: CCGGCAAGTATTGAGGATGT	(ATC)14	139–175	HEX	KY404973	—	—	—
	R: GGTGGAAGAGCTCAAAGACG
SS113	F: CTGATTTCCTTCCTGTTACGTTG	(CTGT)6	132–172	TAMRA	KY404974	—	—	—
	R: CAAATAACCGACTTCAGCTCCTA
SS114	F: TATTCGTGTAAAGATACGTGGGC	(GTGTGA)9	137–167	ROX	KY404975	DNA-directed RNA polymerase II subunit 1	Arabidopsis thaliana	6E-09
	R: TCGGCCATTATTTTAATCACATC
SS120	F: ATATGCCGTCGAGTATCGTCTT	(GCAGTA)4	146–200	ROX	KY404976	ABC transporter G family member 14	Arabidopsis thaliana	0
	R: GAGGAGGTGGTGTACAGGGTAAG
SS122	F: GACGGACTGACTGATACTTGGAT	(TAGCAC)4	125–185	HEX	KY404977	Protein PHLOEM PROTEIN 2-LIKE A1	Arabidopsis thaliana	7E-13
	R: GGAATACTCAAGTTCGCCGTATC

An annealing temperature of 58°C was used for all loci.

Size range values based on 68 individuals.

Forward 5′ label.

The unigenes containing microsatellite loci were searched against the SWISS-PROT database (http://www.expasy.ch/sprot/); — = not found.

Polymorphism assessment

To further evaluate the applicability of these primers, 68 individuals from five representative populations from China, Korea, and Japan (Appendix 1) were used to calculate genetic variation parameters. DNA extraction, PCR amplification, and length assessment of PCR products were performed following the procedures described above. The presence of null alleles and their bias on genetic diversity were evaluated based on the expectation maximization method implemented in FreeNA (Chapuis and Estoup, 2007). Deviation from Hardy–Weinberg equilibrium for each population and linkage disequilibrium for each primer pair were tested using GENEPOP version 4.0.7 (Rousset, 2008). The number of alleles, observed heterozygosity, expected heterozygosity, and polymorphism information content were calculated to assess the genetic polymorphism at each locus using CERVUS version 3.0.3 (Kalinowski et al., 2007).

Two loci (SS20, SS95) with high occurrence of null alleles (>5%) were excluded from the following analysis. No significant deviation from Hardy–Weinberg equilibrium (P < 0.001) was observed for the remaining 15 loci except SS5 in populations CZJ and JFS; SS19 in population KMJ; and SS21, SS100, and SS109 in population JFS, which might be caused by Wahlund effect of specific populations. There was no evidence of significant linkage disequilibrium in any pair of loci. We detected 156 alleles in total, and the number of alleles at each locus ranged from four to 18, suggesting a moderate to high level of polymorphism. The observed heterozygosity, expected heterozygosity, and polymorphism information content for each locus ranged from 0.36 to 0.97, 0.59 to 0.92, and 0.53 to 0.91, respectively (Table 2).

Table 2.

Genetic properties of the 15 newly developed microsatellite loci for Smilax sieboldii. Loci SS20 and SS95 are not included due to a high proportion (>5%) of null alleles.

	CTW (n = 6)				CZJ (n = 14)				CJS (n = 15)				KMJ (n = 16)				JFS (n = 17)				Total (n = 68)
Locus	A	Ho	He	PIC	A	Ho	He	PIC	A	Ho	He	PIC	A	Ho	He	PIC	A	Ho	He	PIC	A	Ho	He	PIC
SS2	4	0.83	0.77	0.65	8	1.00	0.86	0.80	8	0.87	0.85	0.80	9	0.88	0.86	0.82	10	0.82	0.84	0.79	17	0.88	0.92	0.90
SS5	4	1.00	0.78	0.65	10	0.86	0.88	0.83*	7	0.93	0.83	0.78	8	1.00	0.81	0.76	6	0.94	0.82	0.76*	14	0.94	0.90	0.88
SS19	4	0.83	0.74	0.62	7	0.79	0.74	0.67	4	0.60	0.71	0.63	8	0.81	0.87	0.83*	7	0.65	0.79	0.74	13	0.72	0.87	0.85
SS21	2	0.67	0.49	0.35	5	0.71	0.73	0.65	6	0.67	0.81	0.75	1	0.00	0.00	0.00	2	0.00	0.51	0.37*	9	0.36	0.79	0.76
SS33	5	1.00	0.74	0.64	10	0.86	0.89	0.84	4	1.00	0.76	0.68	7	1.00	0.77	0.71	10	1.00	0.88	0.83	18	0.97	0.92	0.91
SS43	5	1.00	0.82	0.70	5	0.71	0.77	0.71	4	1.00	0.72	0.64	5	1.00	0.78	0.71	8	0.71	0.70	0.66	12	0.87	0.87	0.85
SS74	8	1.00	0.91	0.81	5	0.46	0.63	0.55	5	0.40	0.36	0.34	4	0.56	0.60	0.50	9	0.77	0.79	0.74	16	0.60	0.77	0.75
SS100	4	0.67	0.71	0.60	7	0.86	0.77	0.71	6	0.73	0.73	0.67	5	0.86	0.74	0.67	6	0.38	0.76	0.70*	10	0.69	0.82	0.79
SS103	3	0.60	0.69	0.55	4	0.46	0.64	0.54	4	0.73	0.72	0.63	5	0.94	0.80	0.73	4	0.50	0.56	0.48	5	0.66	0.76	0.71
SS108	4	0.83	0.76	0.64	5	1.00	0.68	0.59	2	0.87	0.51	0.37	4	0.94	0.60	0.50	4	0.94	0.64	0.54	7	0.93	0.63	0.55
SS109	4	0.50	0.56	0.48	5	0.71	0.75	0.68	3	0.67	0.67	0.58	5	0.69	0.80	0.74	4	0.47	0.64	0.56*	9	0.62	0.84	0.82
SS113	4	1.00	0.76	0.64	4	0.71	0.55	0.45	5	0.87	0.63	0.56	4	0.44	0.47	0.43	3	0.53	0.42	0.34	9	0.66	0.59	0.53
SS114	3	0.67	0.55	0.45	3	0.79	0.62	0.53	4	0.87	0.66	0.57	3	0.31	0.46	0.40	3	0.77	0.55	0.47	4	0.68	0.69	0.62
SS120	3	1.00	0.67	0.54	6	0.50	0.72	0.65	2	0.40	0.41	0.32	4	0.69	0.73	0.65*	2	0.31	0.27	0.23	8	0.52	0.67	0.64
SS122	4	0.83	0.77	0.65	3	0.33	0.45	0.37	4	0.33	0.41	0.37	4	0.86	0.66	0.57	4	0.67	0.65	0.57	5	0.58	0.70	0.64

Note: A = number of alleles sampled; He = expected heterozygosity; Ho = observed heterozygosity; n = number of individuals sampled; PIC = polymorphism information content.

Voucher and locality information are provided in Appendix 1.

Significant deviation from Hardy–Weinberg equilibrium (P < 0.001).

Transferability evaluation

Transferability of the 15 primers was examined in the accessions of the five related species, i.e., five accessions each for S. californica, S. hispida, S. moranensis, and S. jalapensis and 10 accessions for S. scobinicaulis (Appendix 1). All loci were successfully amplified except two loci (SS21 and SS100) for S. hispida and one (SS33) for S. moranensis (Table 3). Polymorphism was detected in all but two loci (SS21 and SS100) for S. californica, five (SS2, SS19, SS103, SS120, and SS122) for S. hispida, four (SS21, SS74, SS103, and SS114) for S. moranensis, and one (SS100) for S. jalapensis (Table 3). The levels of both cross-amplifiability and polymorphism largely decreased with increasing phylogenetic distance. In total, 12 loci were amplifiable across the other five species in the S. hispida group.

Table 3.

Fragment sizes detected in cross-amplification tests of the 15 newly developed microsatellite markers in the remaining five species of the Smilax hispida group.

Locus	S. scobinicaulis (n = 10)	S. californica (n = 5)	S. hispida (n = 5)	S. moranensis (n = 5)	S. jalapensis (n = 5)
SS2	66–80	72–84	72	72–76	66–84
SS5	98–124	114–116	114–116	114–124	114–116
SS19	116–150	132–136	140	132–138	128–138
SS21	125–127	125	—	131	123–129
SS33	157–179	167–177	167–179	—	167–179
SS43	164–176	168–188	168–170	164–176	166–168
SS74	184–241	193–199	202–217	196	199–214
SS100	152–170	164	—	152–164	164
SS103	263–278	272–278	272	278	257–278
SS108	106–118	109–118	106–118	106–109	94–118
SS109	172–175	127–142	136–142	127–136	127–151
SS113	132–164	160–164	160–164	156–164	132–156
SS114	137–155	137–149	143–149	137	137–155
SS120	146–164	170–182	158	152–170	164–182
SS122	125–179	167–191	173	173–191	167–179

Note: — = amplification failed.

Voucher and locality information are provided in Appendix 1.

CONCLUSIONS

Using high-throughput sequencing, we sequenced and assembled the transcriptome of S. sieboldii without a reference genome. Fifteen EST-SSR markers were successfully developed to evaluate the genetic structure and demography of S. sieboldiii, of which 12 are likely to be useful for all six species of the S. hispida group.

Appendix 1.

Voucher information for Smilax species used in this study.

Species	Population code	Voucher specimensa	Collection locality	Geographic coordinates	n
Smilax sieboldii Miq.	CTW	Xiaoxian Liu, 0812003	Mt. Zhu, Taiwan, China	23.31000N, 120.50000E	6
Smilax sieboldii	CZJ	Yalu Ru, Ru150921001	Mt. Tianmu, Zhejiang, China	30.37809N, 119.42061E	14
Smilax sieboldii	CJS	Yunpeng Zhao, HZU00441	Mt. Longchi, Jiangsu, China	31.24818N, 119.74551E	15
Smilax sieboldii	KMJ	Joongku Lee, GG13	Myeongjisan, Gyeonggi-do, Korea	37.93458N, 127.47325E	16
Smilax sieboldii	JFS	Chengxin Fu & Xinjie Jin, Fu1505092	Fujiyama, Tokyo, Japan	35.50281N, 138.76985E	17
Smilax scobinicaulis C. H. Wright		Pan Li, LP150444	Mt. Wuzhi, Hubei, China	31.08961N, 110.88390E	10
Smilax californica (A. DC.) A. Gray		Pan Li, LP150436	Near Shasta Lake, CA, USA	40.75954N, 122.03657W	5
Smilax hispida Raf.		Yunpeng Zhao, 090834	Croatan National Forest, NC, USA	36.20339N, 86.98333W	5
Smilax jalapensis Schltdl.		Pan Li, US10041	Teopisca, Chiapas, Mexico	16.57310N, 92.50445W	5
Smilax moranensis M. Martens & Galeotti		Pan Li, US10031	Mexico City, Mexico	19.30541N, 99.30743W	5

Note: n = number of individuals sampled.

Vouchers were deposited in the Herbarium of Zhejiang University (HZU), Hangzhou, Zhejiang, China.