Development and characterization of 14 microsatellite markers for Buergeria japonica (Amphibia, Anura, Rhacophoridae)

ABSTRACT

Buergeria japonica is a common frog species distributed throughout almost all islands in Ryukyu Archipelago. Because of their exceptionally wide distribution and higher physiological tolerance comparing to the other anurans, their demographic history and formation of distribution are intrinsic topics in the herpetological fauna of Ryukyu. Microsatellite marker is ideal genetic marker for such studies at inter- and intra-population level. We therefore developed microsatellite markers of B. japonica utilizing Ion PGM™ sequencing. As a result of the screening, we developed a total of 14 polymorphic markers. To test availabilities of these markers, we genotyped four island populations. The total number of alleles and expected hetelozygosities per locus ranged from 4 to 21 and 0.00 to 0.864, respectively. The phylogenetic relationship among the four populations based on the genetic distances of these markers was congruent with general divergence pattern of amphibians and reptiles in Ryukyu area. These markers developed in this study are considered to be useful for future studies about phylogeography and demography of this species.

MAIN BODY

The Ryukyu Kajika frog (Buergeria japonica) is a common frog species widely distributed in most of islands from South-Western Islands of Japan to Taiwan (Ryukyu Archipelago), a chain of numerous islands extending almost 650 km. Although frog species in the Ryukyu Archipelago generally reside only in one of three island areas (North-, Middle-, and Southern-Ryukyu) that are bordered by two deep sea straits (Tokara and Kerama straits), this species is exceptionally distributed in almost all the islands across these island areas (Maeda and Matsui, 1999). In addition to that, this species can breed along the coastal area and also has higher physiological tolerance against heat (Chen et al., 2001; Haramura, 2004; Wu and Kim, 2005). These facts may generate two hypotheses about demographic processes of this species: ancient expansion through land bridges or recent migration across strait by moving through short temporal land bridges or drifting on float lumbers. Thus the phylogenetic relationships and demographic processes of B. japonica are intrinsic topics to be clarified using the molecular markers for inter- and intra-population level. Microsatellite marker is an ideal neutral molecular marker for such fine scale evolutionary studies. However, the sequencing of microsatellite loci with adjacent region and development of the required PCR primers to amplify the microsatellite loci can be costly and time-consuming especially for the non-model organisms. In this respect, next-generation sequencing is the most promising tool to obtain genome sequence data and some studies have already succeeded in developing the microsatellite markers by using Roche 454 platforms (e.g. Schoebel et al., 2013; Rico et al., 2013; Gardner et al., 2011). At present, Life Technologies has developed new platform Ion Torrent and achieved faster and higher throughput in cheaper costs. Therefore, we attempted to isolate the microsatellite loci of B. japonica by using Ion PGM sequencer instead of the 454 platform or traditional screening and cloning method. In addition, we characterized the availability of these loci as molecular genetic markers for B. japonica populations.

During the breeding season of B. japonica, from February to September 2013, we collected total of 64 individuals from Taketomi in Iriomote Is., Kunigami in Okinawa Is., Fukumoto in Amami Is. and Taipei in Taiwan (16 individuals from each population) (Fig. 1). For each population, we caught adult frogs by bare-hands around a single breeding site (within approximately 250 m radius) and then clipped their toes for the genomic DNA extraction.

Fig. 1.

Sampling localities of B. japonica in this study.

We used Ion PGM™ sequencing to obtain sufficient genomic sequence data for isolation of microsatellite loci and design PCR primers to amplify the loci. Genomic DNA for sequencing was extracted from liver tissue of an individual from Taipei, Taiwan which have been stored in Hiroshima University by using NucleoSpin^® Tissue (Macherey & Nagel) and re-suspended in RNase A with the concentration of 20 μg/ml. Genomic DNA library was then constructed using NEBNext Fast DNA Fragmentation & Library Prep Set for Ion Torrent (New England Biolabs). Because we conducted a single sequencing run multiplexing with two individuals of the other species (Hoplobatrachus tigerinus and Odorrana narina), the genomic DNA library was attached by adapter oligo nucleotides using Ion Express Barcode Adaptors 1-16 Kit (Life Technologies). This library was amplified and refined using Ion OneTouch Template OT2 200 Kit (Life Technologies) and sequenced on the PGM™ with Ion 318 chip and Ion PGM Sequencing 200 Kit ver2 (Life Technologies).

