Predicting Genetic Demography of Rear-edge Populations of Hemerocallis Middendori: A Test for Climate Effects in Last Glacial Maximum and Holocene Optimum

Quaternary climate changes signicantly impacted population demography of temperate organisms by shifting their distribution. Notably, the rear-edge populations are considered to be more prone to these changes, but empirical studies showed the southernmost fragmented populations of Japanese woody plants to harbor high genetic diversities due to their southern glacial refugia origin. Therefore, the impacts of Holocene climate warming on rear-edge populations have been rarely demonstrated. For the better interpretation of genetic backgrounds of temperate plants, the association of paleodistributions under both icy and warm climates with species-specic demographic changes is required. A perennial daylily Hemerocallis middendori (Asphodelaceae) is widely distributed in temperate and cool-temperate zones of East Asia. In Japan, larger populations are found in central ranges, while few small populations survive on harsh rock walls in southernmost regions. We focused on these variable populations and aimed to predict the population demography in relation to past climate changes by statistically combining population genetics with paleodistribution modeling.


Abstract
Background Quaternary climate changes signi cantly impacted population demography of temperate organisms by shifting their distribution. Notably, the rear-edge populations are considered to be more prone to these changes, but empirical studies showed the southernmost fragmented populations of Japanese woody plants to harbor high genetic diversities due to their southern glacial refugia origin. Therefore, the impacts of Holocene climate warming on rear-edge populations have been rarely demonstrated. For the better interpretation of genetic backgrounds of temperate plants, the association of paleodistributions under both icy and warm climates with species-speci c demographic changes is required. A perennial daylily Hemerocallis middendor i (Asphodelaceae) is widely distributed in temperate and cool-temperate zones of East Asia. In Japan, larger populations are found in central ranges, while few small populations survive on harsh rock walls in southernmost regions. We focused on these variable populations and aimed to predict the population demography in relation to past climate changes by statistically combining population genetics with paleodistribution modeling.

Results
EST-SSR analysis of 737 individuals from 41 populations revealed 6 regional population groups. Four groups widely dominating the northern-central ranges harbored high genetic diversity, whereas genetic divergence within the groups was low. However, two groups at the southwest edge were geographically and genetically isolated and showed the lowest genetic diversity. Estimated paleodistributions showed a decrease of suitable range during Holocene climate optimum in comparison with LGM, and a sole variable of habitat suitability in the Holocene optimum was able to predict genetic diversity across its range.

Conclusions
We concluded that habitat fragmentation and population decline in relation to the climate warming during the Holocene optimum and interspeci c competition with woody plants resulted in genetic isolation and impoverishment of the rear-edge populations.

