Genetic diagnosis of a rare myrmecochorous species, Plagiorhegma dubium (Berberidaceae): Historical genetic bottlenecks and strong spatial structures among populations

Abstract Distribution of genetic variation over time and space is relevant to demographic histories and tightly linked to ecological disturbances as well as evolutionary potential of an organism. Therefore, understanding the pattern of genetic diversity is a primary step in conservation and management projects for rare and threatened plant species. We used eight microsatellite markers to examine the level of genetic diversity, spatial structure, and demographic history of Plagiorhegma dubium, a rare myrmecochorous herb, populations sampled across northeast Asia and Siberia. We found low within‐population genetic variation associated with historical bottlenecks. Although pairwise F ST values were not much higher than the ones found in similar life form species, STRUCTURE and PCoA revealed a clear broadscale spatial pattern of genetic structure. Bayesian clustering (best K = 6) and PCoA identified three populations that are distinctive from neighboring populations in the Korean peninsula, which suggests potential units for conservation and management plans in Korea. MIGRATE‐N and BAYESASS showed that both contemporary (0.003–0.045) and historical migration rates (2 × e−5−4.6 × e−4) were low. Our findings provide a good example, where genetic considerations should be integrated for conservation and management plans of rare and threatened species.

For example, a microsatellite analysis revealed high inbreeding and recent population bottlenecks associated with habitat destruction in an extremely rare annual herb, Rhinanthus osiliensis (Aavik, Talve, Thetloff, Uuemaa, & Oja, 2017). The magnitude of genetic response to anthropogenic disturbances might depend on life history and ecological traits (Aguilar, Quesada, Ashworth, Herrerias-Diego, & Lobo, 2008;Hylander & Ehrlén, 2013;Young et al., 1996). Negative effects of the disturbances on genetic diversity may be more pronounced in short-lived plant species than in long-lived trees (Young et al., 1996).
Likewise, self-incompatible outcrossing species are more prone to negative genetic effects of the disturbances than selfing plants because obligatory outcrossing plants often use clonal growth under severe disturbances Honnay, Jacquemyn, Bossuyt, & Hermy, 2005).
Demographic processes following ecological disturbances may be inferred through examining patterns of genetic diversity using advanced tools of molecular analysis. For example, a molecular marker analysis revealed significantly lower heterozygosity and higher inbreeding in Pacific jumping mice (Zapus trinotatus), a riparian rodent, after a severe flood (Vignieri, 2010). Inference made through molecular markers is particularly useful for rare species with unknown demographic history as there is no need for destructive sampling.
Furthermore, detailed information of demography and spatial distribution of genetic variation can also be applied to management plans for rare and endangered plant species. Using microsatellite data Straub and Doyle (2009) suggested two genetically distinct management units for the endangered legume (Amorpha georgiana) as the data revealed high genetic divergence between the two units.
Myrmecochory plays an important role for population growth and largely influences spatial distribution of genetic variation in obligatory myrmecochorous plants due to its substantial contribution to the plant migration (Pascov et al., 2015). Because ants disperse plant seeds over short geographic distances (less than 180 m at maximum; Gomez & Espadaler, 2013), gene flow among populations might be extremely limited with the limited migration. In Korea, P. dubium was rarely found in natural habitats and it is listed as an "endangered species" by the Ministry of Environment of South Korea in part due to limited migration ability as a myrmecochorous plant (Ministry of Environment, 2005). The species only inhabits dense deciduous forests mostly along the north side of shady valleys (Ghimire & Heo, 2012;Rhie, Lee, & Kim, 2015). However, the species was removed from the list in 2012 based on the recent population growth observed throughout the nation within the past 10 years (3rd National Environmental Survey, Ministry of Environment unpublished data; Suh & Kim, 2012). There is an urgent need of genetically based diagnosis for the possible biases from the survey of census population size and the expanded distribution in P. dubium.
In the present study, we aim to address (a) the level of genetic diversity in Plagiorhegma dubium; (b) how the species is genetically structured over heterogeneous landscapes; (c) the historical and contemporary gene flow; and (d) whether there have been recent and/or historical population bottlenecks. Using microsatellite markers, we investigated genetic diversity patterns of P. dubium throughout the global distribution (northeast Asia and Siberia). Due to rarity of the species, lower within-population genetic diversity is expected. We hypothesize that populations of P. dubium may be largely structured over varying landscapes with limited migration rates through myrmecochory. Given that the species has not been commonly found until recently (within a few generations assuming generation time of 4 years; Smith, Ronsheim, & Swartz, 1986), we also hypothesized that the population expansion was a recent event after a series of historical bottlenecks.