As the result of the sequencing, we obtained 439,311 reads more than 150 bp in length. Of these raw sequence reads, we found 65 of di-, 22 of tri-, 9 of tetra-, and 4 of hexa-nucleotide repeats by using default setting of MSATCOMMANDER ver 1.0.8 (Faircloth, 2008). We selected total 70 microsatellite loci in order from the loci having larger number of repeats in each repeat type (39 of di-, 21 of tri-, 9 of tetra-, and 1 of hexa-nucleotide repeats) and tested amplification and polymorphism. Primer pairs for each locus were designed by Primer3 ver. 2.2.3 (Rozen and Skaletsky, 2000). All forward primers were attached with M13-tail (Schuelke, 2000) or BStag (Shimizu and Yano, 2011) at their 5’ end. We tested PCR amplification and polymorphism of all of the loci using 16 samples from Taipei population. The PCR for amplification confirmation was carried out in a 10 μl volume containing 5 μl of EmeraldAmp^® MAX PCR Master Mix (TaKaRa), 1 μl of 10 μM primer pairs, and 50 ng of genomic DNA. Genomic DNA for genotyping was extracted from toe tip of individuals in the same method mentioned above. Thermal cycling was performed under the following conditions: 95℃ for 5 min, 35 cycles of 95℃ for 15 sec, 58℃ for 30 sec and 72℃ for 30 sec. The PCR products were visualized by electrophoresis on a 2% agarose gel. Of these 70 loci, 37 were amplified judging by the presence of clear bands on an agarose gel in most of samples and then genotyped by post-labeling PCR and genotyping on 3130xl. M13 and BStag universal primers were labeled with HEX, FAM, NED, or PET fluorescent dyes. The PCRs were carried out in a 10 μl volume containing 5 μl of EmeraldAmp^® MAX PCR Master Mix, 0.1 μl of 1 μM forward and labeled universal primers, 0.2 μl of 2 μM reverse primer and 50 ng of genomic DNA. Thermal cycling was performed under the following conditions: 95℃ for 5 min, 35 cycles of 95℃ for 15 sec, 58℃ for 30 sec and 72℃ for 30 sec, followed by 8 cycles of 95℃ for 15 sec of 53 or 49℃ (M13 or BStag) for 30 sec and 72℃ for 30 sec, and a final extension period of 72℃ for 10 min. The PCR products were electrophoresed on 3130xl Genetic Analyzer (Life Technologies) together with GeneScan LIZ 500 (Life Technologies) as an internal size standard and genotyped using GeneMapper 4.0 (Life Technologies). Finally we screened 14 stably genotyped loci (Tables 1 and 2).

Table 1. Primer informations of 14 microsatellite loci developed in this study

Loci	Accession no.	Repeat motif	Primer sequence (5'-3')	Size range (bp)
Buerj657	AB902925	(AG)16	F: M13 - GCCCTCCACTTTCCTAAATCC	108–132
			R: CCTGTGAGCGGTGTACTAGG
Buerj584	AB902926	(AC)15	F: M13 - ACCTGCTCCTGGACAGAATG	151–169
			R: AAACAATGGGCCAGCAACAG
Buerj852	AB902927	(AG)13	F: M13 - GCAGTTTCATGCATTACAGAGC	157–208
			R: TCCTACTGCTCTTCCCAAGC
Buerj174	AB902928	(AC)11	F: M13 - TCTTGTGGACCCGTGTATGG	128–153
			R: ATAGCAAGGACGGTTCACAG
Buerj343	AB902936	(AAT)8	F: M13 - TGAGGGTAGGGACAGATGTG	108–132
			R: AAGAGCCAGTCAGACTTCCAC
Buerj83	AB902935	(AAT)9	F: M13 - CCCTTTGGAATGACACTGTTG	151–169
			R: AAGACTGCGATGTTAGGGAGAC
Buerj1101	AB902938	(AATT)7	F: M13 - CCAGCATGTGTGTAGGGATC	157–208
			R: CCTGTGCTTGTATCCTGCTAG
