Differentiated historical demography and ecological niche forming present distribution and genetic structure in coexisting two salamanders (Amphibia, Urodela, Hynobiidae) in a small island, Japan

Population genetics

Ecological niche models

Occurrence data

Present records of two groups were collected from our field survey in 2015–2018. To prevent spatial autocorrelation of occurrence data, we randomly selected points within 100 m, with ten replicates from others by using “spThin” package (Aiello-Lammens et al., 2015). Consequently, 37 localities of Group A (Fig. 2A) and 17 of Group B (Fig. 2B) were employed in the final models.

Niche equivalency and similarity comparison

Observed niche overlap values for salamanders on Tsushima were calculated by using ENMtools (Warren, Glor & Turelli, 2010) with Schoeners’s D and Hellinger’s I niche similarity metrics. These indices range from 0 (no overlap) to 1 (identical niche models), which predicted the similarity of ecological niche between species (Warren, Glor & Turelli, 2008).

The “Identity test” (also called equivalency test) was used to test whether the ENM of Group A is equivalent to Group B (Warren, Glor & Turelli, 2008). The test creates a null distribution by pooling the occurrence points for both species, randomizing the species identities of the localities, and creating two new samples of the same sizes as the original samples without consideration of suitable habitat to either species. Then, we compared observed niche overlap values to the null distribution of 100 pseudo-replicate niche overlap values by one-side test with an alpha level of 0.05. We considered that these ecological niches were not equivalent if the observed value fell within the bottom 5% of the null distribution. The test was implemented by using the package ENMTools 1.0 (Warren et al., 2021) in R version 4.1.0 (R Core Team, 2021).

We also applied the “background test” in ENMtools version R (Warren et al., 2021) to test for niche conservatism or divergence of two species. The test compared an ENM of Group A against a null distribution generated from random points selected within the geographic range of Group B. We used 100 replicates for the test. The opposite direction, comparing an ENM of Group B to and ENMs generated randomly within the ranges of Group A, also was implemented. Then, the test compared the observed niche overlap value between Group A and B to the null distribution by a two-sided test and alpha level of 0.05. We determined that when the observed niche overlap value between two species was above the 95% confidence interval of the null distribution, it might support niche conservatism. By contrast, the observed value was below the 95% confidence interval, supporting niche divergence (Warren, Glor & Turelli, 2008). In the case of the null hypothesis was supported, the niche overlap might be explained by regional similarities in the habitat available to each group. If the background test is only significant in one direction but not for the remaining direction, we could reject the null hypothesis that similarity (or divergence) between group is less than expected based on the availability of habitat (Warren, Glor & Turelli, 2008).

In addition, we also applied a “PCA-env” framework of Broennimann et al. (2012) that uses kernel density to estimate the environmental distribution of species and the distribution of available environment to quantify niche overlap in a two-dimensional environmental space. Then, the data were used to implement hypothesis tests similar to the approach in Warren, Glor & Turelli (2008) (called ecospat-identity test, and ecospat-background test). Here, we used functions enmtools.ecospat.id and enmtools.ecospot.bg in ENMTools version R (Warren et al., 2021) to conduct the tests. The functions automatically ran principal components analysis to reduce the predictors to a two-dimensional space (Warren et al., 2021) because we employed more than two predictors in the analysis.

Demographic history