Backgrounds
Climate change caused by glacial-interglacial cycles during the Quaternary period produced signi cant impacts on species distribution and played essential roles in shaping the diversity patterns of contemporary populations [1,2]. These climatic oscillations drove north-south and highland-lowland migration cycles of plant distribution, which strongly affected their demographic history [3][4][5]. Therefore, current genetic patterns within species re ect their responses to expansions and contractions of their habitats associated with the climate cycles. Particularly, rear-edge populations, de ned as those populations residing at the current low-latitude margins of species distribution ranges in cases of Northern Hemisphere [6], are assumed to be in uenced the climate cycles strongly because of the marginality of the species that is characterized by ecological, geographical and genetic dimensions [7].
One of the hypotheses about rear-edge populations is that these populations harbor higher genetic diversity, or are genetically divergent to other populations. This is because these populations have endured and persisted the climatic oscillations for a long time in isolated glacial refugia [1,2] with the accumulation of unique genetic variation, compared to leading northern populations, which have low genetic diversity as a result of founder effects [8]. If these populations maintain stable population size during postglacial periods as a result of the discrepancy between ecological and geographical margins [7], for example, they harbor higher genetic diversity even now. Rear-edge populations also adapt locally to cope with warm and dry conditions, which they have faced [9,10], and these adaptations may be a key for species survival in the underway global warming. Thus, rear-edge populations have been considered to be important for the long-term conservation of genetic diversity and for the investigation of phylogenetic history and evolutionary potential of species [6,11].
As an alternative hypothesis, the rear-edge populations may be genetically impoverished, more prone to population decline, and nal extinction. These populations may be glacial relicts, barely survived in microrefugia [12] under the in uence of ecological, geographical and genetic dimensions such as genetic drifts and low gene ow as a result of isolation or population decline due to the postglacial habitat fragmentation.
Climate changes are still important as determinants of the genetic characteristics of rear-edge populations. Particularly in the case of rear-edge populations, rising temperatures driven by climate warming could amplify drought stress and thus constrain the growth of plants [13], hindering thereby the colonization to new suitable habitats [14] and increasing the risk of population decline. In fact, population decline in rear-edge populations in response to climate change has been reported [15][16][17]. These con icting hypotheses about the rear-edge population indicate that populations demography in relation to climate change can be de ned as immoderate conditions not only as a cold period but also as a warm period.
Population demography of temperate trees in the Japanese Archipelago in relation to the climatic oscillations have been studied intensively in the last two decades [18][19][20][21][22]. These studies showed that populations in southwestern Japan tend to harbor higher genetic diversity because glacial refugia were located in the south. Whereas lower genetic diversity and unstable demography were observed in northeastern populations. These results suggested that the spatial distribution of glacial refugia during the last glacial period strongly in uenced the current genetic diversities of temperate trees. However, the impacts of climate change vary by species and distributed regions [23], and population demography can take idiosyncratic responses even for plants in the same community [24]. Hence, intensive research of plant species with diverse life forms and species characters is essential to test the alternative demographic hypotheses on rear-edge populations.
Hemerocallis middendor i Trautv. et C.A.Mey (Asphodelaceae), is a perennial daylily, characterized by a day ower with an orange-colored perianth undergoing anthesis in the early morning [25]. Pollinator fauna and the presence or absence of self-compatibility of the species have not been investigated, but it may show similar pattern with related species as H. fulva L. var. longituba (Miq.) Maxim which is known to be pollinated usually by swallowtail butter ies [26] and self-incompatibility [27]. The seeds are black and globular and are thought to scatter by gravity. This species occurs widely in temperate and cooltemperate zones of East Asia; central and northeastern China, northern and central Korea, southeastern Siberia, Sakhalin, southern part of the Kuril Islands, and northern and central Honshu and Hokkaido Islands of the Japanese Archipelago [25,[28][29][30] (Fig. 1). In Japan, this species is classi ed into four intraspeci c taxa based on morphological and ecological differences [31,32]; var. middendor i (northern Japan), var. esculenta (widespread in central Japan), var. esculenta f. musashiensis (lowland form in central Honshu), and var. exaltata (insular taxon in Tobishima and Sado Isls). The habitats of H. middendor i are variable, including low elevation slopes facing the sea in colder regions, such as Hokkaido and northern Honshu, high-elevation grasslands, wetlands and rock walls in central Honshu, suggesting that this species has a wide climatic tolerance as long as the habitats harbor full sun conditions. In the southernmost ranges of this species, the populations inhabit dry rock outcrops, which prevent the invasions of competitive temperate trees, and thus they are isolated from neighboring populations. The geographic location, local environments, and extent of spatial isolation from other populations may indicate that these limited southern populations have served as rear-edge populations (e.g., those in Oki Isls., Fig. 1) in the archipelago. Therefore, the wide-ranging H. middendor i is considered to be a suitable species to investigate the historical effects of palaeoclimates on the demography of regional populations, and particularly those of rear-edge populations.
The aim of this study is to determine the respective roles of extreme climate stages of late Quaternary in shaping the population demography of H. middendor i in the Japanese Archipelago. In particular, our investigation focuses on rear-edge populations, which can provide opportunities to test whether they have served as glacial refugial populations as reported from many temperate trees, or whether they represent the glacial relicts in which interglacial climates limited the population diversity and demography. In this study, we rst performed a phylogeographic analysis to de ne regional populations based on geographic distribution and intraspeci c genetic structure. In order to understand the respective impacts of glacial and interglacial climates, population genetic diversities and reconstructed demographic histories of regional populations were then compared among regions and statistically combined with the paleodistributions predicted by an ecological niche modeling.

Results
Phylogeographic structure of H. middendor i complex The log-likelihood of the data L(K) and ΔK were plotted against the number of clusters (K) assumed in the STRUCTURE analysis (Fig. 2). The L(K) continuously increased with increasing number of K, while ΔK showed the highest value at K = 2 and then a secondary peak at K = 6, indicating the hierarchical genetic structure within the genotype data. The clustering pattern at K = 2 showed a north-south differentiation with a phylogeographic break at the central Honshu Island (Fig. 3), which was considered to be the uppermost structure. An increase in the number of populations (K = 3-6) detected further regional groups, separated by seaways or geographic features. Taking into account the geographic distribution of genetic clusters, we assumed K = 6 to be one optimal number of clusters, and the corresponding genetic structure was biologically interpretable (Fig. 4). Relatively large areas of Hokkaido and the regions to the north of the central Honshu Islands were dominated by blue, light blue, yellow, and orange clusters, respectively, while the southernmost populations in Kyoto Pref. and Oki Isl. were almost xed to each red and green cluster (Fig. 4). An insular population (Sado Isl.) was unique in terms of having an admixed genetic origin (Fig. 3,4). In accordance with the STRUCTURE clustering, UPGMA network (Fig. 5), the NJ network ( Fig. 6) and neighbor-net ( Fig. 7) based on D A distances detected geographic structuring of the populations, although bootstrap values for the phylogenetic networks were mostly low.
Considering the concordant genetic structure among the different analysis methods and geographical distribution, we de ned six regional populations (HK, TO, CH, KI, KY, OK in Table 1, S1), except that an admixed population ("Sado") was excluded from any of the groups due to unclear group assignment.

Genetic diversity among regions and taxa
Although no signi cant differences among regional populations were detected for most genetic variation indices (except H S ; P = 0.018), there were geographic trends that diversity values of AR (Allelic richness [33]), Ho (observed heterozygosity), Hs (within population gene diversity) were highest in TO and lowest in KY (Table 1). Relatively high values of pairwise F ST were detected in southwestern groups of KY and OK (0.214 and 0.117, respectively), while F ST (Fixation index) of central-northeast regional population pairs in HK and TO and KI showed lower levels (Table 2). However, geographic trends were more evident for two summary statistics (PAR; private allelic richness [34] and M ratio [35]) that are more sensitive to population demographic change, showing a pattern of lower levels from north-eastern to southwestern regions (Table 1). Notably, M ratio of the regional population OK was equal to the common threshold (M = 0.68) for the genetic bottleneck suggested by Garza and Williamson [35]. Genetic variations were hierarchically partitioned to different layers (i.e., genetic differentiation among regional groups (a) or among taxa (b), among populations, among individuals, within individuals) by AMOVA (Table 3). Comparing the percentage of variations partitioned to the highest layers, the allele frequency differences assigned among regional groups was 10.3% (Table 3a), which was higher than that among taxa (5.7%) ( Table 3b). The results of AMOVA, together with the STRUCTURE clustering patterns, indicate that population genetic differentiation of H. middendor i complex can be explained by geographic isolation, but not taxonomy.  Fig. 8c), KI (p = 0.024, R 2 = 0.217, Fig. 8d). In contrast, a signi cant correlation was rejected for the TO group (p = 0.370, R 2 = 0.009, Fig. 8b), and the levels of genetic differentiation in this group were exceptionally low (less than 0.05), suggesting an extensive gene ow within the region. The number of populations within regional groups of KY and OK was small, so only pairwise distances were shown as dots, which were plotted with regression lines of the other population groups in Fig. 8e. The genetic distances between populations in KY ranged between 0.1 to 0.2, which were relatively high considering the short geographic distances among populations (at most 18 km apart).