| Study species
Plagiorhegma dubium Maxim. (Figure 1) is an herbaceous perennial that belongs to family Berberidaceae (Kim, Lee, Kim, & Kim, in press;Rhie et al., 2015). P. dubium is the only species in the genus Plagiorhegma Maxim. since the first discovery of the species (Maximowicz, 1859). Later, the species was moved to genus Jeffersonia Bart. due to its similarity in floral structure with the F I G U R E 1 A photograph image of Plagiorhegma dubium. The photograph was taken by one of the author, Bo-Yun Kim on the day, April 15, 2015 in Daegu, Gyeongbuk, S. Korea species in Jeffersonia (J. diphylla (L.) Perss; Bentham & Hooker, 1862;Ying, Boufford, & Brach, 2008). The two species share several common features in morphology and are genetically closely related (Kim & Jansen, 1996;Kim, Kim, Kim, & Jansen, 2004). Despite the conflicts of taxonomic position of the species, in the present study, we will refer to it as Plagiorhegma dubium following the plant list (http:// www.theplantlist.org/). Plants normally flower in the early spring around March to April and produce capsule fruits once a year (Kim, in press;Rhie et al., 2015). Flowers of the plant have both male and female parts in a single flower with six stamens and one carpel (Ghimire & Heo, 2012). Pollination of P. dubium is not yet surveyed; however, the species is likely pollinated by insects as the most closely species, J. diphylla is an insect-pollinated species (Smith et al., 1986). There is a hypothesis that the species may be selectively autogamous to ensure seed set when there are only limited pollinators in cold springs (Smith et al., 1986). P. dubium may also vegetatively spread through rhizomes under harsh conditions (Smith et al., 1986).
Due to a showy light violet flower and heart shape leaves in mature plants, about 10-30 cm long, the plants have been recognized as a potential ornamental species (Rhie et al., 2015;Song, 2007). The species also has been well known for its medicinal uses. The plant roots have been used as a folk medicine for stomach aches and as phytochemical extract called "berberine" proven to be effective in lowering cholesterol (Bae, 2000;Kong et al., 2004).  (Table 1). Although empirically tested rates of clonal propagation are unknown, it is likely that the species can propagate clonally under harsh environmental conditions (Smith et al., 1986). To avoid collecting multiple samples from clones, we put at least 1.5 m distances among individuals within each population.

| Sample collection and DNA isolation
We obtained all required permits to collect samples in protected areas from the Korean National Park Service under related rules and policies. Vouchers of collections (KNR2015110-KNR2015250; KNR2016061-KNR2016140) were deposited in Hallym University herbarium. Leaf tissues collected were preserved at room temperature in plastic bags with silica-gel desiccants until DNA extraction.
We extracted genomic DNA from the dried tissues using DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. Quantity and quality of DNA isolates were assessed using a NanoDrop ND1000 (Thermo Fisher Scientific; quality cutoff, OD 260/280 ratio between 1.7 and 1.9) and visualized using a 1.5% agarose-gel electrophoresis.

| Microsatellite genotyping
Fourteen microsatellite loci were used to genotype the total of 198 individuals from 10 populations using primers developed and tested for Plagiorhegma dubium (Kim et al., in press; see Table 1 for the sample information). We performed PCRs in a volume of 25 μl containing 1 ul of 10 ng template DNA, 2 μl dNTPs (20 mM), 2.5 μl 10 × PCR buffer containing MgCl 2 (Takara, Japan), fluorescently labeled (FAM, HEX, and NED) forward and reverse primers, 0.5 μl (10 pmol each). PCR cycling conditions were as follows: 5 min predenaturation at 95°C followed by 30 cycles of 1 min at 95°C, 1 min at 50°C, and 1 min at 72°C, followed by a final 10 min extension step at 72°C. The fluorescently labeled PCR products were pooled with a size standard GS500LIZ (Applied Biosystems, USA). ABI 3730XL automated sequencer (Applied Biosystems) was used to separate out the amplified fragments. Using Gene Marker program (version 2.40, Softgenetics LLC), we examined microsatellite profiles with automated scoring. We then manually checked the result to avoid scoring errors. We tested for presence of null alleles using FreeNA for each of 14 loci within each population (Chapuis & Estoup, 2007).

| Data analysis
We tested whether there is a significant deviation from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) within each population for each of the 14 microsatellite markers using Fisher's exact test (Guo & Thompson, 1992) in Arlequin v. 3.5 (Excoffier & Lischer, 2010). Bonferroni corrections were used to adjust for significant alpha (p) values for HWE and LD. Six of 14 markers were with moderate to high frequency of null alleles in more than two populations. The six markers were either significantly departure from HWE or not independent from one another (not in LD r > 0.5; p < 0.0001) in over two or more populations. Some of statistical analyses such as STRUCTURE depend largely on two strong assumptions, that is, Hardy-Weinberg equilibrium (HWE) and marker independence within populations (Pritchard, Stevens, & Donnelly, 2000). Because significant deviation from HWE is often linked to heterozygote excess or deficiency, the six loci would also have great impact not only on genetic diversity pattern but also genetic distance-based genetic analyses, for example, PCoA and Mantle test. Therefore, we only used eight microsatellite markers for all downstream analyses to avoid marker-based biases. We performed clonality test to avoid presence of multiple clones within each population for ge-   Table 1 for abbreviation of population locations and sample sizes. Genetic barriers are depicted by gray lines with alphabetically ranked orders. Robustness of each barrier was estimated by 100 bootstraps and represented by the thickness of each line (barrier a > 75% support; barrier b > 95% support). The remaining barriers were not present as statistical significances are very weak at a support level of <50% by chance (PgenF IS ) was calculated in each population (cut off for PgenF IS < 0.05).
We calculated the following genetic diversity measures: expected & observed heterozygosity (He & Ho), number of alleles (Na) allelic richness (Ne), and inbreeding coefficient (F IS ) using Arlequin v. 3.5 and GENALEX v. 6.5 (Excoffier & Lischer, 2010;Peakall & Smouse, 2012). We estimated pairwise F ST between population pairs across all 10 populations in Arlequin v. 3.5. For significance testing, we applied 1,000 permutations with replacements. Slatkin (1995) suggested to use R ST for highly variable microsatellite markers assuming the mutation rate is strictly under a generalized stepwise mutation model. Given that microsatellites hardly are under a strict stepwise mutation model in practice, R ST is not reliable and even prone to biases (Balloux & Goudet, 2002;Lugon-Moulin, Brunner, Wyttenbach, Hausser, & Goudet, 1999;Meirmans & Hedrick, 2011). Therefore, we estimated and used F ST in the present study.
Population genetic structure was investigated using a Bayesian model-based clustering approach implemented in STRUCTURE v.
2.3.4 (Pritchard et al., 2000). We used the admixture for ancestry model and the correlated allele frequency model (Falush, Stephens, & Pritchard, 2003). Ten independent runs were repeated with 1,000,000 MCMC iterations following 100,000 steps of burn-in for each K from 1 to 12. To infer the number of clusters (K) best explaining the data, ΔK was calculated (Evanno, Regnaut, & Goudet, 2005) using STRUCTURE HARVESTER Web v. 0.6.94 (Earl & vonHoldt, 2012). Ancestry coefficients of the individual genotypes from 10 repeated STRUCTURE runs were summarized using CLUMPP v. 1.1.2 (Jakobsson & Rosenberg, 2007) with greedy option. We visualized the result in DISTRUCT v. 1.1 (Rosenberg, 2004). Using Monmonier's (1973) maximum difference algorithm implemented in Barrier v 2.2 (Manni, Guerard, & Heyer, 2004), we identified significant genetic discontinuities between groups of populations. We calculated Nei's genetic distances (Nei, Tajima, & Tateno, 1983) for Barrier analysis. We estimated 100 distance matrices by bootstrapping over loci in MICROSATELLITE ANALYSER (MSA) v. 4.05 (Dieringer & Schlotterer, 2003). The robustness of barriers was examined using the bootstrapped distance matrices. To investigate the pattern of clustering among individual genotypes, we also performed principal coordinate analysis (PCoA) on Nei's genetic distance estimated for all 198 individuals in GENALEX v. 6.5. We conducted Analysis of molecular variance (AMOVA) in Arlequin v. 3.5 to hierarchically partition the total molecular variance within and between groups identified by STRUCTURE analysis. We compared the values against 1,000 resampled data sets to determine statistical significance.
We performed a Mantel test to assess isolation by distance (IBD) based on pairwise genetic distance (Slatkin' (1 − F ST )) and log-transformed Euclidean distance for all population pairs (Rousset, 1997) in GENALEX v. 6.5. One thousand random permutations with replacement were used to examine if the correlation coefficient (r) between the pairwise F ST and the Euclidean distance was significantly different from zero.
We estimated historical migration rates among populations and long-term effective population sizes using MIGRATE-N 3.6.11 (Beerli, 2008). MIGRATE-N employs coalescent approach to estimate both effective population size θ (4Neμ, where Ne = effective population size; μ = mutation rate) and asymmetric migration rate, ulations over a long time scale (Beerli & Felsenstein, 2001). We used Bayesian approach and ran two replicates of MIGRATE-N with the Brownian motion mutation model. We chose uniform prior TA B L E 1 Summary of genetic diversity indices and location information of collection sites for Plagiorhegma dubium  see Table 1 for abbreviations). We also estimated contemporary migration rate among the same groups used in MIGRATE-N over the last few generations in BAYESASS 3.0.1 (Wilson & Rannala, 2003). BAYESASS was started with a 1,000,000 burn-in followed by 10,000,000 MCMC iterations using a sampling frequency of 2,000.
We used the default setting for mixing parameters. To compare historical and contemporary migration rates, migration rates estimated from MIGRATE-N were adjusted by mutation rate, 10 −4 /allele/generation found in Marriage et al. (2009)