Partial cyt b gene sequences (413 bp) were determined for 68 specimens (Group A [n = 48], Group B [n = 20]) and deposited in DDBJ (accession numbers: LC638502–LC638569). Twenty-nine haplotypes and nine haplotypes were observed in 73 individuals of Group A and 34 individuals of Group B, respectively (Table 2). Thirty-one and nine polymorphic sites were detected in Group A (n = 73) and Group B (n = 34), respectively. A3, A24, and [A2, A17, A19] haplotypes were the most (34.2%: 25 of 73 individuals), second (8.2%: six individuals), and third (5.5%: four individuals) dominant haplotypes in Group A, respectively. Also, of 34 Group B individuals, B7, [B3, B6], and B1 haplotypes were the most (23.5%: eight individuals), second (20.6%: seven individuals), and third (14.7%: five individuals) abundant haplotypes, respectively. One (A3) and four (B1, B3, B6, and B7) major haplotypes (frequency within each group >10%) were detected in Group A and Group B, respectively. Twenty and four haplotypes were found from one individual in Group A and Group B, respectively. Genetic diversity was higher in Group A (h ± SD: 0.8702 ± 0.0347; π ± SD: 0.00585 ± 0.00356) than in Group B (0.8520 ± 0.0278; 0.00547 ± 0.00343). Mean number of pairwise differences were larger in Group A (2.4140 ± 1.3256) than in Group B (2.2603 ± 1.2739). One and three differences were the most in Group A and Group B, respectively (Fig. 3). Mismatch distribution test of Group A could not rejected significantly the null hypothesis of population expansion (P > 0.05 in both sudden expansion and spatial expansion models), showed a unimodal shape and fitted the curves of both the sudden expansion (Tau = 1.447, SSD = 0.00035184) and spatial expansion models (Tau = 1.328, SSD = 0.00038387) (Fig. 3). These suggests the occurrence of past demographic expansion in Group A. On the other hands, Group B showed a different result; not suggests the occurrence of population expansion in the sudden expansion model (Tau = 2.910, SSD = 0.02556232, P < 0.05), but suggests that in spatial expansion model (Tau = 2.814, SSD = 0.02246269, P > 0.05) (Fig. 3). In the neutrality tests, Group A showed a significant negative value for both Tajima’s D and Fu’s Fs tests (P < 0.01), which suggested to have been occurred a recent population expansion in this group. On the other hands, results of Group B were not significant for both neutrality tests (P > 0.05) (Table 2).

Genetic landscape-shape interpolation analysis illustrated geographic patterns of genetic diversity in each group (Fig. 2). The relatively high peaks (i.e., high genetic diversity) were detected at the high mountainous areas in southern part of Tsushima Island for Group A (Fig. 2A), while the results of Group B showed the high diversity in northern part of the island (Fig. 2B).

Bayesian skyline plot presented different curves between Group A and Group B in change of effective population size. It had been increased constantly at least ca. 300,000 to 50,000 years ago in Group A (Fig. 4A). There was no drastic change for the effective population size of Group B during ca. 200,000 to 100,000 years ago, but the occurrence of a recent expansion was estimated from ca. 75,000 years ago to the present (Fig. 4B). Group A had a larger effective population size than in Group B.

Ecological niches of two salamanders

The Maxent models presented strong ability to generate potential distribution of salamanders on Tsushima Island particularly, the AUC for training Group A = 0.895 ± 0.006, and testing AUC = 0.775 ± 0.068, and for Group B was = 0.908 ± 0.008, and testing AUC = 0.908 ± 0.231. The potential distribution of Group A resulting from our models covered most of Tsushima Island, but the high suitable area showed a concentration on the southern part of Tsushima Island. On the contrary, the predicted distribution of Group B was fragmented and mostly restricted to the northern region. The overlap area between the two groups was fragmented and mainly occurred in the northern part of the island (Fig. 5).

The top three variables contributing the model of Group A account for a total of 63.5%, including Bio12 (24.5%), NDVI_Jun (22.4%), and Landcover (16.6%). For Group B, the variables related to topography were the most important, and the contribution of Elevation and Slope were equal with 28.8%. The Bio12 had the third largest contributor to the Group B model at 15.3% (Table 3).

The shared high contributing variables of both groups were Bio12 (annual precipitation with 24.5% for Group A and 15.3% for Group B) and Elevation (with 10.6% and 28.8% for Group A and Group B, respectively). Interestingly, the response curve of the shared predictors presented distinctive trends. In particular, Group A preferred a higher annual precipitation area, roughly 2,050 mm, whereas the most suitable annual precipitation of Group B was around 1,900 mm (Fig. 6). The elevation showed peak suitability around 100 m and 80 m for Group A, and Group B, respectively. However, the habitat suitability of Group A decreased gently while that of Group B presented a significant decrease after peaking (Fig. 6).

Historical potential distribution

The projected distribution of the two salamanders on paleo-climate reconstructions showed the fluctuation following the time period, and was contrasting compared to the current distribution. Particularly, the potential distributions of both species were increased presently compared to the Last Interglacial, and Last Glacial Maximum (Fig. 7). The distribution of Group A was mostly restricted in the southern part, while the range of Group B concentrated around the isthmus in the center of the island on the Last Interglacial. For the Last Glacial Maximum, the distribution of salamanders tended to move toward the north, Group A focused roughly on the isthmus, whereas Group B moved greatly to the northern tip of the island. On the other hand, the projection on Mid-Holocene increased compared to the present (Fig. 7). Furthermore, it can be seen that the predicted ranges on Mid-Holocene overlapped significantly at the current time for both groups.