Demographic analysis
The log-transformed posterior distribution of the ratio of N C and N A showed the regional differences in population size change in MSVAR analysis (Fig. 9a). The analyses showed that all populations had decreased their population size (i.e., the peak values were below zero), with the population size ratio smallest in KY followed by OK. The posterior distribution of the time since decline represented a sharp peak at approximately 5kya in KY, while no clear peaks were detected in other regional populations (Fig. 9b). Temporal size change estimated by VAREFF showed similar demographic trends (Fig. 9c), that the population size of KY and OK have been smaller than other regional populations. The intense signal of population declines was detected in KI and TO around 5kya. Of note, these predicted ages may not be entirely reliable because of the large size per each time step in the posterior range of population, probably due to the relatively small number of markers used in the analyses.

Paleodistribution modeling
In the MaxEnt analysis, we obtained the mean area under the curve (AUC) value of 0.851 (SD: 0.013).
This is above the value of 0.7 that is considered acceptable for model-based predictions [36]. The Holocene optimum distribution of H. middendor i was predicted based on the assumption that suitable climatic zones existed mainly in Kuril Islands, a coastal area of Hokkaido and Tohoku district, and also higher altitude zone in southern Hokkaido and Honshu Isls. (Fig. 10b), which is similar to the predicted current distribution (Fig. 10a). Meanwhile, the predicted LGM (last glacial maximum) distribution of H.
middendor i showed that suitable habitats mainly existed in more southerly regions, including the central to south-western Japan, where the species are currently absent (Fig. 10c). These predictions indicate that the distribution range shifts would have occurred after the LGM period, which involved the retreats to northern areas and higher altitudinal zones in southerly regions.