for microsatellite markers in
Arabidopsis thaliana (Chiucchi & Gibbs, 2010). We performed a nonparametric hypothesis testing to assess whether difference between the two migration rates is statistically robust using permutation test with 1,000 replicates in R 3.4.4 (R Core Team 2018).
Historical population bottlenecks were examined using Garza- 3.5 (Excoffier & Lischer, 2010). The M-ratio test (Garza & Williamson, 2001) uses the ratio of the number of alleles to the allele size range because a bottleneck is expected to reduce allele number faster than the range of allele sizes. The M-ratio values were statistically tested by 10,000 simulated replicates assuming an effective size, 40 (2 × average sample size per population; 2 was used for diploid) in Critical_M (Garza & Williamson, 2001; available at https://swfsc. noaa.gov/textblock.aspx?Division=FED&id=3298). We also estimated relatively recent reduction in population size within the past few generations using BOTTLENECK v. 1.2.02 (Cornuet & Luikart, 1996;Piry, Luikart, & Cornuet, 1999). BOTTLENECK infers population bottlenecks based on a significant excess or deficit of heterozygosity relative to equilibrium state (Piry et al., 1999). According to a simulation study, the M-ratio is well equipped for identifying a historical bottleneck whereas BOTTLENECK is better suited for detecting recent bottleneck events (Williamson-Natesan, 2005).
BOTTLENECK was run under infinite allele model (IAM) and stepwise mutation model (SMM) with sign and Wilcoxon's sign rank test for statistical significance. We also investigated whether there is mode shift from equilibrium state in allele frequencies.

| RE SULTS
No scoring errors were identified during the genotyping procedures. We found null alleles in six of the 14 microsatellite loci with the frequency of minor to moderate (0.1-0.3 in loci JD_27, JD_66, JD_77, JD_78, JD_81, and JD_93) in more than two populations.
After filtering the six microsatellite markers with null alleles, all 8 markers did not deviate from HWE nor were in LD after Bonferroni correction (p > 0.05). We found significant excess of heterozygotes in those 6 loci (observed heterozygosity > 0.7). Two of the six loci for three populations were with moderate (0.13) to high frequency (0.30) of null alleles. There were a few identical MLGs found in two populations (DG and GJ). However high values of PgenF IS (average = 0.71) indicate that those identical MLGs were not true clones.
Therefore, we treated all individual samples were not clones for the subsequent analysis. The mean expected heterozygosity per population ranged from 0.14 in BS (see Table 1 for population abbreviations) to 0.32 in GS, whereas the mean observed heterozygosity varied from 0.1 (DOM) to 0.46 (VLT) (  Table 2).
Based on delta K from STRUCTURE results, the number of clusters (K, number of randomly mating subgroups) that best explain the data was 6 (Supporting Information Table S2 and Figure S1). To show the clustering pattern with various number of K groups, we provided bar plots of K = 2-8 (Supporting Information Figure S2). The populations US and VLT were the most genetically distinct, followed by GS and GJ (Figure 2a Table 3). AMOVA also revealed that the most genetic variance was partitioned to among group variation (50%) followed by within individuals (39%; Table 3).
Consistent with the STRUCTURE result, the PCoA largely separated out US and VLT from the other 8 populations on the first PC axis (PC1, explained 29% of molecular variance; Figure 3). GS population was forming a cluster that was differentiated from the rest of populations particularly on the second axis (PC1, explain 15% of molecular variance; Figure 3). Genotypes in GJ were scattered throughout both axes with some level of overlap with GS ( Figure 3). Results of Mantel test showed a significant correlation between genetic divergence (F ST ) and geographic distance (Euclidean distance; r = 0.4, p < 0.05; Figure 4).
There was no evidence of asymmetric migration based on the overlapping values between the two directions in 95% confidence intervals for both long-term and contemporary migration rates (CI; Supporting Information Table S3). All contemporary migration rates  Table S3). The estimated θ in MIGRATE-N analysis ranged from 0.33 to 1.26, which are comparable to the ones observed in a long-lived shrub, Rhododendron spp. (Hsieh et al., 2013) and an annual Helianthus spp. (Kane et al., 2009). The longterm mutation-scaled effective population size was the smallest in US (814, 95% CI 250-4,500) and the largest in the cluster2, a group including population GJ and GS (3160, 95% CI 3150-7,500), with the mutation rate.
From the results of the BOTTLENECK analysis, we did not find evidence of recent population bottlenecks for most populations under both IAM and SMM mutation models from the two significant tests except for US (Table 4) (Table 4).