Buerj2469	AB902934	(ACT)10	F: M13 - ATGTTGGGTGTCTTCTGAACC	128–153
			R: CCCATTGGTCAGTCATGGTG
Buerj476	AB902937	(AAT)8	F: M13 - GACACAGTCTCGACAGATGG	126–149
			R: TGGTGCCCGTAATCTGGAAC
Buerj756	AB902931	(AC)9	F: F9TAC - AGAGGAATCAGGGAACAGGC	137–163
			R: AAAGCCACCGGAGTAGCC
Buerj342	AB902930	(AG)10	F: F9GCC - CATCAGCACACACAAAGTTCTG	151–161
			R: GGCAGTGTGATGTCAGCCC
Buerj134	AB902929	(AC)10	F: F9GTC - TGTGTAATTCTCAGCCGACG	125–163
			R: CCGTACCCTGTAGCTGAGTC
Buerj871	AB902932	(AC)9	F: F9GCC -CGCTCACTCCATTAATCCGC	102–141
			R: GTAAGTAAAGCGGCGGAGTG
Buerj1418	AB902933	(AC)9	F: F9GCC - TTTGCCTGTTATCTGCATTGC	180–196
			R: AAGGAGTGCAAGGTGAGCTG

Exponential values on the repeat motifs correspond to the numbers of repeats in the NGS data.

Table 2. Number of alleles and heterozygosities of 14 microsatellite loci in four island populations

Loci	All sites		Amami (N = 16)			Okinawa (N = 16)			Iriomote (N = 16)			Taiwan (N = 16)
	Na	FST	Na	HO	HE	Na	HO	HE	Na	HO	HE	Na	HO	HE
Buerj657	9	0.577	1	0.000	0.000	3	0.188	0.271	3	0.563	0.615	4	0.688	0.713
Buerj584	9	0.376	5	0.688	0.533	4	0.688	0.646	5	0.313	0.375	4	0.375	0.637
Buerj852	10	0.332	5	0.750	0.609	6	0.563	0.500	1	0.000	0.000	4	0.563	0.463
Buerj174	14					5	0.375	0.664	6	0.667	0.771	4	0.375	0.451
Buerj343	6	0.578	2	0.125	0.219	2	0.063	0.061	3	0.313	0.398	3	0.625	0.607
Buerj83	10	0.618	1	0.000	0.000	2	0.125	0.117	3	0.438	0.447	5	0.625	0.686
Buerj1101	7	0.509	2	0.438	0.498	3	0.500	0.529	2	0.250	0.219	2	0.313	0.482
Buerj2469	13	0.435	2	0.625	0.469	1	0.000	0.000	7	0.750	0.805	6	0.750	0.775
Buerj476	17	0.311	9	0.563	0.678	7	0.625	0.613	6	0.563	0.627	4	0.750	0.666
Buerj756	8	0.599	3	0.063	0.600	1	0.000	0.000	2	0.250	0.469	3	0.438	0.361
Buerj342	9	0.193	6	0.750	0.730	5	0.933	0.771	4	0.625	0.643	4	0.400	0.611
Buerj134	4								3	0.375	0.529	1	0.000	0.000
Buerj871	12	0.488	4	0.375	0.326	4	0.400	0.433	4	0.375	0.557	5	0.267	0.613
Buerj1418	21	0.295	7	0.750	0.773	11	0.600	0.864	5	0.500	0.492	4	0.125	0.578
mean	10.643	0.443	3.917	0.427	0.453	4.154	0.389	0.421	3.857	0.427	0.496	3.786	0.449	0.546

N_a: number of alleles, H_O and H_E: observed and expected heterozygosity (values significantly deviated from HWE are in bold number).