Niche equivalency and similarity tests

Our identity test showed that the two salamanders had nonequivalent ENMs (Fig. 8), with Hellinger’s-based I (P < 0.05), and Schoener’s D (P < 0.05), leading to the rejection of the null hypothesis of equivalency test. The background tests indicated that our null hypothesis could not be rejected due to non-significant for both directions of the comparison (P > 0.05 for Hellinger’s-based I; P > 0.05 for Schoener’s D; Fig. 8). The results of ecospat-identity test and ecospat-background test also showed similar trends. The ecospat-identity test rejected the null hypothesis of equivalent niches (P < 0.05) while both directions of ecospat-background test could not reject the null hypothesis of niche similarity (P > 0.05 for Hellinger’s-based I; P > 0.05 for Schoener’s D; Fig. 9).

Demographic history of two salamanders on Tsushima Island

Both haplotype and nucleotide diversities of Group A was higher than those of Group B in our genetic analysis (Table 2), although number of samples was different between two groups (73 specimens in Group A and 34 in Group B), which was reflected by differences on number of localities (37 localities in Group A and 17 in Group B). ENM analyses also inferred that the suitable habitats of Group A were larger than those of Group B at the Last Interglacial (LIG), Mid-Holocene, and current scenario (Fig. 7), and the results of mismatch distribution analysis (Fig. 3) and neutrality test suggested that Group A have undergone a recent population expansion. Group A seemed to invade (or isolate) into Tsushima Island earlier than Group B (Niwa, Kuro-o & Nishikawa, 2021) and the range is separated into two areas (i.e., the northern and southern areas of Tsushima Island) due to split by the intermediate lowland (Fig. 5A). Further, the genetic landscape-shape interpolation analysis showed that high genetic diversity was detected in the southern part of Tsushima for Group A (Fig. 2A). The two areas are not severely isolated but could provide gene flow between the northern and southern populations of Group A. Thus, Group A retains higher haplotype and nucleotide diversity than Group B. Such high intraspecific genetic diversity may enable Group A much more adapted to the variable climate conditions and elevations than Group B. On the basis of ENM analyses, Group A had a restricted distribution in the Last Glacial Maximum (LGM), but a large highly suitable distribution in the Mid-Holocene (Fig. 7). Further, the highly suitable distribution area under current environmental conditions was continuous within Tsushima Island except for the isthmus between the northern and southern areas (Fig. 5A). These suggest that recent population expansion of Group A rapidly occurred on Tsushima during LGM to the Mid-Holocene, which might be caused by ecological adaptation and its genetic basis, e.g., large genetic diversity in Group A.

On the contrary, Group B has a relatively small distribution range in the current climatic condition (Fig. 7) and our ENM analyses suggested the shrinking of their habitat in the LIG and LGM (Fig. 7). In addition, it seems that Group B has not been split into multiple populations (Fig. 5B). Because Group B was closely related to H. nebulosus from Kyushu Island (Niwa, Kuro-o & Nishikawa, 2021), the ancestral population of Group B is suggested to invade to Tsushima from Kyushu Island relatively recently. If such invasion occurred by small-sized populations, genetic diversity in Group B is expected to be small by the founder effect. Such small genetic diversity will prohibit Group B to expand its range as to cover all the island area.

Results of the Bayesian skyline plot indicated a distinctive population dynamics between two groups; Group A was thought to have increased constantly at least ca. 300,000 to 50,000 years ago, while Group B was considered to have undergone no clear changes in the effective population size during ca. 200,000 to 100,000 year ago and then a recent population expansion during last 75,000 years (Fig. 4). Furthermore, the effective population size of Group A was larger than that of Group B (Fig. 4). These would suggest Group A was more suitable to environment and/or geography of Tsushima Island than Group B, which may have been caused by the earlier invasion (or isolation) into Tsushima in Group A than in Group B.

In the mismatch distribution, Group B was suggested to experience the occurrence of population expansion not in both expansion models, and the frequency of number of pairwise difference was the highest at the three differences (Fig. 3). Also, results of Group B were not significant in the neutrality tests, which means that demographic expansion would not have occurred. On the other hand, occurrence of recent population expansion in Group B was inferred by the result of the Bayesian skyline plot (Fig. 4). Interestingly, Group B has four major haplotypes (B1, B3, B6, and B7), in contrast Group A has only one major haplotype (A3). Group B is thought to have diverged from H. nebulosus relative recently (Niwa, Kuro-o & Nishikawa, 2021), thus it is possible that Group B still has retained ancestral haplotypes to the present. Such ancestral variation may affect the results of mismatch distribution and neutrality tests.