| D ISCUSS I ON
Most rare and threatened plants may severely suffer from erosion of genetic diversity, inbreeding depression and limited gene flow derived by reduced population size and/or fragmented habitats (Ellstrand & Elam, 1993;Ottewell et al., 2015;Young et al., 1996).
Conservation and management plans for threatened plant species could benefit from genetic considerations (Hedrick, 2001). There are several ways to take advantage of genetic data for conservation including the followings: (a) diagnose the status of genetic erosion; (b) infer changes in demography; and (c) identify evolutionarily significant units. Our genetic analyses suggested that Plagiorhegma dubium has historically experienced population bottlenecks despite the recent population growth. We found two genetically distinctive groups (a population US and a group with GS, GJ) with a unique pattern of genetic diversity, which might be an important unit for conservation in Korea.
Significant departures from HWE associated with significant excess of heterozygotes were found over half of 10 populations for six of fourteen microsatellite loci. The heterozygote deficiency was minor and can be explained by the presence of null alleles. In fact, we found moderate to high frequency of null alleles in two of those 6 loci for three populations. The excess of heterozygotes was apparent, however, unlike the heterozygote deficiency, only a few explanations have been suggested for the excess: small sexual & self-incompatible populations, over dominance, negative assortative mating and asexual reproduction (Stoeckel et al., 2006). Given that biological and ecological  TA B L E 4 Results of tests for the recent and past bottlenecks in Plagiorhegma dubium populations: G-W index is the Garza-Williamson index known as M-ratio, the ratio of number of alleles to allele size range Consistent with our hypothesis, within-population genetic variation of P. dubium was lower (mean He = 0.22; mean Ho = 0.33; Table 1) than the ones found in short-lived perennials measured from microsatellite data (He = 0.55; Ho = 0.53; Nybom, 2004).
Although genetic variation is highly related to breeding systems in plants (Baker, 1959;Barrett, 2003;Charlesworth, 2003), the reduced genetic variation in our study may not be attributable to inbreeding. F IS values were low and even negative in several populations for a few loci, which suggest no indication of inbreeding. In fact, high heterozygosity and low inbreeding rates are often linked with self-incompatible outcrossing systems (Baker, 1959;Barrett, 2003;Charlesworth, 2003). For example, Ulmus pumila, one of selfincompatible elm species showed low inbreeding rate with high heterozygosity (Zalapa, Brunet, & Guries, 2010). However, unlike many self-incompatible outcrossing plant species, we found reduced heterozygosity in P. dubium. One explanation for the low genetic variation would be recent and/or historical population bottlenecks.
Despite the moderate to large long-term effective population sizes Distribution of alleles and genotypes across space and time is not random but the result from joint effects of several evolutionary forces; that is, mutation, selection, migration, and stochasticity (Loveless & Hamrick, 1984;Spieth, 1974). Therefore, mating system, migration, and ecological factors associated with both elements greatly influence the genetic structure of plants (Loveless & Hamrick, 1984). The average among-population divergence of P. dubium (F ST = 0.24) was comparable but not significantly higher than the ones found in similar life form plant species (F ST = 0.31 estimated from microsatellite data; Nybom, 2004; Table 2). It is somewhat conflicting to our hypothesis given the biology of the species. P. dubium is insect pollinated with the capability of self-fertilization and vegetative propagation (Smith et al., 1986). Seeds are migrated mostly by multiple ant species (Jo, Kim, Seo, & Lee, 2015;Smith et al., 1986).
As seeds are migrated over very limited geographic distances (Giladi, 2006;Horvitz & Schemske, 1986), populations are expected to be highly structured in most myrmecochorous species. Furthermore, the species has been considered rare and endangered until recently, which suggests the populations are likely to be small and isolated.
Given that the indirect measure of migration (F ST estimates) assume rather unrealistic island model with no mutation, no selection and particularly equal population sizes, F ST does not directly reflect migration and gene flow (Whitlock & McCauley, 1999). Indeed, both recent and long-term migration rates we estimated were very low by overlapping or close to zero at 95% CI, indicating limited gene flow among groups (Supporting Information Table S3). The limited gene flow would have resulted in population isolation and creating the observed spatial structure. Alternatively, historical population bottlenecks found in all 10 populations might have influenced the spatial pattern of genetic structure (Loveless & Hamrick, 1984).
The population US and the group with GS and GJ were genetically distinct from the neighboring populations despite the geographic proximity (Figure 2a,b). PCoA and Barrier results revealed a similar pattern that separates the three populations from the neighboring populations (Figures 2 and 3). The abrupt shift of allele frequency distribution in the three populations may be due to environmental barriers that were not formerly recognized. It is more likely that demographic changes experienced in the past may be responsible for the spatial structure that do not reflect geographic distances (Table 4). Despite the complex spatial structures observed, there was a significant IBD pattern among the 10 populations ( Figure 4). Without proper genetic considerations, P. dubium was cleared from the endangered species list in Korea based on recent growth in census population sizes and range expansion (Suh & Kim, 2012). Our research revealed that the species only harbors a limited amount of individual-level genetic variation (Table 1). On top of that, the genetic diversity pattern clearly indicated that all populations have experienced recent and/or historical bottlenecks. Population genetics can offer insights for management plans of threatened plant species (Hedrick, 2001;Hughes et al. 2008;Aavik et al., 2017). Our findings may be applicable for practical conservation and management plans for P. dubium. First, US, GS, and GJ populations with unique alleles should be treated as separate units for conservation in Korea.
Additionally, Chinese populations can be strong candidates for a restoration project in Korea since most Korean populations shared alleles and genotypes with Chinese populations. Although the P. dubium population has been growing and expanding in range with conservation efforts in Korean peninsula, the populations should be genetically monitored to avoid inaccurate evaluation.
However, the pattern of genetic diversity and its distribution examined here may be biased to some extent as the study could not investigate populations in North Korea due to inaccessibility.
There is an urgent need to conduct ecological and evolutionary monitoring for P. dubium populations in North Korea. Also, the fate of the populations might heavily depend on genetic variation that is adaptive particularly in radically changing environments, which neutral genetic markers cannot capture. Adaptive genetic variation in phenotypic traits should be examined for the future studies. Our study provides one good example of how genetic considerations can be appreciated as much as other ecological and political factors in management and conservation plans for rare and threatened species.

ACK N OWLED G M ENTS
We are grateful to Dr.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
YD designed the project and achieved supporting grant. Sample collections, laboratory work, and microsatellite genotyping were conducted by BY. SR wrote the manuscript and conducted genetic and related statistical analyses. All authors edited the manuscript and agreed to submit current version of manuscript.

DATA ACCE SS I B I LIT Y
Microsatellite data are available in the DRYAD Digital Repository (https://doi.org/10.5061/dryad.d6396sd).