Journal of Threatened Taxa | www.threatenedtaxa.org | 26 December 2020 | 12(17): 17263–17275

ISSN 0974-7907 (Online) | ISSN 0974-7893 (Print)

doi: https://doi.org/10.11609/jott.6532.12.17.17263-17275

#6532 | Received 08 August 2020 | Finally accepted 08 December 2020

Genetic and reproductive characterization of distylous Primula reinii in the Hakone volcano, Japan: implications for conservation of the rare and endangered plant

Masaya Yamamoto ¹, Honami Sugawara ², Kazuhiro Fukushima ³, Hiroaki Setoguchi ⁴ & Kaoruko Kurata ⁵

¹ Hyogo University of Teacher Education, 942-1 Shimokume, Kato-city, Hyogo 673-1494, Japan.

^2,3,5 Graduate School of Education, Yokohama National University, 79-2 Tokiwadai, Hodogaya-ku, Yokohama, Kanagawa 240-8501, Japan.

⁴ Graduate School of Human and Environmental Studies, Kyoto University, Yoshida Nihonmatsu, Sakyo-ku, Kyoto 606-8501, Japan.

¹ myamamo@hyogo-u.ac.jp (corresponding author), ² xxhonami.s@gmail.com, ³ kazuhiro-kazuhiro-73@docomo.ne.jp, ⁴ setoguchi.hiroaki.2c@kyoto-u.ac.jp, ⁵ kaoruko@ynu.ac.jp

Editor: Mandar Paingankar, Government Science College Gadchiroli, Gadchiroli, India. Date of publication: 26 December 2020 (online & print)

Citation: Yamamoto, M., H. Sugawara, K. Fukushima, H. Setoguchi & K. Kurata (2020). Genetic and reproductive characterization of distylous Primula reinii in the Hakone volcano, Japan: implications for conservation of the rare and endangered plant. Journal of Threatened Taxa 12(17): 17263–17275. https://doi.org/10.11609/jott.6532.12.17.17263-17275

Copyright: © Yamamoto et al. 2020. Creative Commons Attribution 4.0 International License. JoTT allows unrestricted use, reproduction, and distribution of this article in any medium by providing adequate credit to the author(s) and the source of publication.