Using these 14 loci, we also genotyped the other three populations (Taketomi in Iriomote Is., Kunigami in Okinawa Is. and Fukumoto in Amami Is.) under the same procedure mentioned above. In populations of Okinawa and Amami islands, 1 and 2 loci were not successfully amplified, respectively (Table 2). In Iriomote population, all 14 loci were successfully genotyped including a monomorphic loci. The number of alleles (N_a) of each locus ranged from 4 to 21 in total populations. The observed and expected heterozygosities (H_O and H_E) calculated by GENALEX 6.5 (Peakall and Smouse, 2006) ranged from 0.000 to 0.933 and from 0.000 to 0.864, respectively. The tests for deficiency of Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were performed by Arlequin ver. 3.5 (Excoffier and Lischer, 2010). After Bonferroni correction for multiple testing, no significant LD was observed but the significant deviations from HWE were observed in Buerj476 and Buerj756 in Amami population and Buerj871 and Buerj1418 in Taiwan population. Because none of these loci showed deviation of HWE in the other populations, we didn’t exclude those from the analyses.

The pairwise distance (unbiased F_ST and D_A) matrix among populations calculated by POPTREE2 (Takezaki et al., 2010) is shown in Table 3 and from which neighbor joining tree was constructed (Fig. 2). The bootstrap value was calculated with 1000 replicates. The root of the tree was calculated by the mid-point rooting method, in which the root is placed in the mid-point of the longest path of two taxa (Takezaki et al., 2010). This tree showed closer relationship between Okinawa and Amami populations, being congruent with the general phylogenetic pattern of amphibians and reptiles in Ryukyu Archipelago, i.e., “Kerama gap” (Nishioka et al., 1987; Matsui, 1994; Ota, 1998). Moreover, F_ST and D_A values between Amami and Okinawa populations of B. japonica (0.486 and 0.437 of F_ST and D_A, respectively) were generally lower than those of other frogs endemic to the two islands such as Ishikawa’s frogs (Odorrana ishikawae and O. splendida; mean F_ST and mean D_A were 0.517 and 0.934, respectively, based on 12 SSR markers: Igawa et al., 2011, 2013) and Ryukyu tip-nosed frogs (O. narina and O. amamiensis: mean F_ST and mean D_A were 0.274 and 0.784, respectively, based on 17 SSR markers; Igawa et al., in preparation). Assuming equal mutation rates of these different markers, the lower genetic distances might reflect more recent divergence or higher migration rates between Amami and Okinawa populations in B. japonica than other frogs. In addition, our results revealed that N_a, H_O and H_E in all populations were in the same level (no significant differences were observed by ANOVA), despite the size of each island is extremely different. This suggests these populations have sustained sufficient population size even in small islands. These lower genetic divergences among populations and the same genetic variation within populations of B. japonica are also reasonable when considering their higher physiological tolerance and wide range of potential niche.

Table 3. Pairwise distances (unbiased F_ST and D_A) between each population

	Taipei	Iriomote	Okinawa	Amami
Taipei	–	0.430	0.463	0.425
Iriomote	0.825	–	0.530	0.477
Okinawa	0.821	0.825	–	0.437
Amami	0.790	0.773	0.486	–

Below diagonal: unbiased F_ST, above diagonal: D_A.

Fig. 2.

Phylogenetic relationships of the four island populations of B. japonica inferred from pairwise distance (D_A).

In conclusion, we developed sufficient number of microsatellite markers by utilizing the Ion Torrent PGM. Most of loci were stably amplified and showed polymorphism in all island populations used in this study. The phylogenetic tree based on genetic distances also agreed with divergence pattern of the other amphibian and reptile species in this area. Therefore, these markers are considered to be useful for clarifying dynamic evolutionary and demographic history of B. japonica in association with ecological and physiological feature of this species and geographic history of Ryukyu Archipelago.

ACKNOWLEDGMENTS

We would like to thank Lee Yu for his field assistance in Taipei. We also thank Chie Iwamoto, Nana Fukuda, Mai Fujimi, and Kanako Onizuka for their supports in the experiments. This work was supported by a Grant-in-Aid for JSPS Fellows (No. 25-5065).

REFERENCES

