Identification of management units in threatened populations of Arnica mallotopus Makino (Asteraceae) using novel EST-SSR markers

ABSTRACT

Arnica mallotopus is a perennial herb endemic to the snowy regions of Japan. At the southern edge of its distribution, in Kyoto Prefecture, overgrazing by sika deer and decreased snowfall have resulted in the rapid decline of A. mallotopus populations. Therefore, there is an urgent need for a conservation genetic analysis of the remaining local populations. In this study, we first developed 13 EST-SSR markers to evaluate genetic variation in A. mallotopus. The average number of alleles per locus was 5.33. Genetic analysis using these markers showed that the investigated samples were classified into two groups corresponding to landscape structure. One group isolated from a tributary of the Yura River showed a strong population bottleneck signal, likely resulting from founder effects and subsequent drifts. On the other hand, the genetic diversity of the second group in the main distribution along the Yura River was higher and less inbred. Overall, our assessment suggested recognizing the two genetic groups as management units in conservation programs for the threatened populations.

MAIN

Arnica mallotopus Makino is a perennial herb of Asteraceae endemic to Honshu and its adjacent islands in Japan (Ohashi et al., 2016). The species occurs preferentially on rock outcrops along mountain streams that are well watered by snowmelt and are often maintained by natural disturbances such as avalanches and landslides. As a result, the distribution of A. mallotopus in Honshu strongly corresponds to the snowy regions adjacent to the Sea of Japan. At the southern edge of their range in Honshu, the plants are distributed as fragmented populations confined to mountainous areas that locally enjoy heavy snowfall. One such population is known from the Ashiu Forest of Kyoto Prefecture, where A. mallotopus is ranked as a critically endangered species in the Kyoto Prefecture Red List (http://www.pref.kyoto.jp/kankyo/rdb/bio/db/flo0133.html). In the last few decades, the Ashiu Forest has experienced changes in snowfall pattern: the maximum snow depth fell below 50 cm in five years during 2000–2019, while no such years were observed during 1980–1999 (Nakagawa et al., 2020). The overwintering of sika deer in the area is considered to be a result of the decreased snowfall, along with increases in their population size (Li et al., 1996; Fujiki et al., 2011). Since the beginning of the 2000s, it has been reported that the understory vegetation in Ashiu declined due to deer grazing damage (Kato and Okuyama, 2004). Similarly, the A. mallotopus population in Ashiu has declined due to overgrazing by sika deer (Supplementary Fig. S1). Our thorough surveys conducted along the Yura River in September, 1990 and 1991 showed that the species once formed a large and continuous distribution along the Yura River for several kilometers (Fig. 1A). However, in comparable surveys conducted along the Yura River and one of its tributaries in August, 2015 and 2016, approximately 80% of the patches were lost, leaving only small, fragmented populations (Fig. 1B).

Fig. 1.

Map of the study site in Ashiu, Kyoto Prefecture. (A) Stars represent the locations of populations of A. mallotopus in September, 1990 and 1991. (B) Stars represent the locations of populations of A. mallotopus in August, 2016. Pie charts represent distribution of cluster memberships at the population level in the 11 A. mallotopus subpopulations within Ashiu inferred by STRUCTURE (Pritchard et al., 2000). Note that the existence of Pop. C1 (NNS) was not confirmed in the 1990 and 1991 survey because the locality area was out of survey range.

Given the critical situation of the rear-edge population of A. mallotopus in Ashiu, a conservation strategy is urgently required. As human control of the deer population in Ashiu is still underway, ex situ conservation is considered an appropriate way to establish a new population of A. mallotopus from wild seeds in safe gardens. In this process, an assessment of the genetic diversity of the remaining populations is needed, because the composition of the ex situ populations should reflect the genetic variation in natural populations. In addition, the distribution of A. mallotopus in the Ashiu Forest is spatially divided into two areas: the C1 population is isolated along a tributary river and the other is along the main Yura River; Fig. 1B. The latter Yura River population has been further divided into 10 subpopulations (Pop. C2–C11), which were rapidly fragmented within the last three decades. To determine the management units of the A. mallotopus plants in the wild, the genetic structure of the remaining populations should be revealed. However, there are currently no useful genetic markers for this species to obtain population genetic information.

In this study, therefore, we first developed novel polymorphic simple sequence repeats in expressed sequence tag (EST-SSR) markers for A. mallotopus. The established markers were checked for their variability using three populations on Honshu, which showed that 13 markers displayed allelic polymorphisms. We then applied these markers to the threatened populations of A. mallotopus in the Ashiu Forest to elucidate the spatial genetic structure and genetic diversity of the wild populations.

To develop the EST-SSR markers, total RNA was extracted from leaf tissues of A. mallotopus (collected in Kyoto Prefecture; a representative specimen, TNS 774946) using the Agilent Plant RNA Isolation Mini Kit (Agilent Technologies, Santa Clara, CA, USA). A cDNA library was prepared following the Illumina protocol (Illumina, San Diego, CA, USA), and 100-bp paired-end short-read sequencing was conducted using the Illumina HiSeq 2000. The 23,718,722 raw reads were quality-trimmed (quality parameter limit = 0.03; five terminal nucleotides from the 5′ end were removed due to their relatively low quality), and the trimmed reads were assembled with a default setting using CLC Genomics Workbench version 7.5.1 software (CLC bio, Aarhus, Denmark) (DRA Accession No. DRR169028). The assembly resulted in 62,737 contigs with an N50 of 965 bp. Microsatellite regions (≥ 8 dinucleotide repeats, ≥ 8 trinucleotide repeats) were searched using MSATCOMMANDER (Faircloth, 2008). In total, we identified 328 SSR motifs and selected 96 primer pairs as candidates based on repeat numbers (≥ 8 dinucleotide repeats, ≥ 10 trinucleotide repeats). One of three different tag sequences (5′-CACGACGTTGTAAAACGAC-3′, 5′-TGTGGAATTGTGAGCGG-3′, 5′-CTATAGGGCACGCGTGGT-3′ (Takara Bio, Kusatsu, Japan)) was added to the 5′ end of each locus-specific forward primer. Reverse primers were tagged with the PIG-tail sequence 5′-GTTTCTT-3′ (Brownstein et al., 1996).

To assess the amplification and polymorphism of the microsatellite regions, 19, 22 and 23 individuals of A. mallotopus from three populations were used (Pop. A, B and C1, respectively; Supplementary Table S1). Sixteen individuals of A. unalaschcensis Less. var. unalaschcensis (Pop. D) were used to check the cross-species transferability of the primers. The vouchers of each species and locality are listed in Supplementary Table S1. Silica gel-dried leaf materials were pulverized with TissueLyser (QIAGEN, Hilden, Germany). After polysaccharides were removed from this powder using the wash medium described in Setoguchi and Ohba (1995), total DNA was extracted using the cetyltrimethylammonium bromide method (Doyle and Doyle, 1987). The total PCR volume of 6 μl contained approximately 0.5 ng DNA, 3 μl of 2× QIAGEN Multiplex PCR Master Mix, 0.01 μM forward primer, 0.2 μM reverse primer and 0.1 μM fluorescently labeled universal primer. PCR amplification was performed with an initial denaturation at 94 ℃ for 5 min followed by 33 cycles at 94 ℃ for 1 min, 55 ℃ for 45 s and 72 ℃ for 45 s, and then by a final extension at 72 ℃ for 8 min. Amplified products were loaded onto an ABI 3130 auto-sequencer (Applied Biosystems, Foster City, CA, USA) using the GeneScan LIZ-600 size standard (Applied Biosystems), POP7 polymer (Applied Biosystems) and 36-cm capillary array. Fragment size was determined using GeneMapper software (Applied Biosystems).

To evaluate the polymorphisms of the EST-SSR markers, genetic diversity indices (N_A: number of alleles; H_O: observed heterozygosity; and H_E: expected heterozygosity) were calculated using GenAlEx ver. 6.503 (Peakall and Smouse, 2006). Significant deviation from Hardy–Weinberg equilibrium was tested with 1,000 randomizations using FSTAT 2.9.4 (Goudet, 1995). The test for the presence of null alleles was performed using MICRO-CHECKER version 2.2.3 (Van Oosterhout et al., 2004).

From the 96 primer pairs tested, 13 loci were successfully amplified and showed allelic polymorphisms in A. mallotopus (Tables 1 and 2). For the A. mallotopus populations, the number of alleles (N_A) and mean values of H_O and H_E are shown in Table 2. When all populations were combined, the N_A ranged from three to eight with an average of 5.33; the mean values of H_O and H_E were 0.246 and 0.645, respectively. For each population, the N_A ranged from one to seven with an average of 2.70; the average values of H_O and H_E were 0.252 and 0.322, respectively. The MICRO-CHECKER analysis indicated the existence of null alleles at one locus (ar2522) in Pop. A, but no such signatures were detected in the other two populations of A. mallotopus. In all samples (n = 64), significant deviation from Hardy–Weinberg equilibrium was shown at all loci, indicating a population structure effect. On the other hand, in each population, the observed heterozygosity did not significantly deviate from Hardy–Weinberg equilibrium at any locus (P > 0.05). The versatility of the markers was confirmed by the finding that twelve loci were successfully amplified in the other species, A. unalaschcensis var. unalaschcensis (Supplementary Table S2).

Table 1. Characteristics of the 13 EST-SSR markers developed for Arnica mallotopus

Locus name	Primer sequences (5’-3’)	Repeat motif	Allele size range (bp)	Ta (℃)	Fluorescent label	BLASTX top hit description	E-value	DDBJ accession no.
Ar10332	F:TGTGGAATTGTGAGCGGCATGTTCCTCCATTGAAGAAGC R:GTTTCTTTCAACGCTTGTGACTCTTGC	(ATC)12	369–401	55	VIC	uncharacterized protein LOC110880153 [Helianthus annuus]	1.00E-19	ICRF01000001
Ar18013	F:CACGACGTTGTAAAACGACAGTGGCCCAAAGAAACGAAAG R:GTTTCTTGCTTTCATCCTTCCCGTTCG	(AT)14	276–302	55	6-FAM	ribosomal protein S2 [Olea europaea subsp. cuspidata]	2.00E-167	ICRF01000002
Ar25439	F:TGTGGAATTGTGAGCGGGGAGTTAGCCTTCAAGATTCGG R:GTTTCTTAATGCGAAGAAACGGTGGAC	(AG)9	289–313	55	VIC	hypothetical protein E3N88_07687 [Mikania micrantha]	3.00E-131	ICRF01000003
Ar10397	F:CACGACGTTGTAAAACGACAGGATCTGTGAAGAGACGGG R:GTTTCTTCTCACTAGACCGCCAGACC	(AG)19	241–269	55	6-FAM	probable serine/threonine-protein kinase At1g54610 [Helianthus annuus]	0	ICRF01000005
Ar1548	F:TGTGGAATTGTGAGCGGACGTTAAATGGCAAGTTGATCG R:GTTTCTTTGATCACAACCATATACCGGC	(AT)9	351–367	55	VIC	No significant similarity found.		ICRF01000006
Ar16393	F:CTATAGGGCACGCGTGGTCAGGGTCTTTGCCTTCTTCC R:GTTTCTTATGGGTGCAAGAAGGTTCAC	(AAT)16	169–226	55	NED	uncharacterized protein LOC110920925 [Helianthus annuus]	0.00E+00	ICRF01000007
Ar21935	F:TGTGGAATTGTGAGCGGCGCAGGTGAAGTAAACGAGG R:GTTTCTTAGGTATGGAAATCAGGGAGTCC	(AT)20	199–258	55	VIC	protein ELF4-LIKE 4-like [Lactuca sativa]	2.00E-63	ICRF01000008
Ar265	F:CTATAGGGCACGCGTGGTTCTCCTTCGTCAGCATCCTC R:GTTTCTTGTCAACCATCCGCATTTCCC	(AAT)13	212–257	55	NED	hypothetical protein E3N88_10590 [Mikania micrantha]	0.00E+00	ICRF01000009
Ar2522	F:CTATAGGGCACGCGTGGTGCATTGTCATCCTCACAGCC R:GTTTCTTTGACGTGAGACAGCCATTTG	(AG)13	141–173	55	NED	No significant similarity found.		ICRF01000010
Ar10317	F:TGTGGAATTGTGAGCGGGGATCTCTGTGCTGCGTTTC R:GTTTCTTGCTGTTCATCATACAACGACAC	(ATC)12	122–164	55	VIC	SNF1-related protein kinase catalytic subunit alpha KIN10-like [Helianthus annuus]	0	ICRF01000011
Ar15054	F:CTATAGGGCACGCGTGGTACCACCCGCTACCAAGTAAC R:GTTTCTTTTGCGCGAAGATCAAAGAGC	(AT)14	253–281	55	NED	alcohol dehydrogenase 7 [Artemisia annua]	2.00E-27	ICRF01000012
Ar25155	F:CACGACGTTGTAAAACGACGAATTAATCCAGACGGCGGC R:GTTTCTTACCAACACCTTCATGAACCG	(ATC)12	340–376	55	6-FAM	uncharacterized protein LOC110868322 [Helianthus annuus]	7.00E-63	ICRF01000013
Ar45093	F:TGTGGAATTGTGAGCGGTTGCTTCATCGACACCATCC R:GTTTCTTCTGCAAGTCAACTCTCCAGC	(ATC)12	198–243	55	VIC	hypothetical protein E3N88_05109 [Mikania micrantha]	5.00E-87	ICRF01000014

Note: Both PIG-tails for reverse primers and universal sequences for forward primers are underlined. Only proteins whose E-values were lower than E-05 are shown. T_a, annealing temperature.

Table 2. Genetic diversity statistics for three populations of Arnica mallotopus based on newly developed EST-SSR markers

Pop	Pop. A (n = 19)			Pop. B (n = 22)			Pop. C1 (n = 23)			ALL (n = 64)
Locus name	NA	HO	HE	NA	HO	HE	NA	HO	HE	NA	HO	HE
ar10332	3	0.158	0.317	3	0.435	0.521	1	0.000	0.000	3	0.194	0.518a,b
ar18013	1	0.000	0.000	1	0.000	0.000	1	0.000	0.000	3	0.000	0.665a,b
ar25439	5	0.737	0.659	3	0.043	0.198	2	0.429	0.495	6	0.371	0.672a,b
ar10397	2	0.316	0.488	2	0.043	0.043	2	0.353	0.360	5	0.224	0.612a,b
ar1548	2	0.263	0.229	3	0.304	0.453	1	0.000	0.000	4	0.218	0.478a,b
ar16393	4	0.526	0.569	3	0.087	0.124	2	0.294	0.389	7	0.276	0.760a,b
ar21935	3	0.273	0.566	6	0.391	0.566	2	0.200	0.180	5	0.264	0.589a,b
ar265	3	0.455	0.368	7	0.522	0.697	2	0.389	0.375	8	0.471	0.781a,b
ar2522	4	0.091	0.492a,b	3	0.261	0.396	2	0.227	0.201	7	0.200	0.732a,b
ar10317	3	0.421	0.525	2	0.087	0.159	1	0.000	0.000	4	0.156	0.598a,b
ar15054	3	0.579	0.508	3	0.261	0.468	2	0.000	0.083	4	0.266	0.573a,b
ar25155	2	0.211	0.465	4	0.174	0.508	2	0.409	0.491	6	0.286	0.750a,b
ar45093	3	0.368	0.481	4	0.348	0.571	2	0.174	0.159	5	0.297	0.636a,b
Mean	3.00	0.313	0.410	3.44	0.232	0.333	1.67	0.210	0.222	5.33	0.246	0.645

Note: N_A, number of alleles; H_O, observed heterozygosity; H_E, expected heterozygosity. ^a Significant possibility of presence of null alleles (99% confidence level) detected by MICRO-CHECKER (Van Oosterhout et al., 2004). ^b Significant deviation from Hardy–Weinberg equilibrium (99% confidence level).

Utilizing the newly developed EST-SSR markers, we performed a conservation genetic analysis of the fragmented population of A. mallotopus in the Ashiu Forest, Kyoto. The same methods as described above were used for DNA extraction, PCR amplification and fragment analysis. To investigate the genetic structure of A. mallotopus within Ashiu, we conducted a STRUCTURE analysis (Pritchard et al., 2000) and principal component analysis (PCA). The STRUCTURE analysis was performed under the assumptions that each individual had admixed ancestral origins and that different gene pools retained their correlated allele frequencies (Falush et al., 2003), and the sampling locations were designated as prior information. Twenty independent simulations were run for each K (K = 1–12), with 100,000 burn-in steps followed by 100,000 Markov chain Monte Carlo steps. We computed the corresponding ΔK statistic (Evanno et al., 2005) using STRUCTURE HARVESTER 0.6.94 (Earl and vonHoldt, 2012). The PCA was performed by the R package “adegenet” (Jombart, 2008) in R 3.5.2 (R Core Team, 2018). The results of the clustering analyses showed that the A. mallotopus population in the Ashiu Forest could be differentiated into two groups (see below). Therefore, we calculated the F_ST values between the two groups. In addition, we conducted STRUCTURE analysis and PCA including the samples from all populations (Pop. A, B and C) to clarify the degree of genetic differentiation among localities. The genetic diversity (N_A, H_O, H_E), coefficient of inbreeding (F_IS), allelic richness (A_R) and private allelic richness (PA_R) per locus and the mean values of these terms were also calculated for each group. The genetic diversity, F_ST and F_IS indices were calculated using GenAlEx 6.503; A_R and PA_R were estimated by HP-Rare (Kalinowski, 2005). Moreover, we manually calculated the M ratio (Garza and Williamson, 2001) to estimate genetic bottlenecks.

As a result of the STRUCTURE analysis, we found that while L(K) continuously increased with increasing K, ΔK showed the highest single peak at K = 2 (Supplementary Fig. S2). The genetic clusters at K = 2 corresponded to the two groups of Pop. C1 and Pop. C2–C11 (Fig. 1B). Similar results were observed in the PCA (Supplementary Fig. S3), in which the C1 and C2–C11 populations were separated on the first PCA axis. The individuals from one patch in Nanase (NNS: Pop. C1), along a tributary of the Yura River, and ten subpopulations along the main Yura River (YUR: Pop. C2–C11) were clearly differentiated into the two genetic groups (F_ST = 0.409). The results of clustering analyses including the populations from other regions demonstrated significant differentiation between the two groups in Ashiu, despite the short geographic distance between them (Supplementary Fig. S4, S5). These two groups tended to have lower genetic diversity than the other populations (Table 3, Supplementary Table S3). Focusing on the two groups in Ashiu, the YUR group showed higher genetic diversity than the isolated NNS group (H_E: 0.210 in NNS vs. 0.340 in YUR; H_O: 0.190 vs. 0.210; A_R: 1.60 vs. 2.17; PA_R: 0.27 vs. 0.75). Similarly, the number of monomorphic markers was larger in NNS than in YUR (four in NNS vs. one in YUR, Supplementary Table S3). The inbreeding coefficients were positive in both groups, but higher in the YUR group (F_IS: 0.133 in NNS and 0.348 in YUR). The M ratio of NNS (M = 0.528) was considerably lower than that of YUR (0.769) (Table 3).

Table 3. Mean genetic diversity indices of two groups in Ashiu

Group	n	NA	HO	HE	FIS	AR	PAR	M
NNS	23	1.69	0.190	0.210	0.133	1.60	0.27	0.528
YUR	30	2.62	0.210	0.340	0.350	2.17	0.75	0.769
Pop. A	19	3.00	0.313	0.410	0.217	2.57	1.42	0.628
Pop. B	22	3.44	0.232	0.333	0.320	1.92	0.86	0.654

Note: N_A, number of alleles; H_O, observed heterozygosity; H_E, expected heterozygosity; F_IS, coefficient of inbreeding; A_R, allelic richness; PA_R, private allelic richness; M, M ratio.

These results suggest that the two groups (YUR and NNS) can be recognized as management units for genetic conservation in the Ashiu Forest. The genetic divergence between them was unexpectedly high (F_ST = 0.409), considering the short geographic distance between NNS and YUR (ca. 1 km). In contrast, the subpopulations within YUR, separated by similar geographic distances, showed little genetic differentiation (Fig. 1B and Supplementary Fig. S3). Therefore, the high F_ST value is plausibly the result of factors specifically operating for the NNS population. First, NNS is found as a population isolated from the main distribution areas of the Yura River (Fig. 1), and no other suitable habitats along the tributary river are found today. Hence, the NNS population can be considered to have been established by seed dispersal from the YUR group; after its settlement, it is likely that gene flow between the YUR and NNS groups was negligible due to the absence of suitable habitat between them. Genetic bottlenecks during colonization, and accumulating effects of genetic drift, likely reduced the genetic polymorphisms and ultimately fixed these to one allele at four markers in NNS (Supplementary Table S3). Notably, the genetic bottleneck in NNS appears to have been very severe, such that its M ratio of 0.528 is much lower than the common threshold for bottlenecks (M = 0.68) (Garza and Williamson, 2001).

On the other hand, until at least 30 years ago, the YUR group maintained a continuous distribution for several kilometers along the Yura River (Fig. 1A). The population density was once high, and the distance between subpopulations was short enough to maintain gene flow between neighboring subpopulations in YUR. Such distribution characteristics would have contributed to the higher genetic variation maintained in the population unit of YUR. Although the current distribution in the YUR group is fragmented due to deer grazing, the remaining subpopulations still constitute a single genetic cluster (Fig. 1B). This finding is reasonable, because only a short time has passed since the recent population decline (only about five generations, assuming a generation time of five years (K. M., unpublished)), and therefore the genetic effects of fragmentation may not yet be apparent.

Overall, our genetic assessment distinguished two genetic groups of A. mallotopus in the Ashiu Forest, which probably arose through historical population dynamics in the complex mountain landscape. The level of genetic differentiation between the groups was high, and each group showed particular genetic characteristics. Thus, they can be reasonably recognized as separate management units in conservation programs for threatened wild populations. Furthermore, the genetic variation of each group characterized by EST-SSR markers can be referred to as baseline information in setting up and optimizing ex situ populations in the future.

ACKNOWLEDGMENTS

We are grateful to Ms. K. Funatsu, Dr. K. Horie, Ms. E. Katori, Ms. K. Magota, Dr. D. Takahashi and Mr. S. Yatagai for their efforts in collecting the plant materials. We also thank the editor and two anonymous reviewers who provided many helpful suggestions for improving the manuscript. Research permissions were provided by Fukushima Prefecture and Oze National Park (for Pop. A), Tottori Prefecture and Daisen Oki National Park (for Pop. B), Ashiu Forest Research Station, Kyoto University (for Pop. C), and Hokkaido Prefecture and Daisetsuzan National Park (for Pop. D). This work was supported by the TAKARA Harmonist Fund for S. S., the Nippon Life Insurance Foundation Fund for M. I. I., and the Pro Natura Foundation Japan for the Ashiu Biological Conservation Project and Environmental Research.

REFERENCES

•  Brownstein,  M. J.,  Carpten,  J. D., and  Smith,  J. R. (1996) Modulation of non-templated nucleotide addition by Taq DNA polymerase: primer modifications that facilitate genotyping. Biotechniques 20, 1004–1010.
•  Doyle,  J. J., and  Doyle,  J. L. (1987) A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 19, 11–15.
•  Earl,  D. A., and  vonHoldt,  B. M. (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361.
•  Evanno,  G.,  Regnaut,  S., and  Goudet,  J. (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620.
•  Faircloth,  B. C. (2008) MSATCOMMANDER: detection of microsatellite repeat arrays and automated, locus-specific primer design. Mol. Ecol. Resour. 8, 92–94.
•  Falush,  D.,  Stephens,  M., and  Pritchard,  J. K. (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587.
•  Fujiki,  D.,  Kishimoto,  Y., and  Sakata,  H. (2011) Present distribution of a sika deer population and its impact on vegetation around Mt. Hyonosen, Hyogo Prefecture. Japanese J. Conserv. Ecol. 16, 55–67.
•  Garza,  J. C., and  Williamson,  E. G. (2001) Detection of reduction in population size using data from microsatellite loci. Mol. Ecol. 10, 305–318.
•  Goudet,  J. (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. J. Hered. 86, 485–486.
•  Jombart,  T. (2008) adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405.
•  Kalinowski,  S. T. (2005) HP-Rare 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol. Ecol. Notes 5, 187–189.
•  Kato,  M., and  Okuyama,  Y. (2004) Changes in the biodiversity of a deciduous forest ecosystem caused by an increase in the Sika deer population at Ashiu, Japan. Contr. Biol. Lab. Kyoto Univ. 29, 437–448.
•  Li,  Y.,  Maruyama,  N.,  Koganezawa,  M., and  Kanzaki,  N. (1996) Wintering range expansion and increase of sika deer in Nikko in relation to global warming. Wildl. Cons. Jpn. 2, 23–35.
•  Nakagawa,  H.,  Hasegawa,  A.,  Hayashi,  D.,  Furuta,  M.,  Kishimoto,  Y.,  Miyagi,  Y.,  Ohashi,  K.,  Okabe,  Y.,  Yamauchi,  T., and  Ishihara,  M. I. (2020) Long-term monthly climate data at the forest stations of Kyoto University. Ecol. Res. 35, 733–741.
•  Ohashi,  H.,  Kadota,  Y.,  Murata,  J.,  Yonekura,  K., and  Kihara,  H. (2016) Wild Flowers of Japan. Heibonsha, Tokyo (in Japanese).
•  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.
•  Pritchard,  J. K.,  Stephens,  M., and  Donnelly,  P. (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. R Foundation for Statistical Computing, Vienna, Austria.
•  Setoguchi,  H., and  Ohba,  H. (1995) Phylogenetic relationships in Crossostylis (Rhizophoraceae) inferred from restriction site variation of chloroplast DNA. J. Plant Res. 108, 87–92.
•  Van Oosterhout,  C.,  Hutchinson,  W. F.,  Wills,  D. P., and  Shipley,  P. (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538.