Population Structure and Evolution after Speciation of the Hokkaido Salamander (Hynobius retardatus)

The Hokkaido salamander (Hynobius retardatus) is endemic to Hokkaido Island, Japan, and shows intriguing flexible phenotypic plasticity and regional morphological diversity. However, to date, allozymes and partial mitochondria DNA sequences have provided only an outline of its demographic histories and the pattern of its genetic diversification. To understand the finer details of the population structure of this species and its evolution since speciation, we genotyped five regional populations by using 12 recently developed microsatellite polymorphic markers. We found a clear population structure with low gene flow among the five populations, but a close genetic relationship between the Teshio and Kitami populations. Our demographic analysis suggested that Teshio and Erimo had the largest effective population sizes among the five populations. These findings regarding the population structure and demography of H. retardatus improve our understanding of the faunal phylogeography on Hokkaido Island and also provide fundamental genetic information that will be useful for future studies.


Introduction
One of the major goals of molecular ecology is to understand the current genetic structure and pattern of genetic diversity of organisms. Moreover, these data are useful for conservation management, because genetic diversity is an essential factor in morphological and behavioral phenotypic variation. Salamanders have been used as model species for studies of genetic structure and diversity [1] because gene flow among populations is restricted by their low dispersal ability [2].
Another Hynobius species, the Hokkaido salamander (H. retardatus) is endemic to Hokkaido Island, the northernmost island of the Japanese archipelago. This species has a number of interesting characteristics, including neoteny [15], temperature-dependent sex differentiation [16], and predator-or prey-induced phenotypic plasticity [17][18][19][20][21]. Furthermore, the reported diversification of several life-history traits among regional populations [22] has motivated phylogenetic and phylogeographic studies based on intensive genetic analyses.
Efforts have begun to elucidate genetic diversity in this species. Matsui et al. [23] showed by allozyme electrophoresis of three H. retardatus populations that their genetic differences did not correlate with morphological differences but with geographic location. In addition, Azuma et al. [24] recently identified three mitochondrial haplotype groups, each with a different geographic distribution. However, current genetic structure among local populations and the evolutionary history of those populations remain obscure owing to an overall paucity of genetic analysis data.
Insight into the biogeographic history of Hokkaido Island has been provided by phylogeographic studies of the Hokkaido salamander. The divergence of the focal species from other salamanders is estimated to have occurred 14-18 million years ago (Ma) [6]. Because this speciation time approximately corresponds to the emergence of the Japanese archipelago during the Miocene, Pliocene, and Pleistocene epochs, this species has experienced unstable geological and climate conditions in its current habitat. Hokkaido Island remained connected to Sakhalin Island and the Eurasian Continent by land bridges until the end of the last glacial period, about 10,000 years ago [25]. Most vertebrates in Hokkaido migrated from continental Eurasia via these land bridges. For example, the dark red-backed vole (Myodes rex), which is endemic to Hokkaido Island, immigrated into Hokkaido from the north by the Middle Pleistocene [26], and the brown bear (Ursus arctos) may have migrated to Hokkaido on at least three separate occasions before the disappearance of the land bridges [27]. Another example is the wood mouse (Apodemus speciosus), which colonized Hokkaido Island when the sea level dropped during glacial periods of the Quaternary [28]. Unlike the potential habitable area of these terrestrial organisms, the habitable area of the Hokkaido salamander is restricted to lentic freshwater environments; as a result, the mobile ability of the species is low, a characteristic that it shares with other species of freshwater biota. Therefore, by studying the phylogeography of this species, insight can gained into the biogeographic history of Hokkaido lentic freshwater biota.
In this study, we examined gene flow and the population structure of the Hokkaido salamander by using recently developed polymorphic microsatellite markers. We found low gene flow between populations and a distinct genetic structure. We expected the genetic information obtained by this study to improve our understanding of the phylogeography of this salamander and the overarching biogeography of the lentic freshwater biota on Hokkaido Island. To gain sufficient volume of DNA, we cultivated these salamanders' eggs in the laboratory. Each salamander egg cluster was placed into a different stock tank (33.4 cm × 20 cm × 10 cm high) filled with 2 L of aged tap water and kept in the laboratory at 4°C. After the eggs hatched, we randomly selected hatchlings for genotyping to avoid choosing ones with the same parents. Genomic DNA for genotyping was extracted from tail-tip tissues of approximately 20 individuals from each locality with a DNA suisui-F kit (Rizo, Tsukuba, Japan) following the manufacturer's instructions.