•  Chen,  T. C.,  Kam,  Y. C., and  Lin,  Y. S. (2001) Thermal physiology and reproductive phenology of Buergeria japonica (Rhacophoridae) breeding in a stream and a geothermal hotspring in Taiwan. Zool. Sci. 18, 591–596.
•  Excoffier,  L., and  Lischer,  H. E. L. (2010) Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567.
•  Faircloth,  B. C. (2008) MSATCOMMANDER: detection of microsatellite repeat arrays and automated, locus-specific primer design. Mol. Ecol. Resour. 8, 92–94.
•  Gardner,  M. G.,  Fitch,  A. J.,  Bertozzi,  T., and  Lowe,  A. J. (2011) Rise of the machines – recommendations for ecologists when using next generation sequencing for microsatellite development. Mol. Ecol. Resour. 11, 1093–1101.
•  Haramura,  T. (2004) Salinity and other abiotic characteristics of oviposition sites of the rhacophorid frog, Buergria japonica, in coastal habitat. Curr. Herpetol. 23, 81–84.
•  Igawa,  T.,  Okuda,  M.,  Oumi,  S.,  Katsuren,  S.,  Kurabayashi,  A.,  Umino,  T., and  Sumida,  M. (2011) Isolation and characterization of twelve microsatellite loci of endangered Ishikawa’s frog (Odorrana ishikawae). Conserv. Genet. Resour. 3, 421–424.
•  Igawa,  T.,  Oumi,  S.,  Katsuren,  S., and  Sumida,  M. (2013) Population structure and landscape genetics of two endangered frog species of genus Odorrana: different scenarios on two islands. Heredity 110, 46–56.
•  Maeda,  N., and  Matsui,  M. (1999) Frogs and Toads of Japan (Revised edition). Bun-Ichi Sogo Shuppan, Tokyo (in Japanese).
•  Matsui,  M. (1994) A taxonomic study of the Rana narina complex, with description of three new species (Amphibia: Ranidae). Zool. J. Linn. Soc. 111, 385–415.
•  Nishioka,  M.,  Ueda,  H., and  Sumida,  M. (1987) Intraspecific differentiation of Rana narina elucidated by crossing experiments and electrophoretic analyses of enzymes and blood proteins. Sci. Rep. Lab. Amphib. Biol. Hiroshima Univ. 9, 261–303.
•  Ota,  H. (1998) Geographic patterns of endemism and speciation in amphibians and reptiles of the Ryukyu archipelago, Japan, with special reference to their Paleogeographical implications. Res. Popul. Ecol. 40, 189–204.
•  Peakall,  R., and  Smouse,  P. E. (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol. Ecol. Notes 6, 288–295.
•  Rico,  C.,  Normandeau,  E.,  Dion-Côté,  A. M.,  Rico,  M. I.,  Côté,  G., and  Bernatchez,  L. (2013) Combining next-generation sequencing and online databases for microsatellite development in non-model organisms. Sci. Rep. 3, 3376.
•  Rozen,  S., and  Skaletsky,  H. J. (2000) Primer3 on the WWW for general users and for biologist programmers. In: Bioinformatics Methods and Protocols: Methods in Molecular Biology (eds: S. Misener and S. A. Krawetz), pp. 365–386. Humana Press, Totowa, New Jersey.
•  Schoebel,  C. N.,  Brodbeck,  S.,  Buehler,  D.,  Cornejo,  C.,  Gajurel,  J.,  Hartikainen,  H.,  Keller,  D.,  Leys,  M.,  Ríčanová,  S.,  Segelbacher,  G.,  Werth,  S., and  Csencsics,  D. (2013) Lessons learned from microsatellite development for nonmodel organisms using 454 pyrosequencing. J. Evol. Biol. 26, 600–611.
•  Schuelke,  M. (2000) An economic method for the fluorescent labeling of PCR fragments. Nat. Biotechnol. 18, 233–234.
•  Shimizu,  T., and  Yano,  K. (2011) A post-labeling method for multiplexed and multicolored genotyping analysis of SSR, indel and SNP markers in single tube with bar-coded split tag (BStag). BMC Res. Notes 4, 161.
•  Takezaki,  N.,  Nei,  M., and  Tamura,  K. (2010) POPTREE2: Software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. Mol. Biol. Evol. 27, 747–752.
•  Wu,  C. S., and  Kam,  Y. C. (2005) Thermal tolerance and thermoregulation by Taiwanese Rhacophorid tadpoles (Buergeria japonica) living in geothermal hot springs and streams. Herpetologica 61, 35–46.