Funding: This work was financially supported by ProNatura Foundation Japan, Support Society of Faculty of Education and Human Science of Yokohama National University, and Grants-Aid for Science Research from Japan Society for Promotion of Science (#19K16219).

Competing interests: The authors declare no competing interests.

Author details: Masaya Yamamoto is an assistant professor in Hyogo University of Teacher Education. Honami Sugawara is a school teacher in Tokyo, Japan. Kazuhiro Fukushima is also a school teacher in Shizuoka Pref., Japan. Hiroaki Setoguchi is a professor in Kyoto University. Kaoruko Kurata is an associate professor in Yokohama National University.

Author contribution: KK conceived the project; MY, HS and KF collected the data; MY led the writing of the manuscript. HS and KK contributed critically to the drafts and give final approval for publication.

Acknowledgements: Fieldworks are conducted under permits from Kanto Regional Environment Office, Ministry of the Environment. We are grateful to Misato Wakai and Saki Sugihara, undergraduate and graduate students in the Kurata laboratory (Yokohama National University) for help in field surveys.

Abstract: Genetic and ecological evaluation are crucial in effective management of rare and endangered species, including those exhibiting complex breeding systems such as distyly. We studied a threatened distylous herb Primula reinii in the Hakone volcano, central Japan, to obtain baseline information of reproductive and genetic status towards conservation. In two representative populations inhabiting a central cone and somma of the volcano, population size, floral morph ratio, stigmatic pollen deposition, and fruit-set were measured. Using microsatellite markers, we evaluated genetic diversity, structure and differentiation of populations. Population bottlenecks and historical changes in population size were also estimated from genotype data. We found significant deviation from equal morph ratios in the central cone population, which also exhibited skewed mating success together with a high frequency of pollination within the same morph. These trends were not detected in the somma population. From genetic insights, the central cone population showed slightly lower genetic diversity, whereas no significant deviation from Hardy-Weinberg equilibrium was found in either population. The estimated moderate genetic differentiation and admixed genetic structure suggest recent lineage divergence and/or gene flow between populations. While robust evidence for a recent bottleneck was not obtained in our analyses, a clear signature of historical population contraction was detected in the central cone population. Our findings suggest that the skewed morph ratio strongly influenced the reproduction of small and isolated populations in the short-term, highlighting the vulnerability of distylous plant populations under ongoing anthropogenic pressure.

Keywords: Distyly, morph bias, reproduction, genetic diversity, volcanism.

INTRODUCTION

Worldwide, numerous plants are already threatened by human-caused stress (e.g., habitat destruction) and climate changes (Jackson & Kennedy 2009). Among these plants, species having a sophisticated entomophilous breeding system such as distyly (heterostyly) are likely to be the most vulnerable to the detrimental effects of isolation and unreliable pollination service due to anthropogenic environmental alteration (Washitani et al. 2005).

Distyly is a floral polymorphism, where populations have two floral morphs (a long- and short-styled morph; hereafter, referred to as the L- and S-morphs) that differ reciprocally in the heights of stigmas and anthers in flowers. Besides the morphological differences, distylous plants usually have a heteromorphic incompatibility system that prevents selfing and intra-morph mating (Barrett 2002); only cross-pollination (i.e., legitimate pollination) between L- and S-morphs results in seed setting. In theory, such morphologically and physiologically disassortative mating between floral morphs generally leads to an equilibrium with equal morph ratios, as a result of negative frequency-dependent selection and simple inheritance of distyly (Heuch 1979; Barret & Shore 2008). Accumulated evidence in distylous plants, however, has provided numerous examples of variation in population morph ratios (e.g., Kéry et al. 2003; Brys et al. 2008; Meeus et al. 2012). It has been advocated that floral morph bias can be governed by several factors, such as stochastic and deterministic events (Matsumura & Washitani 2000; Kery et al. 2003), maternal fitness differences between morphs (Hodgins & Barret 2006), and a combination of weak heteromorphic incompatibility and pollen limitation (Barret 1989; van Rossum et al. 2006; Brys et al. 2008).

Skewed morph ratios have often been found in small isolated distylous plant populations (e.g., Matsumura & Washitani 2000; Kery et al. 2003). Furthermore, it is well known that the deviation of morph frequencies from a 1:1 ratio can have negative reproductive and genetic consequences for populations. Indeed, skewed morph ratios result in the limited availability of compatible mates, which can contribute to reduced reproductive success (Kery et al. 2003; Wang et al. 2005; Pedersen et al. 2016) and increase the effects of genetic drift (Byers & Meagher 1992). Moreover, a combination of the loss of effective pollinators and the absence of a strict heteromorphic incompatibility system can increase inbreeding (selfing and biparental inbreeding) in morph-biased populations (Barret 1989; Wilcock & Neiland 2002; van Rossum & Triest 2006). Another consequence of skewed morph ratios is an increase in genetic drift and inbreeding, which can lead to a loss of genetic diversity (van Rossum & Triest 2006; Meeus et al. 2012) and may eventually result in the breakdown of distyly (Washitani 1996). Therefore, conservation studies on threatened distylous species which integrate ecological and genetic approaches are indispensable for assessing current status and predicting future extinction probability, as well as for planning effective conservation and/or restoration strategies (Washitani et al. 2005).

The present paper investigated Primula reinii Franch. et Sav. (Primulaceae), an endemic primrose that inhabits mountainous regions of Japan. As in most primroses, the plant is a self-incompatible distylous species (Richards 2003). Because of its attractive and relatively large flowers with dwarf foliage, the primrose has suffered from anthropogenic activity (e.g. horticultural exploitation) in the wild. Based on the rarity and serious reductions in numbers and populations, the species is listed in the ‘Vulnerable’ category of the latest Japanese Red List (Ministry of the Environment 2019). Despite required effective management, especially for small populations exposed to ongoing human activities, little is known about their population status, reproductive success, or remnant genetic diversity. Simultaneously, this might provide an opportunity to study the immediate responses of a distylous plant population to demographic changes. In this study, we assessed the genetic and reproductive status of two neighborhood populations of P. reinii to provide baseline information pertinent to the conservation and preservation of this rare and endangered primrose. To discuss factors affecting reproductive success in natural habitats, we measured population size, morph frequency, stigmatic pollen deposition, and fruit-set within a population, and evaluated the genetic diversity and structure within and among populations.

MATERIALS AND METHODS

Plant species and study site

Primula reinii is a diploid (2n = 24) perennial herb occurring as a chasmophyte on wet shaded rocky cliffs in the mountains (Image 1). The natural populations are rare and usually small and isolated from each other because of their narrow edaphic niche and low dispersal ability arising from the nature of the species. A single ramet produces one to two pink flowers from mid-April to early May. Although flower visitors are not well-known in this species, its narrow corolla tube and recent pollinator observations in their related species (Yamamoto et al. 2018) imply that long-tongued flying insects, such as bumblebees and bee flies can be effective pollinators for the species. Under cultivation conditions, the generation time (between seed germination and first flowering) of the species is estimated to be 2–3 years (Yamamoto et al. 2017).

Fieldwork was conducted from April 2013 to October 2014. We selected P. reinii populations from two sites in the Hakone volcano within a special protected zone of the Fuji-Hakone National Park in central Japan (Fig. 1a). These sites were severely isolated by a volcanic landform, i.e., caldera (distance: ca. 7 km) (Fig. 1b). One was on Mt. Kintoki (KIN population, 35.289’N, 139.004’E, 1,212m) and the other was on Mt. Komagatake (KOM population, 35.228’N, 139.021’E; 1,275m). In both sites, the primrose occurs on rock outcrops near the mountain top in which populations are composed of scattered patches within an area of approximately 20 × 20 m². Although a few smaller patches also located alongside the studied populations, these rim populations were sparse and separated geographically from the studied center population (>100 m). Thus, each studied population was considered to be one reproductive unit.

Formation and strong eruptive activities of the Hakone volcano initiated 650–350 ka and continued until 3ka (Nagai & Takahashi 2008). The mean annual precipitation and average air temperature at this area are 2,132mm and 8.8⁰C, respectively (http://en.climate-data.org/location/769594/, accessed 7 February 2019). Currently, the Hakone volcano is an attractive destination to tourists. There is more than 10,000 people living within the caldera, and an approximate average of 60,000 tourists visit the area every day. Including Primula reinii, approximately 80% of endangered plants found in the volcano are presumed to have decreased due to habitat destruction and horticultural exploitation (Osawa & Inohara 2008).

Measurements of basic population traits, pollen deposition, and fruit-set

During the flowering season, the number of all flowering individuals and flowers within each population, along with flower morph (L- or S-morph), were recorded. Whether the two morphs were equally frequent within a population (i.e., deviations from a 1:1 morph ratio) was investigated with Chi-square goodness-of-fit tests.

To elucidate the role of pollen limitation and morph-ratio variation on female reproductive success, we measured the stigmatic pollen load. In the natural fields, 20 fully-opened flowers were collected from each population in the afternoon and carefully transported to the laboratory the same day. In the laboratory, flowers were dissected, and stigma removed and mounted on a microscope slide in aniline blue staining solution (0.1% aniline blue in 0.1 M K₃PO₄). Under a compound microscope (Olympus), we directly counted the number of legitimate (pollen from the opposite morph) and illegitimate (pollen from the same morph) pollen grains on stigmas of each floral morph based on pollen size differences (L-morph: 18.0 ± 0.1 μm; S-morph: 28.6 ± 0.1 μm; Y. Ojima, unpublished data).

At the beginning of the fruiting season (August), the population mean for fruit-set per flower was measured with the exception of some flowers that were used for the measurements of stigmatic pollen loads. Even if the fruit was set, the reproductive success will strongly depend on fruit predation (e.g., Matsumura & Washitani 2000; Yamamoto et al. 2013). Thus, we continued observations until October that was immediately before avulsion of the capsules.

All field surveys described above were conducted in 2013 and 2014 for KOM and KIN population, respectively.

Population sampling, DNA extraction, and microsatellite genotyping

In each population, 32 plants were sampled randomly (without distinction of floral morph types) from the entire area occupied by each population for genetic analysis. Leaf materials were collected and dried in silica gel. Before DNA extraction, leaves were homogenized with a disposable homogenizer (Biomasher 2; Nippi Co., Tokyo, Japan) to a fine powder. Total DNA was extracted from 40 to 80 mg silica-dried leaf tissue using the grass-fiber filter method (Takakura 2011). The extracted DNA was dissolved in a TE solution and stored at 4⁰C until use.

After the preliminary marker screening, the genotypes of each individual were characterized following seven microsatellite markers that were originally developed for Primula sieboldii E. Morren: ga0161, ga0218, ga0580, ga0691, ga1140, Pri0141, and 2ca135 (Ueno et al. 2003, 2006, 2009; Kitamoto et al. 2005). PCR amplifications and allele-size determination of fragment analysis were performed in accordance with the methods described by Yamamoto et al. (2017).

Population genetic analysis

For all seven microsatellite loci, the absence of linkage disequilibrium (LD) and the presence of null alleles were tested using Genepop v4.2 (Raymond & Rousset 1995). The LD test was verified using a Markov chain method with 1,000 dememorization steps, and 1,000 iterations per batch. Null allele frequencies were estimated by maximum-likelihood estimator based on the expectation-maximization algorithm (Dempster et al. 1977) with the default setting.

The following measures were calculated for each population: number of alleles (A), effective number of alleles (A_E), number of private alleles (A_P), expected heterozygosity (H_E), and inbreeding coefficient (F_IS). Deviations from the Hardy-Weinberg equilibrium were determined by the exact test and permutations. All measurements were calculated using GenoDive v2.0 (Meirmans & van Tienderen 2004). GenoDive was also used to compute the population’s genetic differentiation pairwise F_ST and G’_ST indices (Hedrick 2005), and F_ST was tested for significance using 10,000 permutations.

To estimate genetic structure of P. reinii populations in Hakone volcano, we used the model-based clustering method STRUCTURE 2.3.4 (Pritchard et al. 2000) and non-model-based method principal component analysis (PCA). STRUCTURE analysis was conducted for all samples across the two populations. Under an admixture model with correlated allele frequency, 20 independent simulations were run for each K (K = 1–5) with 5 × 10⁵ Markov chain Monte Carlo (MCMC) steps and a burn-in period of 10⁵ interactions. The most likely value of K was determined by the ΔK method (Evanno et al. 2005) with STRUCTURE HARVESTER 0.6.94 (Earl & vonHoldt 2012). CLUMPAK (Kopelman et al. 2015) was used to average the outputs from multiple STRUCTURE runs and produce the graphical results. The F value, the amount of genetic drift between each cluster and a common ancestral population, was also calculated for each cluster. The PCA analysis was performed using the package adegenet 2.0.1 (Jombart 2008) in R 3.5.2 (R core Team 2018).

To detect a genetic imprint of past population bottlenecks, we first used the heterozygosity excess method (Cornuet & Luikart 1996) implemented within the program BOTTLENECK v1.2 (Piry et al. 1999). This method is suitable to detect very recent and less severe bottlenecks, and has low false positive error rates (Williamson-Natesan 2005). All simulations were done with mutation-drift equilibrium conditions (2,000 replicates) under the stepwise mutation model (SMM), infinite allele model (IAM), and two-phase mutation model (TPM: 70% SMM and 30% IAM). A two-tailed Wilcoxon signed-rank test was used to determine a significant excess of heterozygosity.

We also calculated the M-ratio (Garza & Williamson 2001) for each population using Arlequin (Excoffier et al. 2005). The M-ratio test is considered to have a greater detection power for ancient and moderate-to-severe population declines in comparison with the former method (Williamson-Natesan 2005). M-ratio represents the number of alleles relative to the range in allele sizes. After a severe bottleneck, the number of alleles should reduce faster than the allelic size range, which results in a reduced M-ratio (Garza & Williamson 2001). Thus, the magnitude of the decrease reflects the severity and duration of the reduction in population size, and generally an M-ratio <0.68 is indicative of the presence of a bottleneck (Garza & Williamson 2001).

Finally, we conducted a Bayesian demographic analysis using the R package, Vareff (Nikolic & Chevalet 2014). In contrast to the first two moment-based methods, this coalescent-based approach can examine temporal changes in the effective population size (Ne). The function Vareff simulates prior changes in the effective population size from microsatellite data by resolving coalescent theory and using an approximate likelihood MCMC (Nikolic & Chevalet 2014). After a series of preliminary runs, we used the prior parameter settings for each population (Table S1), following recommendations from Nikolic & Chevalet (2014). We set the estimated mutation rate to 5 × 10^–4 (Estoup et al. 2002) for all loci, and ran each analysis under a two-phased mutation model with a proportion of 0.22 for multiple mutations (Peery et al. 2012), for 10⁵ MCMC steps (NumberBatch = 1,000,000, LengthBatch = 10), sampling every 10 steps (SpaceBatch = 10) with an acceptance ratio of 0.25 (AccRate = 0.25), after burning of 10,000 steps. Estimations of sizes were searched for from sampling time to 5,000 and 500 generations ago.

RESULTS

Population traits and morph ratio

Each population trait is summarized in Table 1. A total of 72 flowering individuals and 99 flowers were found in the KIN population, whereas the KOM population had fewer (52 individuals and 69 flowers). The morph ratio in KIN did not deviate significantly from a 1:1 ratio (L-morph ratio = 0.54), even in the number of flowers (L-morph ratio = 0.52). In contrast, the number of flowering individuals of the L-morph was significantly higher than that of the S-morph in KOM (L-morph ratio = 0.65) and even higher for the number of flowers (L-morph ratio = 0.70).

Pollen deposition

Stigmas of both floral morphs received pollen grains in each population, but the numbers varied greatly between the individuals, ranging from zero to 321. In the KIN population, no differences in stigmatic pollen loads were detected between morphs (Fig. 2a). In addition, the proportion of deposited legitimate pollens was not significantly different between both morphs (Fig. 2c), while the proportions varied greatly among the L-morph stigmas in comparison with the S-morph. This is complemented by the result that the S-morph stigmas received significantly more legitimate pollen grains than the L-morph stigmas (Fig. 2b), implying that S-morph stigmas were pollinated more effectively than the opposite morph.

In contrast, in the KOM population, L-morph stigmas received a significantly greater number of pollen grains than the S-morph stigmas (Fig. 2a). After classifying pollen grains, however, we found no legitimate pollen grains loaded on the L-morph stigmas (Fig. 2b); that is, most L-morph stigmas were covered with a large quantity of illegitimate pollen. Although several S-morph stigmas were legitimately pollinated, similar to other populations (Fig. 2c), there was no significant difference in the number of legitimate deposited pollen grains between the two morphs (Fig. 2b).

Fruit-set

At the population level, fruit-set ratio was much higher in the KIN population than in the KOM population (32.3 % and 14.3 %, respectively) (Fig. 2c). Within a population, both L- and S-morph scored comparable values in the KIN population (37.1 % and 26.7 %, respectively). In contrast, fruit-set of L-morph in the KOM population was less than half of the opposite morph (10.5 % and 27.3 %, respectively). We continued monitoring until October, but no evidence of fruit predation was found in either population, namely fruit-set was almost unchanged throughout the fruiting season.

Genetic diversity

LD between locus pairs was not significant. Although the frequencies of the majority of the null alleles were lower than 0.1, higher frequencies of null alleles were detected on 2ca135 and ga1140 loci in the KOM population (0.227 and 0.114, respectively). As the presence of null alleles may affect the estimation of genetic diversity or differentiation, we excluded the two loci and repeated several analyses to compare results between seven and five microsatellites. This trial revealed no clear difference in the results based on 5 vs. 7 loci (Table S2). Thus, seven loci were used in all analyses described below.

Genetic diversity parameters for the two populations are presented in Table 2. In total, 63 alleles were amplified by seven microsatellite markers, with an average 9.0 alleles per locus. All diversity measurements were slightly higher in the KIN population (A = 6.3, A_E = 3.3, A_P = 2.9 and H_E = 0.652) than in the KOM population (A = 6.1, A_E = 2.6, A_P = 1.9 and H_E = 0.544). The inbreeding coefficient value (F_IS) was positive and comparable between the populations, but each value did not deviate significantly from zero.

Genetic structure and evidence a recent bottleneck

A moderate genetic differentiation was detected between populations (F_ST = 0.115, P < 0.001; G’_ST = 0.286). The STRUCTURE analysis based on the ΔK method indicated that ΔK was 462.5 for K = 2 and ΔK were <3 for other values of K. Therefore, the optimal ΔK for K = 2 showed that the best-fit model for the 64 sampled individuals of P. reinii revealed two clusters (Fig. 3a). Although several admixed individuals were found in each population, all samples formed a clear genetic structure between the two populations. The F values of clusters produced by STRUCTURE analysis were higher in the KOM population (F = 0.186) than in the KIN population (F = 0.086), indicating that the KOM population had undergone a larger genetic drift compared to that of the KIN population. In the PCA (Fig. 3b), the first two axes explained 17.0% and 10.7% of the variances in the experimental data, respectively. The results also distinguished the two populations, suggesting the existence of two genetic units corresponding to each population.

In BOTTLENECK analyses, the two-tailed Wilcoxon signed-rank test provided statistical support (P < 0.025) to the presence of a recent bottleneck in the KOM population under the SMM and TPM, whereas no evidence was found in KIN (Table 2). On the contrary, the M-ratio test indicated that both populations experienced a reduction in population size. The M-ratio values were 0.312 and 0.338 for the KIN and KOM populations, respectively (Table 2). A clear signature of historical population contraction was detected only in KOM via the third method, the Bayesian population demographic analysis. The bottleneck began approximately 1,000 generations ago (Fig. 4a), whereas a gradual decline was settled at least 100 generations ago (Fig. 4b). In contrast, the KIN population seems to have historically had a large constant population size (Fig. 4a); however, recent changes were unclear due to the broad confidence levels (Fig. 4b).

DISCUSSION

Reproductive status and genetic diversity

The observed low reproduction in KOM is congruent with reports that morph-biased populations experience reduced reproduction (Byers & Meagher 1992; Kery et al. 2003; Wang et al. 2005; Pedersen et al. 2016). Given that almost all stigmas were covered with L-morph pollen grains (Fig. 2), it is plausible that frequent self- or intra-morph (i.e., illegitimate) pollination had occurred among the KOM L-morphs. Therefore, our ecological data indicate that the low fruiting success in KOM L-morphs was caused by stigmatic clogging (Yeo 1975) as a consequence of the skewed morph ratio. Because L-morph flowers generally produce greater amounts of pollen grains than S-morph flowers (Richards 2003), it is apparent that the total pollen pool within KOM was occupied by a large amount of L-morph pollen. Similar to our results, previous studies in distylous plants showed higher female reproductive success in the relatively less abundant morph than the dominant morph (e.g., Wyatt & Hellwig 1979; Thompson et al. 2003; Wang et al. 2005; García-Robledo & Mora 2007). Thus, these results may demonstrate negative frequency-dependent patterns of reproductive success in the distylous primrose.

The indices of genetic diversity were relatively high and comparable between the two populations (Table 2), despite the skewed morph ratio observed. In addition, each population exhibited low F_IS levels with no significant deviation from the Hardy–Weinberg equilibrium. These results allow for the conclusion that Primula reinii growing in the volcano had maintained sufficient genetic diversity as a result of outbreeding.

Overall, this study suggests the persistence of distylous self-incompatibility system in the P. reinii populations. Nevertheless, determination of the exact causes of floral morph bias in KOM was not possible based on the limited ecological and genetic data currently available. Because skewed morph ratios are often explained by several biotic and abiotic factors as discussed in the Introduction, there is a need for future studies investigating the ability of selfing and intra-morph mating, maternal fitness differences between morphs, pollinator assemblage, and population demography.

Genetic differentiation and structure

Our molecular analysis showed that genetic differentiation was moderate between the two populations (F_ST = 0.115). Additionally, signs of genetic admixture between the populations were detected in PCA and STRUCTURE analyses (Fig. 3). There are at least two non-exclusive explanations for this: recent lineage divergence and gene flow. According to accumulated geographical surveys, Mt. Kintoki (locality of KIN pop.) and Mt. Komagatake (locality of KOM pop.) formed approximately ca 350–300 ka (Nagai & Takahashi 2008) and ca 27–20 ka (Kobayashi 1999; Nagai & Takahashi 2008), respectively. Formation of the central cone (i.e., Mt. Komagatake) clearly corresponded to the period of the last glacial maximum (LGM; ca 25–15 ka), suggesting that the KOM population was established at least after the last glacial period. The observed high F value (STRUCTURE analysis) and low private alleles in KOM may support a migration scenario that the population experienced a founder effect arising from a post-glacial refugial isolation and subsequent migration from the lowland of the caldera to the high-altitude areas of the central cone during the late Pleistocene and Holocene. Hence, it is plausible that the detected genetic admixture between populations suggests incomplete lineage sorting (i.e., sharing ancestral polymorphism between populations) due to recent lineage divergence.

Given the geographically close relationship between the populations (Fig. 1b), the presence of contemporary gene flow will also be taken into consideration. Because the two populations are severely isolated by a volcanic landform, gene flow mediated by pollen would be a plausible hypothesis. Moreover, in the flowering season we found claw marks, a useful indicator for the pollination services provided by bumblebees (Washitani et al. 1994), on the petals of each population. This may suggest that the bumblebees have a key role in pollination within the Primula reinii populations. Although bumblebees are known as strong-flying insects (e.g., Rao & Strange 2012), previous observations in other Primula species have demonstrated that pollen transfer by bumblebees generally occurs within short distances (e.g., Ishihama et al. 2006). Therefore, we determined that the pollen flow between the populations might occur contemporarily but on very rare occasions. Nevertheless, deciding among the possible explanations for the genetic composition of the primrose in the Hakone volcano is difficult due to the weak evidence based on an insufficient number of loci.

Recent and historical demography

The two tests for a recent bottleneck yielded mixed results (Table 2). Based on the BOTTLENECK analysis, only the KOM population exhibited excess heterozygosity. In contrast, the M-ratio test supported a recent population size reduction in both populations. As mentioned above, however, because these inconsistent results might be attributed to the low statistical power of our sample size (e.g., number of loci or individuals), our results should be interpreted with caution. Nevertheless, such conflicting results often indicate the severity or timing of the reduction in population size (Williamson-Natesan 2005; Marshall et al. 2009; Padilla et al. 2015; Tóth et al. 2019), and were expected due to the differences in power detecting a bottleneck (Peery et al. 2012).

Considering the robust results in KOM, it is likely that the morph-biased population may have undergone more recent and severe bottlenecks in comparison with another population. In theory, the BOTTLENECK analysis can demonstrate population bottlenecks over a period of 0.2–4.0 Ne generations (Cornuet & Luikart 1996). Assuming for KOM population of Ne = 100 (Fig. 4) and a generation time of 2–3 years, it translates into approximately 50–1000 years before the present. On the other hand, a clear sign of recent (within 100 generations) population bottleneck was not found in the Bayesian demographic analysis (Fig. 4). Therefore, based on results from a series of demographic analyses, it is difficult to draw a definitive conclusion on whether recent bottlenecks occur or not, and thus, we defer a final conclusion until more genetic data are available in the future.

Contrary to this, the Bayesian demographic analysis provided strong evidence in support of a historical population bottleneck in KOM inhabiting the central volcanic cone. The first signs of population decline would have occurred 2–3 ka (assuming a generation time of 2–3 years). This timeframe post-dates a climatic warming, known as the Jomon optimum transgression, that occurred approximately 6ka, implying that historical population bottlenecks were likely due to volcanic activities as opposed to climatic events. According to geological records, the last major eruption of the volcano was from the central cone in 2.7–2.9 ka (Kobayashi 1997; Kobayashi et al. 2006), and intermittent phreatic eruptions continued until present-day. Although speculative, these evidences may support the idea that the historical population declines experienced by the KOM could have been associated with repeated eruptive activities in the central cone. Perhaps, the detected recent bottlenecks in KOM are caused by eruptive activities rather than human activities.

On the other hand, the estimated effective population size in the KIN population inhabiting the somma mountains was large and constant in the long term, suggesting that the population has been maintained without suffering from volcanic eruptions occurring in the central cone. Further studies for the lineage divergence and demographic history of P. reinii in this region, using more informative datasets (e.g., single nucleotide polymorphisms), will be valuable because volcanism is one of the key abiotic factors in the plant’s diversification and distribution in Japan (e.g., Yoichi et al. 2017; Nagasawa et al. 2020), located in the Pacific Ring of Fire.

Implication for conservation

Our study suggests that morph imbalances are striking effects on the reproduction of P. reinii population in the short-term. Accordingly, a measure of morph ratio should be given top priority in conservation management of the species, and enhancement of habitat monitoring should be considered as in situ managements to protect remnant individuals and to maintain optimum morph frequencies from horticultural exploitation. Considering the observed negative frequency-dependent patterns of reproductive success, if heteromorphic self-incompatibility is totally strict in P. reinii, the skewed morph ratio in KOM may be improved in the future when regeneration is successful. However, the exact breeding system of the species remains poorly understood. Therefore, in addition to other examinations (e.g., the germination requirements and the effect of storage time of seeds) towards a future ex situ conservation strategy, the levels of within morph fertility and selfing ability should be resolved immediately to evaluate the medium- to long-term risk of extinction in the remnant populations across species distribution ranges.

The two surveyed populations in the Hakone volcano were distinguished by two genetic clusters, suggesting that each population should be divided into a different management unit to maintain evolutionary distinctiveness and ecological viability (Moritz 1994; Frankham et al. 2002). The moderate genetic differentiation and the presence of large amounts of private alleles between the populations highlight this suggestion; thus, artificial inter-population crossing should be avoided in this case. Nevertheless, the lack of samples from other parts of the volcano will influence the estimated genetic structure. Thus, an exhaustive population sampling, including other remnant small population, is required to elucidate the genetic structure and demographic history of P. reinii occurring in the Hakone volcano as is also needed for planning conservation strategies.

To our knowledge, this is the first conservation genetics study on threatened plants in the Hakone volcano, which harbors approximately 1,800 plant species (Tanaka 2008). Thus, the results discussed here will be useful for designing both in situ and ex situ conservation strategies for P. reinii as well as other plants inhabiting the volcano and shed light on the instability of plant populations due to the impacts of volcanism and human activities. Our study highlights the importance of studies in conservation, integrating ecological and genetic approaches to accurately assess the population status of endangered species and draw up effective conservation strategies.

Table 1. Number of blooming plants during the flowering season in each population. Results of χ² goodness-of-fit tests for the similarity between the two morphs.

Number of	L-morph	S-morph	Total	χ^2	P
KIN pop.
Flowering individuals	39	33	72	0.50	0.48
Flowers	51	48	99	0.09	0.76
KOM pop.
Flowering individuals	34	18	52	4.92	0.03
Flowers	48	21	69	10.57	0.001

Table 2. Genetic diversity and detection of a recent population bottleneck of the two Primula reinii populations.

Pop. code	Genetic diversity measurements					P values of Wilcoxon test			M-ratio
	A	AE	AP	HE	FIS	IAM	SMM	TPM	Mean	SD
KIN	6.3	3.3	2.9	0.652	0.068	0.109	0.297	0.813	0.312	0.221
KOM	6.1	2.6	1.9	0.544	0.073	0.296	0.007	0.015	0.338	0.199
total	9.0	2.8	-	0.598	0.070	-	-	-	-

A, mean number of alleles; A_E, mean number of effective alleles; A_P, mean number of private alleles; H_E, expected heterozygosity; F_IS, coefficient of inbreeding; IAM, infinite allele model; SMM, stepwise mutation model; TPM, two-phase mutation model.

Table S1. Prior parameter settings for each population using VarEff software.

Parameters	KIN	KOM	Description
JMAX	4	4	Number of when the effective size has changed
DMAX	10	7	the maximal distance between alleles
NBAR	100	100	prior value for the effective population size
RHOCORN	0	0	coefficient of correlation between effective population size in successive intervals
VARP1	3	3	variance of prior log-distribution of effective population size
VARP2	3	3	variance of prior log-distribution of time intervals
GBAR	10000	5000	number of generations since the assumed origin of the population
Diagonale	0.5	0.5	a smoothing parameter

Table S2. Genetic diversity measurements and population differentiation between the two populations based on the selected five loci (ga0161, ga0218, ga0580, ga0691 and Pri0141).

Pop. code	A	AE	HE	FIS	FST	G’ST
KIN	6.0	3.3	0.637	0.020	0.111	0.249
KOM	5.0	2.1	0.472	0.029

For Image & Figures - - Cliick Here

REFERENCES

Barret, S.C. (1989). The evolutionary breakdown of heterostyly, pp. 151−169. In: Bock J.H & Y.B. Linhart (eds.). The Evolutionary Ecology of Plants. Westvier Press, USA, 614pp. https://doi.org/10.1201/9780429310720

Barrett, S.C. (2002). Evolution of sex: the evolution of plant sexual diversity. Nature Reviews Genetics 3: 274.

Barrett, S.C. & J.S. Shore (2008). New insights on heterostyly: comparative biology, ecology and genetics, pp. 3−32. In: Franklin-Tong V. (ed.). Self-incompatibility in Flowering Plants. Springer, Germany, XLI+313pp. https://doi.org/10.1007/978-3-540-68486-2

Brys, R., H. Jacquemyn & T. Beeckman (2008). Morph-ratio variation, population size and female reproductive success in distylous Pulmonaria officinalis (Boraginaceae). Journal of Evolutionary Biology 21: 1281−1289.

Byers, D.L. & T.R. Meagher (1992). Mate availability in small populations of plant species with homomorphic sporophytic self-incompatibility. Heredity 68: 353−359.

Cornuet, J.M. & G. Luikart (1996). Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144: 2001–2014.

Dempster, A.P., N.M. Laird & D.B. Rubin (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (methodological) 39: 1–22.

Earl, D.A. & B.M. vonHoldt (2012). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources 4: 359–361.

Estoup, A., P. Jarne & J.M. Cornuet (2002). Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. Molecular Ecology 11: 1591–1604.

Evanno, G., S. Regnaut & J. Goudet (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology 14: 2611–2620.

Excoffier, L., G. Laval & S. Schneider (2005). Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evolutionary Bioinformatics 1: 47–50.

Frankham, R., D.A. Briscoe & J.D. Ballou (2002). Introduction to conservation genetics. Cambridge university press, UK.

García-Robledo, C. & F. Mora (2007). Pollination biology and the impact of floral display, pollen donors, and distyly on seed production in Arcytophyllum lavarum (Rubiaceae). Plant Biology 9: 453–461.

Garza, J.C. & E.G. Williamson (2001). Detection of reduction in population size using data from microsatellite loci. Molecular Ecology 10: 305–318.

Hedrick, P.W. (2005). A standardized genetic differentiation measure. Evolution 59: 1633–1638.

Heuch, I. (1979). Equilibrium populations of heterostylous plants. Theoretical Population Biology 15: 43–57.

Hodgins, K.A. & S.C. Barrett (2006). Female reproductive success and the evolution of mating-type frequencies in tristylous populations. New Phytologist 171: 569–580.

Ishihama, F., S. Ueno, Y. Tsumura & I. Washitani (2006). Effects of density and floral morph on pollen flow and seed reproduction of an endangered heterostylous herb, Primula sieboldii. Journal of Ecology 94: 846-855.

Jackson, P.W. & K. Kennedy (2009). The global strategy for plant conservation: a challenge and opportunity for the international community. Trends in Plant Science 14: 578–580.

Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24: 1403–1405.

Kéry, M., D. Matthies & B. Schmid (2003). Demographic stochasticity in population fragments of the declining distylous perennial Primula veris (Primulaceae). Basic and Applied Ecology 4: 197–206.

Kitamoto, N., M. Ueno, Y. Tsumura, I. Washitani & R. Ohsawa (2005). Development of microsatellite markers in Primula sieboldii E. Morren, a threatened herb. Japanese Journal of Conservation Ecology 10: 47–51.

Kobayashi, M. (1997). ¹⁴C ages of pyroclastic-flow deposits from central cones on the western slope of old somma of Hakone volcano, central Japan. Bulletin of Volcanological Society of Japan 42: 355–358.

Kobayashi, M. (1999). Tephrochronology and eruptive activity on Hakone volcano in the past 50 ka. The Quaternary Research 38: 327–343.

Kobayashi, M., K. Mannen, M. Okuno, T. Nakamura & K. Hakamata (2006). The Owakidani Tephra group: A newly discovered post-magmatic eruption product of Hakone Volcano, Japan. Bulletin of Volcanological Society of Japan 51: 245–256.

Kopelman, N.M., J. Mayzel, M. Jakobsson, N.A. Rosenberg & I. Mayrose (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Molecular Ecology Resources 15: 1179–1191.

Marshall, J., B. Kingsbury & D. Minchella (2009). Microsatellite variation, population structure, and bottlenecks in the threatened copperbelly water snake. Conservation Genetics 10: 465–476.

Matsumura, C. & I. Washitani (2000). Effects of population size and pollinator limitation on seed-set of Primula sieboldii populations in a fragmented landscape. Ecological Research 15: 307–322.

Meirmans, P.G. & P.H. van Tienderen (2004). GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes 4: 792–794.

Meeus, S., O. Honnay, R. Brys & H. Jacquemyn (2012). Biased morph ratios and skewed mating success contribute to loss of genetic diversity in the distylous Pulmonaria officinalis. Annals of Botany 109: 227–235.

Ministry of the Environment (2019). The Japanese red lists 2019. Retrieved from http://www.env.go.jp/press/files/jp/110615.pdf

Moritz, C. (1994). Defining ‘evolutionarily significant units’ for conservation. Trends in Ecology & Evolution, 9, 373–375.

Nagai, M. & M. Takahashi (2008). Geology and eruptive history of Hakone volcano, central Japan. Research Report of the Kanagawa Prefectural Museum Natural history, 13: 25–42.

Nagasawa, K., H. Setoguchi, M. Maki, H. Goto, K. Fukushima, Y. Isagi, Y. Suyama, A. Matsuo, Y. Tsunamoto, K. Sawa & S. Sakaguchi (2020). Genetic consequences of plant edaphic specialization to solfatara fields: Phylogenetic and population genetic analysis of Carex angustisquama (Cyperaceae). Molecular Ecology 29(17): 3234–3247. https://doi.org/10.1111/mec.15324

Nikolic, N. & C. Chevalet (2014). Detecting past changes of effective population size. Evolutionary Applications 7: 663–681.

Osawa, T. & S. Inohara (2008). Understanding current status and factors of degradation of threatened plants in the Hakone region of the Fuji-Hakone-Izu National Park −Investigation using park volunteer research data. Japanese Journal of Conservation Ecology 13: 179–186.

Padilla, D.P., L.G. Spurgin, F.A. Fairfield, J.C. Illera & D.S. Richardson (2015). Population history, gene flow, and bottlenecks in island populations of a secondary seed disperser, the southern grey shrike (Lanius meridionalis koenigi). Ecology and Evolution 5: 36–45.

Pedersen, J.L., S.E. Macdonald & S.E. Nielsen (2016). Reproductive ecology of the distylous species Houstonia longifolia: implications for conservation of a rare species. Botany 94: 983–992.

Peery, M.Z., R. Kirb, B.N. Reid, R. Stoelting, E.L.E.N.A. Doucet-Bëer, S. Robinson, C. Vásquez-Carrillo, J.N Pauli & P.J. Palsbøll (2012). Reliability of genetic bottleneck tests for detecting recent population declines. Molecular Ecology 21: 3403–3418.

Piry, S., G. Luilart & J.M. Cornuet (1999). BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. Journal of Heredity 90: 502–503.

Pritchard, J.K., M. Stephens & P. Donnelly (2000). Inference of population structure using multilocus genotype data. Genetics 155: 945–959.

R Core Team (2018). R: A Language and Environment for Statistical Computing. Retrieved from https://www.r-project.org/

Rao, S. & J.P. Strange (2012). Bumble bee (Hymenoptera: Apidae) foraging distance and colony density associated with a late-season mass flowering crop. Environmental entomology 41: 905–915.

Raymond, M. & F. Rousset (1995). Genepop (Version-1.2) –population-genetics software for exact tests and Ecumenicism. Journal of Heredity 86: 248–249.

Richards, J. (2003). Primula. Timber Press, USA.

Takakura, K. (2011). Improved method of plant DNA extraction using glass-fiber filter. Bunrui, 11: 139–149.

Tanaka, N. (2008). Plants in Hakone [translated from Japanese]. Newsletter of the Kanagawa Prefectural Museum of Natural History, 14: 21.

Thompson, J.D., S.C. Barrett & A.M. Baker (2003). Frequency–dependent variation in reproductive success in Narcissus: implications for the maintenance of stigma–height dimorphism. Proceedings of the Royal Society of London. Series B: Biological Sciences 270: 949–953.

Tóth, E.G., F. Tremblay, J.M. Housset, Y. Bergeron & C. Carcaillet (2019). Geographic isolation and climatic variability contribute to genetic differentiation in fragmented populations of the long-lived subalpine conifer Pinus cembra L. in the western Alps. BMC Evolutionary Biology 19: 190.

Ueno, S., Y. Tsumura & I. Washitani (2003). Development of microsatellite markers in Primula sieboldii E. Morren, a threatened Japanese perennial herb. Conservation Genetics 4: 809–811.

Ueno, S., N. Kitamoto, R. Ohsawa, Y. Tsumura & I. Washitani (2006). Nine additional microsatellite markers for Primula sieboldii E. Morren. Conservation Genetics 6: 1063–1064.

Ueno, S., Y. Yoshida, Y. Taguchi, M. Honjo, N. Kitamoto, I. Washitani, R. Ohsawa & Y. Tsumura (2009). Development of 120 microsatellite markers for Primula sieboldii E. Morren for linkage mapping. Conservation Genetics 10: 1945–1952.

van Rossum, F. & L. Triest (2006). Within-population genetic variation in the distylous Primula veris: Does floral morph anisoplethy matter in fragmented habitats?. Perspectives in Plant Ecology, Evolution and Systematics 7: 263–273.

van Rossum, F., S.C. De Sousa & L. Triest (2006). Morph-specific differences in reproductive success in the distylous Primula veris in a context of habitat fragmentation. Acta Oecologica 30: 426–433.

Wang, Y., Q.F. Wang, Y.H. Guo & S.C. Barrett (2005). Reproductive consequences of interactions between clonal growth and sexual reproduction in Nymphoides peltata: a distylous aquatic plant. New Phytologist 165: 329–336.

Washitani, I. (1996). Predicted genetic consequences of strong fertility selection due to pollinator loss in an isolated population of Primula sieboldii. Conservation Biology 10: 59–64.

Washitani, I., M. Kato, J. Nishihiro & K. Suzuki (1994). Importance of queen bumble bees as pollinators facilitating inter-morph crossing in Primula sieboldii. Plant Species Biology 9: 169–176.

Washitani, I., F. Ishihama, C. Matsumura, M. Nagai, J. Nishihiro & M. Nishihiro (2005). Conservation ecology of Primula sieboldii: synthesis of information toward the prediction of the genetic/demographic fate of a population. Plant Species Biology 20: 3–15.

Wilcock, C. & R. Neiland (2002). Pollination failure in plants: why it happens and when it matters. Trends in plant science 7: 270–277.

Williamson-Natesan, E.G. (2005). Comparison of methods for detecting bottlenecks from microsatellite loci. Conservation Genetics 6: 551–562.

Wyatt, R. & R.L. Hellwig (1979). Factors determining fruit set in heterostylous bluets, Houstonia caerulea (Rubiaceae). Systematic Botany 4: 103–114.

Yamamoto, M., K. Horita, D. Takahashi, Y. Murai & H. Setoguchi (2018). Floral morphology and pollinator fauna of sister species Primula takedana and P. hidakana in Hokkaido Island, Japan. Bulletin of the National Museum of Nature and Science Series B 44: 97–103.

Yamamoto, M., K. Kurata & H. Setoguchi (2017). Conservation genetics of an ex situ population of Primula reinii var. rhodotricha, an endangered primrose endemic to Japan on a limestone mountain. Conservation Genetics 18: 1141–1150.

Yamamoto, M., M. Yasui, H. Setoguchi & K. Kurata (2013). Conservation of an endangered perennial herb, Primula reinii var. rhodotricha. Journal of Nature Restoration and Conservation 6: 23–29.

Yoichi, W., X.F. Jin, C.I. Peng, I. Tamaki & N. Tomaru (2017). Contrasting diversification history between insular and continental species of three-leaved azaleas (Rhododendron sect. Brachycalyx) in East Asia. Journal of Biogeography 44: 1065–1076.