Ethic statement
The research reported herein was approved by Institutional Animal Care and Use Committee of National University Corporation Hokkaido University. Because subjects of regulation by the committee are mammalians, birds, and reptiles, approval number was not issue for our experiments. The focal species was not enlisted as endangered and protected species in any districts in their distribution area at present, and thus no permission was required for the procedure of this study. All salamanders used in this study were deeply euthanized with 0.05% benzocaine in fresh water and stored individually in ethanol, and used for DNA extraction. This was carried out in accordance with the recommendations in the Hokkaido University Manual for Implementing Animal Experimentation. Genotyping We genotyped 12 microsatellite loci of H. retardatus following to the method described by Matsunami et al. [29]. From the resultant genotypes of individuals, we calculated observed and expected heterozygosities (H O and H E ), the number of alleles (N A ), and the inbreeding coefficient (F IS ) with GENALEX 6.5 software [30]. We performed tests for deviation from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) with Genepop ver. 4.2 software [31][32]. Existence of null alleles often leads to overestimation of F ST [33][34]. We used MICRO-CHECKER V2.2.3 [35] to check for the presence of null alleles.

Population structure analysis
To infer the phylogenetic relationships among the five populations, pairwise genetic distance and a phylogenetic tree based on the D A distance of Nei et al. [36] with 1000 bootstrap iterations were calculated with POPTREE2 software [37]. Genetic population structures were analyzed with STRUCTURE version 2.3.2 [38], which uses a Bayesian clustering method to assign individuals to genetic population units. To explore the genetic structure of each population, we conducted multiple analyses while varying the number of Bayesian clusters (K) from 2 to 9. We also used admixture models for Markov Chain Monte Carlo (MCMC) inference with prior information on the locality of samples (LOCPRIOR) and correlated allele frequency. We ran 1,000,000 MCMC repetitions after discarding the first 100,000 iterations as burn-in, and took 10 repeated simulations for each K estimation. To estimate the optimal K value, we analyzed our results according to the method of Evanno et al. [39], which is implemented in the STRUC-TURE HARVESTER web tool [40]. In this method, log-likelihood values of the 10 repeated runs for each K value and their variances are used to calculate ΔK. The average of each replicate cluster analysis was calculated by CLUMPP version 1.1.2 software [41]. The results of the calculations were visualized with DISTRUCT version 1.1 software [42].
To infer subgroups among the five local populations, we estimated the effective population size of each group and the migration rate among the groups. We used two programs to estimate gene flow: BAYESASS+ version 1.2 [43] estimates gene flow by a genetic assignment method, and MIGRATE version 3.6.11 [44] estimates gene flow and effective population sizes by a coalescent method. Of relevance here is that the genetic assignment method adopted by the BAYE-SASS+ typically measures recent migration rates, whereas the coalescent method adopted by MIGRATE estimates long-term migration rates [45]. BAYESASS+ uses fully Bayesian MCMC resampling to estimate the migration rate (m), allele frequency (P), and inbreeding value (F) as variable. We ran a total of 3,000,000 MCMC iterations and sampled the chain every 2000 iterations, discarding the first 1,000,000 iterations as burn-in. MIGRATE uses coalescent simulations to estimate parameter M, which is proportional to the pairwise migration rate (M = m/μ, where m = migration rate per generation and, μ = mutation rate per generation) and parameter Θ, which is proportional to the effective population size N e (Θ = 4N e μ) of each population. We ran 10 short chains with 50,000 generations, sampled every 20 generations, and 8 long chains with 5000 generations, sampled every 20 generations.

Results Genotyping
We genotyped a total of 103 individuals from five local populations using the 12 microsatellite loci (Table 1; S1 Table). The total number of alleles ranged from 3 to 20 with a mean of 8.667 per locus (S2 Table). An Hret02 allele was fixed in Nopporo population, and alleles of Hret05, Hret06, Hret07, and Hret08 were fixed in the Hakodate population. Microchecker program indicated null alleles in five loci (Hret01, Hret02, Hret03, Hret07, and Hret09) in some populations. Observed heterozygosity ranged from 0.208 to 0.664 with a mean of 0.440 per locus. We conducted 54 HWE tests representing every polymorphic locus/population combination (α = 0.00019 after Bonferroni correction) and detected significant deviation from HWE in four tests of four loci (Hret01 in Teshio, Hret03 in Kitami, Hret07 in Teshio, and Hret09 in Teshio). Because none of these loci showed significant HWE deviation in the other populations, we did not exclude them from the analysis. We also conducted 279 LD tests representing all pairwise comparisons of loci in each population, and detected no significant deviation of LD (α = 0.000036 after Bonferroni correction). To confirm the accuracy and reproducibility of our genotyping, we genotyped randomly selected samples twice; in each population, we found only minor differences in the number of alleles (N A ), observed heterozygosity (H O ), expected heterozygosity (H E ), and inbreeding coefficient (F IS ). We subsequently conducted analyses based on loci having no null alleles. Then, we confirmed that the results were qualitatively similar.