Paleodistributions control the extant population genetic diversity
Arcsin-transformed (asin) He (expected heterozygosity or gene diversity [37]), AR, and PAR were modeled using generalized linear mixed model (GLMM) as a function of habitat suitability during the LGM and Holocene optimum periods. In all genetic diversity indices, the AIC with a combination of a sampling scale of 5 km-6kya and 5 km-21kya were the lowest among the candidate models (Table S2 in Supplementary le 1). As the results of model selection, the model that considered the habitat suitability of 5 km-6kya as a sole explanatory variable was the best among the candidate models for the statistics of asin He and AR (Table S3 and S4 in Supplementary le 1). The estimated coe cients of the xed term of climatic suitability in Holocene were both signi cantly positive with 9.96E-03 (SE: 1.73E-03) and 1.55E-03 (SE: 3.57E-04) for asin He and AR (P = 1.55E-08 and 1.80E-05), respectively (see also Fig. 11).
Regarding the modeling of PAR, a null model attained the lowest AIC among the candidates, meaning that the habitat suitability was unable to predict PAR of the current populations.

Discussion
Genetic differentiation of H. middendor i complex The results of AMOVA (Table 3) and STRUCTURE clustering (Fig. 4) showed that population genetic differentiation of H. middendor i complex mostly re ects geographic factors and disagrees with taxonomy, which is in agreement with previous reports [25,38]. Possible explanations for the discrepancy with intraspeci c taxonomy include that (1) generation of the taxonomic entities was so recent that the genetic markers failed to detect allele frequency differences, which may apply for the case of f. musashiensis (a taxon recognized by a subtle habitat difference [38]), and (2) lineage admixture generated a novel phenotypic population, which is assumed for an insular taxon of var. exaltata (in Sado Isl.) with a complex ancestry (Fig. 4). Therefore, except for one particular case of var. exaltata, our molecular data support a hypothesis that geographic isolation was the primary force to drive the genetic differentiation of populations in the Japanese Archipelago.
Central-peripheral pattern revealed in H. middendor i in the Japanese Archipelago.
A combination of population genetic analysis, demographic analysis, and paleodistribution models revealed the regional variation in population responses to past climate changes. Below, we discuss the histories of H. middendor i populations by focusing on each of three geographic regions, that are central, northern leading-edge, and rear-edge populations.
Central populations such as TO, CH, and KI were located in the center of distribution in both the LGM and following Holocene warming (Fig. 10). The relatively stable demography of these populations (higher M values in Table 1, Fig. 9a) indicates that the central regions would have functioned as long-term refugia for this species, allowing them to maintain the higher genetic diversities (Table 1). One mechanism to shape long-term refugia in this region may be the presence of high-elevation mountain ranges (> 1,500 m a.s.l.), in which plants can shift their distributions upwards and downwards to stay in the same mountain areas [4,39]. The current genetic differentiation was the lowest in the central regions (Table 1) and was similar in TO group with F ST = 0.036 and non-signi cant IBD pattern (Fig. 8b). The nding suggests that the effect of random drifts is mostly cancelled by counteracting gene ow among neighboring populations in the region, where the species is commonly distributed in various habitats [30].
Genetic diversity of H. middendor i in the Japanese Archipelago was higher in central and lower in northern and southern peripheral regions; a pattern follows the prediction from the "central-peripheral hypothesis" [40,41]. The species showed such pattern because the peripheral populations, particularly the rear-edge populations, have lost genetic diversity in relation to population size decline (Fig. 9), which was driven by habitat fragmentation and demographic stochasticity [42][43][44]. The northern populations are currently seen in coastal and lowland areas of Sakhalin and northern Hokkaido Isl. [30]. Palynological records showed that the territories were once covered by boreal coniferous woodland and dry steppe vegetation during the LGM [45], and it was during the post-glacial period when the temperate plants became dominant [46]. In line with the vegetation history, the ENM (ecological niche model) predictions suggested that the expansion of suitable habitats of H. middendor i took place after the LGM periods, indicating the current northern populations are the leading edges. Slightly lower genetic variations in the HK group (Table 1), compared with the southerly TO, may re ect recent demographic stochasticity associated with post-glacial range expansions toward the north. However, since Hokkaido was connected to the continent by land during the last glacial period, there may have been gene ow between populations in both regions during the last glacial period. More extensive sampling will be needed to verify whether the founder effect really worked. The two regional populations in OK and KY can be considered as the rear-edge populations in H. middendor i in Japan, as they are both located in southernmost ranges as isolated populations, having lower genetic diversity and experienced demographic bottlenecks ( Table 1, Fig. 9). Considering the decrease of suitable range and current geographic isolation predicted by ENM, the rear-edge populations may be glacial relicts of the ancestral populations that once expanded the distributions southward during the glacial period (Fig. 10), as suggested in other cool-temperate plants (e.g., a relictual conifer population of Thuja standishii in Oki Islands) [47].

