Development of microsatellite markers for partially and putative fully mycoheterotrophic varieties of Pyrola japonica sensu lato (Ericaceae)

We developed microsatellite markers to compare the genetic variation and reproductive biology between the partially mycoheterotrophic Pyrola japonica var. japonica and the putative fully mycoheterotrophic P. japonica var. subaphylla . Fif-teen primer pairs were developed for P. japonica sensu lato and they were tested on 77 ramets from three populations of the two varieties. Thirteen loci were polymorphic in at least one of the two var. japonica populations, whereas only four loci were polymorphic in the var. subaphylla population. The considerably lower genetic variation of the var. subaphylla population may be attributed to frequent selﬁng and/or inbreeding. The markers developed in this study will be useful for comparing the genetic diversity of P. japonica s. l. populations and measuring gene ﬂow within and between populations and varieties.


INTRODUCTION
Modification of the mating and breeding system is a major component of the evolution of plant mycoheterotrophy (Waterman et al., 2013), which is the ability to obtain carbon sources from fungi. Many fully mycoheterotrophic plants, which are nonphotosynthetic and mainly occur in shaded habitats, produce flowers that are specialized for autogamy (Waterman et al., 2013), possibly due to adaptation to their habitats, which typically lack effective pollinators. Tracking the changing mating system will provide insights into the processes that occur during the evolution of full mycoheterotrophy.
Pyrola japonica Klenze ex Alef. sensu lato, which includes the two varieties var. japonica and var. subaphylla (Maxim.) Andres, is widely distributed in East Asia (Takahashi, 1993;Liu et al., 2010). This species includes a total of three morphotypes in Japan; while var. japonica and var. subaphylla are mainly characterized by green scapes with well developed leaves (hereafter, the GL morphotype) and reddish scapes with rudimentary scale-like leaves (the RS morphotype), respectively, a third morphotype that has reddish or green scapes with intermediate leaf sizes and does not correspond to the two varieties is also found (the RL morphotype) (Shutoh et al., 2016). Consistent with a previous report that leaf size can be a rough predictor of the degree of mycoheterotrophy (Hynson et al., 2009), the GL and RS morphotypes exhibit different degrees of mycoheterotrophy: the GL morphotype exhibits partial mycoheterotrophy and the RS morphotype exhibits nearly full or full mycoheterotrophy (Takahashi, 1993;Matsuda et al., 2012;Shutoh et al., 2016). Because they have distinct chloroplast haplotypes, these three morphotypes can be genetically differentiated (Shutoh et al., 2016). Although most P. japonica s. l. populations are composed of a single morphotype, some populations are composed of two or more (Shutoh et al., 2016). Interestingly, the floral morphology is nearly identical in the species (Takahashi, 1986(Takahashi, , 1993, and there are no reports of pollinators. Therefore, little is known about reproductive isolation and any differences in reproductive biology among the three morphotypes. In the present study, we developed microsatellite markers for P. japonica s. l. to investigate the genetic diversity and genetic structure that may indicate differences in breeding systems among plant morphotypes with different degrees of mycoheterotrophy.

MATERIALS AND METHODS
Isolation of microsatellite loci We collected fresh leaf samples from one individual for each of the three P. japonica s. l. morphotypes (Table 1). Genomic DNA was extracted using the DNeasy Plant Mini Kit (Qiagen, Germantown, MD, USA). DNA fragment libraries were constructed using the Ion Xpress Plus Fragment Library Kit (Thermo Fisher Scientific, Waltham, MA, USA), amplified using the Ion PGM Template OT2 400 Kit (Thermo Fisher Scientific), and then sequenced using the Ion PGM Sequencing 400 Kit (Thermo Fisher Scientific) and an Ion 318 Chip v2 (Thermo Fisher Scientific). The total number of reads and the average read length obtained from the GL, RS, and RL morphotypes were 429,045 and 216 bp, 392,577 and 222 bp, and 481,396 and 211 bp, respectively. Subsequently, we used MSATCOMMANDER v. 0.8.2 (Faircloth, 2008) to identify microsatellite regions that contained more than eight di-or six trinucleotide repeats. As a result, 160, 83 and 50 potential microsatellite loci were identified in the RS, GL and RL morphotype, respectively. A total of 293 primer sets were designed using Primer3 v. 4.0.0 with default settings (Rozen and Skaletsky, 1999).
Preliminary amplification tests for all 293 primer pairs were conducted using two ramets from one GL population, two from one RS population, and four from two RL populations (Table 1). An M13 tag sequence (5′-CACGACGTTGTAAAACGAC-3′ or 5′-TGTGGAATT-GTGAGCGG-3′) was added to each forward primer for fluorescent labeling (Boutin-Ganache et al., 2001), and PCR amplifications were performed using the Qiagen Multiplex PCR Kit (Qiagen). Each reaction contained 10 ng of extracted DNA, 2.5 μl Multiplex PCR Master Mix, 0.01 μM fluorescently labeled forward primer, and 0.2 μM reverse primer in a final volume of 5 μl. The amplification process consisted of an initial denaturation at 95 °C for 15 min; 30 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 90 s and extension at 72 °C for 1 min; and a final extension at 60 °C for 30 min. Fragment sizes were determined using an ABI PRISM 3130 Genetic Analyzer and Gene Mapper software (Applied Biosystems, Foster City, CA, USA). Finally, we selected 15 primer pairs that yielded clear peak patterns in all eight samples tested (Table 2). Among the 15 primers, 12 were developed from RS, two from GL, and one from RL. Using the protocol described above, we also performed cross-amplification tests using four ramets each from P. incarnata (DC.) Freyn and P. nephrophylla Andres.
Characterization of microsatellite markers We used these 15 loci to examine the genetic variation of 52  (Goudet, 1995).

RESULTS AND DISCUSSION
In the two GL populations, 13 of the 15 loci were polymorphic, and all 52 ramets had distinctive multilocus genotypes. The range of H o and H e was 0.00-0.70 (mean = 0.33) and 0.00-0.78 (mean = 0.36), respectively (Table 3). Significant linkage disequilibrium (P < 0.05 GL, RS, and RL: loci developed from sequences of var. japonica, var. subaphylla, and intermediates between the two varieties, respectively; T a : annealing temperature; R A : allele size range; Tc: cross-amplification test; Pi: Pyrola incarnata; Pn: P. nephrophylla; M and P: monomorphic and polymorphic loci consistently amplified in four samples of each species, respectively. after Bonferroni correction) was not observed. Significant deviation from Hardy-Weinberg equilibrium was not observed, indicating a low frequency of null alleles. Meanwhile, in the RS population, only four loci were polymorphic, and the 25 RS ramets contained only six unique multilocus genotypes (Table 3). Significant deviation from Hardy-Weinberg equilibrium was observed for two loci, Pj-RS043 and Pj-RS097 (P < 0.05), and no significant linkage disequilibrium was observed (P < 0.05 after Bonferroni correction).
In the cross-amplification tests in P. incarnata, only one of the 15 loci was amplified, and it was polymorphic ( Table 2). On the other hand, in P. nephrophylla, six loci were amplified; one of the six was polymorphic and the others were monomorphic (Table 2). These loci were consistently amplified or unamplified in every sample of each species.
Our results demonstrate that P. japonica var. japonica and var. subaphylla exhibit different levels of genetic variation (Table 3). In the GL populations, relatively high genetic diversity was observed, suggesting that these populations are maintained by outcrossing. Conversely, in the RS population, low levels of genetic variability were found as shown in fewer polymorphic loci, and low H o and high F is values (Table 3). Although cross-amplification of microsatellite markers to related species generally indicates lower genetic variation and higher null allele frequency compared to that in the origi-nal species (Ellegren et al., 1997), the genetic diversity in the RS population was consistently low even at the 12 loci that were developed based on the RS sequences and successfully amplified in the GL populations. Therefore, the remarkably low genetic variation in the RS population is likely explained by selfing and/or inbreeding, not by null alleles. The putative difference in reproductive characteristics between the GL and RS populations is consistent with other plant species, in which the mating systems differ between full and partially mycoheterotrophic varieties (e.g., Bidartondo, 2005;Suetsugu, 2015;Suetsugu et al., 2015).
Overall, we developed 15 microsatellite markers for P. japonica s. l., which includes both partially and putative fully mycoheterotrophic varieties. At the 15 loci that were successfully amplified in both varieties, the var. subaphylla population exhibited lower genetic variation, which suggests that var. subaphylla switched to selfing and/or inbreeding for reproduction during the evolutionary process of acquiring full mycoheterotrophy. These 15 markers will be useful for comparing the genetic diversity of other P. japonica s. l. populations and examining gene flow within and between populations and varieties.
We are grateful to Masako Yamaguchi, Motoko Igari, Goshikinuma Kuruminomori, Fukushima-shi Kotorinomori, and Nichirei Corporation for facilitating fieldwork; Yumi Yamashita for help with the genetic analysis; and the associate editor and two anonymous reviewers for their constructive comments. This