Interestingly, the effective population sizes of both two groups did not decrease around the LGM (ca. 22,000 years ago; Fig. 4), although suitable distribution areas of them shrank at that time (Fig. 7). These mean that population density of two groups was high at the LGM. Salamanders of Hynobius are widely distributed in East Asia and some species occurs on the subarctic zone (e.g., Hynobius retardatus on Hokkaido Island of Japan; Hynobius leechii in the north-eastern China) (Sparreboom, 2014) and up to high mountains at 2,600 m in elevation (Hynobius nigrescens in Honshu Island of Japan; Matsui et al., 2020). Salamanders of Hynobius are cold-resistant animals judged by these facts. Therefore, two groups on Tsushima would have not decreased in the effective population sizes even in the LGM.

Distributional pattern of two salamanders on Tsushima Island

Our results on ENMs of two salamanders from Tsushima Island suggested that both salamanders have their own geographic ranges and unique ecological niches (Figs. 5–7). The potential distribution of Group A was larger (Fig. 5), and suitable range in environmental variables also was broader than those of Group B (Fig. 6), indicating that the Group A had a wider ecological niche of environmental variables compared to Group B. In other words, it means that the Group A was more tolerant than Group B.

The contribution of environmental predictors to models of each group also showed distinctive (Table 3), and the response curves of the shared high contribution variables presented considerably different trends (Fig. 6), indicating that they had different correlation with the set of available conditions. The distribution of Group A was affected by precipitation, vegetation, and land cover, which suggests that Group A has been adapted to climate and biological conditions in Tsushima Island. On the contrary, the distribution of Group B was affected by elevation and slope, which suggests that Group B has selected a given habitat based on topography (i.e., low altitudinal and flat lands). Group B may have been less adapted to the climate and biological conditions in the island because of the relatively recent invasion to the island.

Group A and Group B were found syntopically at nine localities (Fig. 1B), where their larvae inhabited the same mountain streams. Furthermore, two groups were phylogenetically close each other and genetic distance of them was not so large (uncorrected p-distance = 9.2%), although they were not sister species (Niwa, Kuro-o & Nishikawa, 2021). Therefore, ecological habits (e.g., prey item) may be similar with each other. Some biological factors such as competition also may affect the formation of current distributions of two groups. However, the knowledge about competition between two groups is lacked, further studies will be required to clarify factors affecting their distribution patterns.

Our models on historical distribution for the two salamanders suggested a significant contraction during LIG and LGM compared to the current distribution (Fig. 7). Our finding revealed that the distribution of these salamanders might have been affected by the climate change during the Quaternary. The climate in the LIG and LGM were colder and drier than the Mid-Holocene and current conditions (Tsukada, 1983; Takahara & Kitagawa, 2000). Particularity, the global cooling during LGM caused by reduction of sea level and atmospheric CO₂ likely led to the smaller potential distribution for both groups than that of current models (Fig. 7). In which, we also found that the suitable area of Group B tends to move to the northern part of Tsushima Island in the LGM. It can be explained by a latitudinal gradient temperature on Japanese Archipelago at this period (Tsukada, 1983). By contrast, the Mid-Holocene model showed a suitable area larger than the predicted by present model for both groups (Fig. 7), supporting by wetter and cooler to warmer climate (Takahara & Kitagawa, 2000; Lutaenko et al., 2007).

Historical refugia of species enable us to expand knowledge on ecological resilience, migration rates in response to shifting climates, and enhance our understanding of how population may react to future climate change (Wielstra et al., 2010). In our result, the genetic landscape-shape interpolation analysis was relatively consistent with the historical suitable area of ENMs, and which suggested refugia used during glacial ages. In Group A, the high genetic diversity and high suitable distribution at the LIG and LGM were projected in southern part (Figs. 2A, 7), suggesting that the past refugia for Group A have existed on the southern part of Tsushima. On the other hand, the relatively high genetic diversity for Group B could be observed in genetic landscape-shape interpolation analysis in northern part of the island (Fig. 2B). The result of the ENM based on LIG and LGM climate projection also showed that suitable areas of Group B concentrated in the northern part (Fig. 7). The location of historical refugia existed within the present distribution of both species suggested that the current ranges of these salamanders were promoted from their refugia during historical ice ages.