Holocene optimum controlled the demography of the rearedge populations
Which of the climatic extremes of the glacial maximum (cold and dry condition) or climate optimum (warm and wet condition) controlled the population genetic demography has been a controversial topic in understanding the rear-edge plant populations [6]. In this study, the demographic inference revealed that the declines of rear-edge populations occurred in the Holocene epoch (Fig. 9). Also, GLMM model selection showed that the current genetic diversity (He and AR) of H. middendor i can be explained by habitat suitability in the warmest stage, but not by suitability in the glacial climate (Fig. 11). The positive correlation between genetic diversities and 6kya habitat suitability suggests that the lowest levels of genetic diversity of rear-edge populations were likely shaped by population bottlenecks due to climate warming. This pattern is contrary to many other trees in temperate zones in the Japanese Archipelago, in which current genetic diversities re ect the refugial distributions during glacial period(s) and rear-edge populations are characterized by higher genetic diversity [18][19][20][21][22].
The contrasting genetic consequences can be explained by several reasons. The rst consideration is differential life history traits among plant species. Woody plants, which generally have longer life span or gene ow distance than herbs including H. middendor i, are less susceptible to the drift effects [48]. Moreover, perennial herbs are generally more cold hardy than woody perennials, because they can endure winter seasons with their winter buds in the ground [49]. The fact that the current distribution of H. middendor i extends to high latitude zones (~ N54°), where most temperate trees are absent is indicative that the species' northern limits may be controlled by temperature during the growth season. Hence, under the glacial climates, the species may have retained larger-size populations in the northern regions than other temperate woody species did. Second, the rear-edge population of perennial herbs is considered more prone to competitive exclusion than woody plants [50]. In addition to the direct effects of climate warming, inter-speci c competition with surrounding warm-temperate trees in the rear-edge regions would have been in uential for the survival of the light-demanding H. middendor i, which is currently con ned to rock outcrops in these regions. The increased genetic differentiation levels in KY and OK groups would be genetic evidence for severely limited gene ow due to habitat fragmentation in the rear-edge populations (Fig. 8). Therefore, we propose that combined effects of (1) glacial persistence in northern regions such as in northern Hokkaido Isl. and (2) rear-edge population declines due to climate warming and competitions with warm-temperate trees would have led to the unprecedented demographic response of H. middendor i to climate changes.