Population structure analysis
Pairwise F ST and Nei's D A distances among populations ranged from 0.162 to 0.348 and from 0.613 to 1.166, respectively (S3 Table). The pairwise F ST and Nei's D A distances between Kitami and Hakodate were largest, suggesting that these populations are highly diverged, whereas Nei's D A distances between Kitami and Teshio were smallest. A neighbor-joining phylogenetic tree was reconstructed with these Nei's D A distances (Fig 2) shows the Kitami and Teshio populations clustered with a high bootstrap probability (88%); the Nopporo and Erimo populations are also clustered, but the statistical support is low (58%).
We also inferred subgroupings of the five populations with distinctive multi-locus allele frequencies by calculating ΔK from K = 2 to K = 9 with the STRUCTURE program (Fig 3). In these calculations, ΔK supports that K = 5 as the most suitable population structure model, indicating that each population is genetically distinct. In addition, the populations were clearly divided into two major groups (Erimo-Nopporo-Hakodate and Teshio-Kitami) at K = 2; thus the population structure analysis result is consistent with the phylogenetic tree reconstruction (Fig 2).
The migration rates calculated with BAYESASS+ program were uniformly very small, ranging from 0.0032 to 0.0050 for all population pairs (S4 Table). Although some of the migration rates estimated with the MIGRATE program were also low, rather high migration rates were estimated among the Erimo, Kitami, and Nopporo populations (Fig 4).

Discussion
We reconstructed a cladogram among local populations of Hokkaido salamander by using microsatellite markers. Although some parts of the phylogeny were consistent with the genealogy based on mitochondrial DNA [24], there were also inconsistencies. In our microsatellite analysis, only the Kitami and Teshio populations were clustered with high confidence, whereas the mitochondrial data analysis placed the Nopporo and Teshio populations into the same haplotype subgroup. The Ishikari-Yufutsu lowland frequently appears to act as a barrier in species distributions, including for taxa such as insects [46] and the Japanese brown bear (Ursus arctos) [47]. Nopporo is located at this lowland (Fig 1). During the last interglacial (130,000-60,000 years ago), the Ishikari-Yufutsu lowland was submerged, preventing migration [48]. Therefore, populations of the Hokkaido salamander probably became divided at that time. The Kitami and Teshio populations are clearly clustered in our STRUCTURE analysis for K = 2, and our microsatellite data are more informative than the mitochondrial data reported previously. Thus, we conclude that the genetic relationship indicated by our microsatellite data reflects the true phylogenetic relationships among all five populations. By considering the two independent estimations of gene flow in combination, we concluded that migration occurred rather frequently in the distant past and much less frequently in the recent past. In particular, relatively high gene flow was indicated between the Nopporo and Kitami populations by the MIGRATE program, which estimates long-term migration rates, even though they are separated by high mountains (Kitami, Yubari, and Hidaka ranges). Landscape genetics studies of amphibians have shown that topography affects gene flow and that differences in altitude (such as mountains and valleys) specifically act as a barrier to gene flow between populations (e.g. [49][50][51][52]). Our result for these H. retardatus populations, however, implies that the varying topography did not act as barrier to gene flow, at least in some cases. Hynobidae species tolerate relatively low body temperatures and are abundant in mountainous areas of the Japanese archipelago [53]. The Hokkaido salamander in particular is well adapted to the cold weather on Hokkaido Island and easily moves from the lowlands into the mountains. Therefore, it may be able to migrate across different landscapes.
The estimated divergence time of H. retardatus is about 18-14 Ma [6], but the populations apparently did not diverge until about 0.72-0.32 Ma [24]. There are several possible explanations for this difference. Populations of the Hokkaido salamander experienced many environmental changes and two population bottlenecks [24]. The consequent demographic fluctuations, including population extinctions and recolonization, may have led to lineage sorting or to shallow divergences or both. Our estimated gene flow among populations may therefore reflect migration and colonization occurring in the distant past, but long after speciation.

Conclusion
We clarified the population structure of the Hokkaido salamander and gene flow among populations. Populations of this species show distinct genetic differences and generally low gene flow. Thus, each population has different genetic features and evolved independently. These findings are applicable to the conservation management of this species. Because the populations are genetically clearly different, it is important to protect all of the local populations. In the future, a more detailed population structure and the evolutionary history of this species should be examined by using more comprehensive data, such as genome-wide polymorphic markers. The causal relationships between genetic background and the phenotypic differences among populations also need to be resolved.
Supporting Information S1