Ecological niche differentiation of two salamanders on Tsushima Island

As expected, we found evidence of nonequivalence for the ecological niche of the two groups of the salamanders on Tsushima Island from both identity test and ecospat-identity test (Figs. 8, 9), presenting a lack of exchangeability ecology between them. Otherwise, the result of the background test and the ecospat-background test showed that the null hypothesis was not rejected, meaning the comparison pairs do not show either niche divergence or conservatism based on the background test. The background test could not reveal the divergence/conservatism of the ecological niches between Group A and Group B (Figs. 8, 9), which tells that the niche difference is not so great between them. In fact, the two groups were found syntopically at nine localities from the central to northern parts of Tsushima (Fig. 1B). Probably, another factor may involve determining the distributional difference in the two salamanders. Furthermore, the environmental variables within the selected background sites for these groups were relatively similar, especially for climate data when the area of Tsushima is only a small island (area ~ 696 km²), probably leading to the accepted null hypothesis in background tests. Moreover, nonbiological factors (e.g., level of resolution, or methods of selection of environmental predictors) might also lead to a weak power in the statistic test (Blair et al., 2013). In the present study, it is worth comparing ecology, breeding habits, and life history among syntopic and allopatric areas of the two salamanders. Niwa, Kuro-o & Nishikawa (2021) reported Group B breed earlier than Group A and such difference in breeding habit enable them to occur syntopically. In a preliminary survey by one of the authors (K Niwa, 2015–2018, personal observation), Group A tends to breed in the fast-flowing stream but Group B does in the slow-flowing one, although sometimes they breed in the same stream. Group A might be more adapted to the stream than Group B, which is supported by the wider niche in Group A than Group B and the longer divergence from lentic-breeding ancestors in Group A than Group B in phylogeny (Niwa, Kuro-o & Nishikawa, 2021). The similar result also was detected in other newt species, such as an endemic newt (Lissotriton boscai) in Iberian Peninsula (Peñalver-Alcázar, Jiménez-Valverde & Aragón, 2021), or Lissotriton italicus and L. vulgaris meridionalis in Italian peninsula (Iannella, Cerasoli & Biondi, 2017), or European plethodontid salamanders (genus Hydromantes) (Ficetola et al., 2018). Niche divergence of two groups could be explained by the heterogeneous habitat in environmental space available to each group (Warren, Glor & Turelli, 2008; Blair et al., 2013).

Generally, the body size of closely related organisms correlates with their distribution, species with larger body size possess the widely distributional range due to their competitive relationship in food and/or optimal habitat (Costa et al., 2008; Penner & Rödel, 2019). However, our results showed a different trend when Group A has significantly wider distribution on Tsushima Island, but the body size of two species is mostly equivalent (K Niwa, 2021, unpublished data). Thus, it can be seen that the ecological niche of the two salamanders on the small island were not clearly correlated with their body size. One possible reason to explain the result is the adaptation capacity of each species due to the different time invaded to Tsushima Island. Group A was isolated on Tsushima Island earlier than that of Group B, ca. 3.5 to 3.2 MYA and ca. 1.5 to 1.4 MYA, respectively (Niwa, Kuro-o & Nishikawa, 2021) and occupied mainly the island, including high mountainous area. However, studies on the interspecies interaction and behavior plasticity in both larva and adult period should be conducted to explore the exact mechanisms of ecological relationship between two groups.

Why do two salamanders have the differential demographic/historical patterns?

The present study showed the differential demographic/historical response between two salamanders, although they were mostly sympatric in the small island. One possible reason to explain these differential response between the two groups could be attributed to the different times of invasion to Tsushima Island. Group A is considered to have invaded the island earlier than Group B (Niwa, Kuro-o & Nishikawa, 2021). The different times of invasion may have promoted the different levels of adaptations to the mountainous/lotic environment of the island (i.e., niche differentiation) and those of demographic variation between Group A and Group B. Group A invaded earlier thus to be better adapted in the island than Group B. However, there are no diverse habitats (e.g., breeding site) available to be completely separated by each group in the small (696 km²) and mountainous (a few plains) (Nagaoka, 2001) island. Therefore, it is possible that the ecological niches were not split between two groups and the current sympatric (syntopic) distribution of two salamanders has been retained.