Conclusions
We conclude that habitat deterioration caused by climate warming after the glacial period resulted in genetic isolation and impoverishment of the rear-edge populations of H. middendor i. Although our genetic analysis showed that the species persisted in the Holocene optimum in the southernmost regions, it is uncertain whether the edge populations can similarly respond to and survive the anthropogenic climate warming. Population collapses [51,52] and genetic diversity decreases [53] driven by on-going climate warming are being reported in the world. In the case of H. middendor i, most of the rear-edge populations are already limited to rocky walls, in which extreme summer droughts defoliate the plant populations. Moreover, the populations in the KY region are suffering from herbivory damage by deer, the population size of which has suddenly increased due to reduced snowfall and relaxed hunting pressure in the last two decades [54]. Genetic assessment of the rear-edge populations performed in this study provides insights into the management decisions for the effective conservation strategy of genetically unique populations.

Population sampling
We  Table S1 in Supplementary le 1. This species is known to have very small populations in southwest Japan, and regional populations KY and OK used in this study are all populations that remain (The distribution of regional population is derived from the results of STRUCTURE analysis, see results

DNA extraction and genotyping methods
The leaf samples were immediately dried with silica-gel beads and kept in the dark at room temperature until DNA extraction. The leaf materials (ca. 1.0 cm 2 ) were then stored at -80 °C in a deep freezer (MDF-C8V1-PJ, Panasonic, Japan) and were pulverized with a TissueLyser (Qiagen. Hilden, Germany). After polysaccharides removal with wash medium described in Setoguchi and Ohba [58] from the resulting powder, total DNA was extracted using the CTAB (Cetyltrimethylammonium bromide) method [59]. For genetic analysis, we used 12 EST-SSR markers developed for H. middendor i [60]. The total PCR (polymerase chain reaction) reaction volume was 6 µL, which contained approximately 0.

Data analysis
Genetic structure in the Japanese Archipelago We employed a model-based Bayesian algorithm implemented in STRUCTURE 2.3.4 [61] to detect genetic structure within the H. middendor i samples. The STRUCTURE analysis was performed according to the assumption that each individual had admixed ancestral origins and different gene pools retained their correlated allele frequency (along with the "correlated allele frequencies model" described by Falush et al. [62]) and that the sampling locations were designated as prior information to obtain better parameter estimates [63]. Twenty independent simulations were run for each K (K = 1-20), with 100,000 burn-in steps followed by 100,000 Markov chain Monte Carlo steps. We computed the corresponding ΔK ad hoc statistics [64] with STRUCTURE HARVESTER 0.6.94 (http://taylor0.biology.ucla.edu/struct_harvest) [65].
Replicates were combined with CLUMPP [66] and bar plots of assignment probabilities were built using DISTRUCT 1.1. [67] The genetic structure among populations was explored by constructing neighbor-joining (NJ) [68] and UPGMA networks based on D A distance [69] among populations, and the signi cance of interior branches of the tree was evaluated based on bootstrap values using a software POPTREE2 [70]. We also calculated pairwise D A distance among populations using GenAlEx 6.51 [71]. This matrix was then used to construct Neighbor-net using Splits Tree4 [72].
Because both model-based and distance-based clustering analyses revealed six population groups corresponding to geographic distributions (see results), all populations were classi ed to six regional groups; that is HK (Hokkaido), TO (Tohoku district), CH (Chubu district), KI (Kinki district), KY (Kyoto prefecture), and OK (Oki Isls). The following genetic diversity and demographic analyses were performed separately for the population groups.

Isolation by distance analysis
To determine the extent of gene ow within a regional population, the presence of signi cant isolation-bydistance patterns within the six regional populations was investigated by applying the Mantel test with 999 permutations to the pairwise relationship between the log-transformed geographic distances

Calculation of genetic diversity
The levels of within-population genetic diversity were evaluated by calculating four genetic diversity statistics using GenAlEx 6.51: average number of alleles per locus (NA), Ho, and xation index (F IS ). AR and PAR were calculated by HP-Rare [74]. By the calculation of allelic richness using the rarefaction method, which can be used to compare allele diversity between populations with different sizes, the standard sample size was set to more than six individuals in this study. M ratio, which evaluates the effect of a population bottleneck, was calculated for the populations with sample size more than 10 individuals following the method described by Garza and Williamson [35]. We also compared the levels of genetic variation (AR, Ho, He, F IS, and F ST ) among regional populations as de ned above. For testing the effect of population grouping, regional genetic diversity and two-sided P-values after 1,000 permutations were assessed using FSTAT 2.9.4 [75].

Genetic differentiation among regions and taxa
Previous genetic studies of H. middendor i were performed to reveal phylogenetic relationships between H. middendor i populations, using universal molecular markers, i.e. chloroplast DNA [25,76], Random Ampli ed Polymorphic DNA (RAPD) [38] and anonymous single nucleotide polymorphisms (SNPs) from genotyping-by-sequencing [77]. These molecular phylogenetic studies showed that some intraspeci c clades were weakly associated with geography. However, because these studies analyzed only one or a few samples per population, no information is available on genetic diversity and demographic changes of regional populations.
To evaluate the contribution rate of genetic differentiation among regions and taxa, we performed a hierarchical analysis of molecular variance [78]. Calculation of regional F-statistics (F RT , F SR , F ST , F IS, and F IT ) and their signi cance using the permutation approach were carried out using the program GenAlEx 6.51. Genetic diversity was hierarchically partitioned into four levels: (1) among taxa (a) or among regions (b), (2) among populations, (3) among individuals, and (4) within individuals. In the analysis of the region model, we set up seven regional groups, i.e., the six regional groups de ned by STRUCTURE analysis and one group consisted of the Sado Isl. population (which was genetically "admixed"; see results).

Demographic analysis
We used a full-likelihood Bayesian method MSVAR 1.3 [79], which is e cient at detecting population declines and expansions from microsatellite data, which provided neither too weak nor too recent event.
MSVAR is also superior to moment-based methods (the M ratio test and Bottleneck) for detecting changes in population size, independent of time and the severity of the event [80]. For the analysis based on coalescent simulations, we employed 10 microsatellite markers that appeared to follow the stepwise mutation model (excluded two markers: hem19141 and hem15494). We estimated the current and ancestral effective population sizes of regional populations (N C and N A ), and the time since the population size started to change. Prior distributions were set as follows: starting value for current and ancestral population size was 10 3 for all loci, mutation rate was 10 − 4 , starting values for time since decline/expansion for all loci were 10 4 , log 10 scale for the prior means and variance of current size, ancestral size, mutation rate and time since decline/expansion were (4, 1), (4, 1), (-4, 2), and (4.5, 2), respectively. Generation time was assumed to be 10 years, by considering 4 years as the time to maturation and 6 years as continuing reproduction period.
To infer the changes in each regional population size, we used R package VAREFF, an approximate likelihood method with a Monte Carlo Markov Chain approach [81]. VAREFF is especially useful to provide

Paleodistribution modeling
Ecological niche modeling (ENM) is a powerful tool for revealing distribution changes in the past times.
ENM can also provide a source of evidence independent from genetic data, which is particularly pertinent for species that are poorly represented in the fossil records. Once current occurrence data and environmental variables are collated in an ecological niche model, historical species distributions can be approximated by projecting species parameters onto a simulated past climate [82].
Because the results of our demographic analyses showed that effective population size changes happened since or after the last glacial period (see results), our paleodistribution modeling targeted the two historical times that represented either climatic extremes: The Last Glacial Maximum (LGM, about 21,000 years ago) and the Holocene optimum, (about 6,000 years ago). For the generation of ENM for H. middendor i, a maximum entropy algorithm implemented in MaxEnt version 3.3.3 k [83] was used to relate occurrence records to the current climate variables. These records were collected either from our original eld surveys or specimen records downloaded from the GBIF website (https://www.gbif.org/).
After the removal of doubtful records, 271 records were used for model construction. We considered nine bioclimatic variables as predictors, which represented in uential factors in terms of temperature and precipitation controls of the species and did not show signi cant clamping effects in model projections (annual mean temperature, temperature seasonality, mean temperature of the coldest quarter, mean temperature of the warmest quarter, annual precipitation, precipitation of the warmest quarter, precipitation of the coldest quarter, precipitation of the wettest quarter and precipitation of the driest quarter). Current climate and paleoclimate layers during the 21kya-CCSM (LGM) and 6kya-CCSM (Holocene optimum) in 30 arc-sec resolution were downloaded from the Worldclim database (https://www.worldclim.org/) [84]. In MaxEnt analysis, 10,000 background points were randomly generated for the study area (the Japanese Islands and Sakhalin), and other parameters were set as defaults, except regularization multiplier as 2, to avoid model over tting. To evaluate the optimal quality of the prediction, 100 replicated bootstrapping was performed.

Relationship between paleodistribution and genetic diversity of the extant population
The values of genetic diversity indices estimated from microsatellite data and paleodistribution reconstructed from current climatic variables and species distribution records were estimated independently. Therefore, the combination of these two approaches can signi cantly improve our understandings of the past processes, which were shaped through population genetic diversity [85][86][87].
To test which of the glacial or interglacial climates were the limiting factor of the current population genetic diversity of H. middendor i, we examined the correlation between genetic diversity estimates and climate suitability around the focal populations, as predicted by ENM. We analyzed three different genetic diversity measures (AR: allelic richness, PAR: private allelic richness, and asin He: arcsine-transformed gene diversity) of the current populations as a function of habitat suitability during the LGM (21kya) and the Holocene optimum (6kya). The distribution probability of H. middendor i populations during the LGM and Holocene optimum were summed around each population locality, and hereafter, this value will be called "habitat suitability". The projected distribution probability given for each 30-sec pixel within the study area was sampled and summed within multiple-ring buffers using ArcGIS 9.2 Desktop (ESRI, California, USA). Multiple-ring buffers were generated from all sampled locations except for the population number of individuals less than six using QGIS 2.14 (QGIS Development Team). We changed buffer radius from 5 to 300 km (5,10,25,50,100,150,200, 250, 300) as an estimated effect of the xed term may vary according to the spatial scale with which habitat suitability was evaluated [20]. We used a generalized linear mixed model (GLMM) to analyze the genetic diversity indices with two xed terms of habitat suitability at 21kya and 6kya, in which the loci as random intercepts, because diversity values at these loci are likely to be more related to each other than to the values from different loci. To explore the combination of spatial scales best explaining genetic variations in each genetic diversity index, we selected the most parsimonious model with minimum Akaike's information criterion (AIC) value among the candidate models. We used R-package lme4 [88] to obtain AIC values among combinations of the different sizes of sampling radius of 21kya and 6kya (9 × 9 scales). Model selection was performed by Rpackage MuMIn [89]. Parameters of the best models were estimated by R-package lmerTest [90]. All statistical analyses were conducted in R 3.5.2 (R Development Core Team 2018). Ethics approval and consent to participate Not applicable.