Introduction

Islands are important centers of biodiversity and high plant endemism, with species evolved after colonization (Kier et al. 2009). Analyses and modeling of populations on islands are important for elucidating evolutionary roles of geographic isolation, smaller population size and further ecological processes (MacArthur and Wilson 1963). Populations on islands are often isolated and small comparing to mainland populations, it is a resultant of colonization from the mainland (Stervander et al. 2015). Their sustenance is largely influenced by the balance between restricted colonization due to isolation from other landmasses and extinction depending on habitat areas, sometimes small on islands (MacArthur and Wilson 1963; Warren et al. 2015). After establishment of initial founding populations, subsequent migration between islands and the mainland is an important process underlying genetic divergence and evolution on islands. The small areas of islands and their isolation from the mainland are also considered to play important roles in evolution, different to the mainland (James et al. 2016). Their isolation restricts migration, so colonizing populations generally have substantially lower genetic diversity than the source populations. Thus, we can expect differences in evolutionary and survival history between populations on islands and those in the mainland to be reflected in the genetic variation of the populations. Comparison of island populations to mainland source populations provide valuable insights regarding the survival, evolution and endemism of island populations and more general evolutionary processes (McGlaughlin et al. 2014; Nakamura et al. 2010). Further, environments on islands and the mainland often differ, which promotes adaptive evolution or speciation on islands (Gillespie 2004). High environmental heterogeneity such as temperature and precipitation are regarded as a major driver of such divergent selection on islands (Juan et al. 2000; Weigelt et al. 2013).

The Fossa Magna element plants are a group of 30 endemic species or varieties of plants endemic to the Izu Islands and the central part of Honshu Island of Japan (Maekawa 1949; Takahashi 1971), located on the boundary of four tectonic plates (Taira et al. 1998). In contrast to typical oceanic islands, such as the Galapagos, Hawaiian and Canary Islands, the Izu Islands, all of which are volcanic, are isolated by deep oceanic seaways (Fig. 1) but geographically close to Honshu (treated as the mainland here), thus many species are present on both the Izu Islands and Honshu, and plants on the Izu Islands do not have high endemism due to moderate isolation. Further the Izu Islands are situated roughly along north-south line extending ca. 25–480 km south of the mainland, so their geographic isolation and species composition gradually changes farther away from the mainland (Takahashi 1971). Some Fossa Magna element plants (such as Alnus sieboldiana, Prunus lannesiana and Swertia noguchiana) are endemic to the central part of Honshu and the Izu Islands, while others (such as Calanthe izu-insularis and Lilium auratum var. platyphyllum) are solely endemic to the Izu Islands. The endemism restricted to the Izu Islands is thought to have formed due to isolation from Honshu (Yamada and Maki 2012; Yamamoto et al. 2017), but migration or colonization pattern for endemic species in both the Izu Islands and Honshu are not known. More specifically, it is not known whether this endemism formed through evolution on Honshu followed by colonization of the Izu Islands or evolution on the Izu Islands promoted by isolation followed by ‘re-colonization’ (Bellemain and Ricklefs 2008; Stervander et al. 2015) of Honshu. Thus, elucidation of processes involved in the islands’ colonization, and comparison of changes in the size of populations of species present on both the Izu Islands and Honshu, which have widely differing areas, could provide valuable insights regarding evolutionary processes and their effects on islands (Warren et al. 2015).

Fig. 1: Species range and locations of the eight analyzed populations of Rhododendron tsusiophyllum.
figure 1

Large and small red dots indicate locations of analyzed populations and distribution records of the species, respectively. The three-letter strings are population codes listed in Table 1.

Full size image

Rhododendron tsusiophyllum Sugim. is a dwarf shrub (<50 cm) and a member of the Fossa Magna element plants. The species is restricted to slopes and ridges in upsides of mountains including volcanos of both Honshu and the Izu Islands in the Japanese archipelago, which are sunny, well-drained and wind-swept, with sparsely distributed trees. As such habitats are restricted on Honshu, it has an isolated island-like distribution even on the mainland (Fig. 1). It produces fine seeds (ca. 4 mg/100 seeds, W. Yoichi, personal observation), so the species presumably colonized the Izu Islands by wind. Rhododendron tsusiophyllum is rare, with restricted distribution and listed as vulnerable in the red list maintained by the Ministry of the Environment of Japan (2019). Especially, the species’ habitat area is <1 km2 on Kozushima and Mikurajima Islands of the Izu Islands. Moreover, all of the populations on the Izu Islands grow on volcanic slopes, and a population on Miyakejima Island is thought to have been extinguished by an eruption on AD 2000 (Kamijo and Hashiba 2003). Thus, the populations on the Izu Islands are expected to be more vulnerable than those on Honshu, because of their small size and volcanic disturbance, and we expect this demographic difference between populations on the two regions to be imprinted within its genome. Although the species is rare, it can be found in various habitats and ranges from cool temperate and subalpine sites at 1000–2000 m altitudes on Honshu, to warm temperate sites at 500–800 m altitudes on the Izu Islands (Tanaka 1999). At its coldest limit, it shares habitats with alpine or subalpine plants such as Rhododendron aureum, Vaccinium vitis-idaea and Picea species, while at its warmest limit broad-leaved evergreen trees or shrubs such as Castanopsis sieboldii and Eurya japonica occur (Tanaka 1999). Its closest relative, R. tschonoskii, which is distributed in northern Honshu and Hokkaido, inhabits areas with cold climates, including sites along the northern limit of R. tsusiophyllum (Yoichi et al. 2019). Morphologies of its branches and leaves slightly differ between Honshu and the Izu Islands (Takahashi and Katsuyama 1992). Thus, we expect populations on the two regions to be isolated by climatically not only by geographically, and further to evolve against climatic differences (Weigelt et al. 2013).

In this study, we reconstructed evolutionary history and characterized demographic history of R. tsusiophyllum, focused on its colonization into islands using double digest restriction-site associated DNA (ddRAD) sequencing. We then tested the role of environmental variation in divergent adaptation between populations on Honshu and the Izu Islands. The aim of this study was to answer the following questions: (i) Does R. tsusiophyllum colonize from mainland to islands, or reversely? (ii) What is genetic structure across islands and the mainland? (iii) What is difference in genetic diversity and population size with its historical changes? (iv) Does there exist habitat differentiation between islands and the mainland? Our results will provide new knowledge about colonization and evolution of plant species, and conservation of genetic diversity on islands.

Materials and methods

Plant material and DNA experiment

We collected sample leaves from up to 16 individuals, with taking distance >20 m avoiding substructure within populations, from eight localities of R. tsusiophyllum and six individuals from three localities of R. tschonoskii, which is distributed in northern Honshu and Hokkaido, for use as an outgroup (Yoichi et al. 2019), and immediately dried using silica gel (Fig. 1, Table 1, Supplementary Table 1). These collections were conducted in 2014–2017.

Table 1 Population code, locality and genetic diversity estimates for each sampled population.
Full size table

Genomic DNA was extracted from the samples using the modified CTAB method (Murray and Thompson 1980) or a DNeasy Plant Mini kit (Qiagen) after treatment with sorbitol buffer (Wagner et al. 1987). Genome-wide SNPs were detected by ddRAD sequencing, using Peterson’s protocol with some modifications (Peterson et al. 2012), as described in Yoichi et al. (2019). Genomic DNA was digested with BglII and EcoRI, then digested fragments of 350–400 bp in libraries were sequenced with 51 bp single-end reads using a HiSeq2000 platform (Illumina, San Diego, USA).

SNPs detection

After removing reads containing low-quality bases and adaptor sequences from the raw data using Trimmomatic v. 0.33 (Bolger et al. 2014), SNPs were called by Stacks v. 1.48 (Catchen et al. 2013) including reads from the two species with the following settings: minimum depth option for creating a stack, -m 3; maximum distance between stacks within an individual, -M 2; number of mismatches between loci, -n 1; gapped option for creating stacks, --gapped --max_gaps 6; and default settings for other options. The SNP genotypes were called for eight populations of R. tsusiophyllum (set 1) and both eight populations of R. tsusiophyllum and three populations, treated as one group, of R. tschonoskii (set 2). For this, the Populations program in Stacks was used with the following restrictions: less than four minor alleles, --min_maf 0.013 or 0.011; observed heterozygosity less than 0.5, --max_obs_het 0.5; genotyping rate >60% of individuals within populations, -r 0.6; and present in all of the populations, -p 8 or 9. Pairwise R2 values for each SNP pair were calculated with PLINK 1.90b (Purcell et al. 2007), and loci with values >0.8 were removed using the whitelist option in Populations to exclude SNPs with high linkage. In addition, SNP genotypes for three regions in R. tsusiophyllum (set 3), which were identified by population analyses (see Results), were called with the following restrictions: include all minor alleles, --min_maf 0.00; observed heterozygosity <0.5, --max_obs_het 0.5; genotyping rate >80%, 90% or 100% of individuals within regions, -r 0.8, 0.9 or 1.0; present in all of the regions, -p 3; and pairwise R2 values for each SNP pair < 0.8.

Phylogenetic and population genetic analyses

Phylogenetic relationships among individuals, including or excluding the outgroup species, were evaluated independently by constructing a neighbor-net based on p-distances using SplitTree4 (Bryant and Moulton 2004; Huson and Bryant 2006). In addition, phylogenetic relationships among populations, including or excluding outgroup species, were evaluated by constructing neighbor-joining trees based on DA distances (Nei et al. 1983). The significance of nodes was evaluated from bootstrap probabilities based on 1000 replicates of the neighbor-joining method using POPTREE2 (Takezaki et al. 2010).

The individual-based genetic structure in R. tsusiophyllum was inferred by the sparse non-negative matrix factorization method using sNMF (Frichot et al. 2014). Different numbers of ancestries (K) in the range 1–10 were tested, and ranked using the cross-entropy criterion.

Coalescent-based phylogenetic relationships were inferred by Bayesian framework with visualization of uncertainty in the majority-rule topology using SNAPP v. 2.4.7 implemented in BEAST v. 2.4.7 (Bouckaert et al. 2014; Bryant et al. 2012) from SNP set2 which include the outgroup species for rooting. Because SNAPP analysis requires high computational intensity when using a large number of SNPs and individuals, we extracted two samples with the highest genotyping rate from each population of R. tsusiophyllum and four samples of R. tschonoskii to reduce computational intensity (Stervander et al. 2015). We used default settings for priors including forward (u) and reverse (v) mutation rates. Two independent Markov chain Monte Carlo runs, each consisting of 1,000,000 iterations (saving every 1000 steps), were conducted to confirm stationarity. The convergence with effective sample size (ESS) for each parameter over 200 was checked using TRACER v. 1.6 (Drummond and Rambaut 2007) and two runs were combined after removing 10% as burn-in. A maximum-clade credibility tree was generated using TreeAnnotator v. 2.4.7, and topologies of trees and the consensus tree were visualized using DensiTree v. 2.2.6 (Bouckaert 2010). Relative migrations among regions based on Nm (Alcala et al. 2014) were estimated using divMigrate function of R package diveRsity v. 1.9.90 (Sundqvist et al. 2016) in R v.3.5.2 (R Development Core Team 2019). Asymmetry in migrations between regions was assessed by comparing 95% confidence intervals based on 1000 bootstrap replicates.

Demographic analysis

Historical changes in population size in each region (Chichibu, Hakone, and the Izu Islands) were reconstructed using Stairway plot v. 2 (Liu and Fu 2015). This is a model-flexible method for calculating changes of effective population size over time from expected composite likelihood based on site frequency spectrum (SFS). Folded SFS for each region was calculated using Arlequin v. 3.5 (Excoffier and Lischer 2010). As missing data affect folded SFS estimates, SFSs were calculated for the data with no missing values (-r 1.0 for Population implemented in Stacks). As such datasets include less SNPs than sets containing moderate levels of missing data (-r 0.8 and 0.9), we compared distributions of expected heterozygosity for the three types of datasets with different genotyping rates (-r 0.8, 0.9 and 1.0) and evaluated whether the small number of SNPs with no missing values in each region had sufficient information for the demographic analysis. In calculation, as minor alleles enveloping sequence errors cause deviations in estimates from true demographic changes, we conducted simulations both with and without singletons. Two hundred of bootstrap iterations were performed for each estimation, with mutation rates of 7.1 × 10−9 per site per generation, following Ossowski et al. (2010), to scale generations and 7.1 × 10−10 per site per year (assuming a generation time of 10 years for the species) to scale theta (θ) to the effective population size (Ne).

Climatic differentiation and association analyses

For evaluating differentiation of habitats or climatic niches of R. tsusiophyllum, distributions of individuals, with removal of duplications within each 30 arc-sec cell, were collected from herbarium records (TI and TNS) and field observations. Basic climatic variables at 30 arc-sec resolution were also obtained from WorldClim (http://www.worldclim.org) and terrestrial data at 1 arc-min resolution from the ETOPO1 global relief model (https://www.ngdc.noaa.gov/mgg/global/). Six temperature and precipitation variables (annual mean temperature, mean temperature of the warmest quarter, mean temperature of the coldest quarter, annual precipitation, precipitation of the warmest quarter, precipitation of the coldest quarter) and altitude were chosen for evaluating habitat differentiation likely to affect the species distribution. The climatic differentiation among populations was evaluated by principal component analysis (PCA) in R.

Evidence of natural selection was tested by analyzing associations between allele frequencies and environmental variables of populations by latent factor mixed modeling using LFMM v. 1.4 of R package LEA v. 1.99 in R (Frichot and François 2015, Frichot et al. 2013, 2015). LFMM uses each environmental variable as a fixed effect within each model, and includes a number of latent factors (K) as a covariate to take account of effects of population structure. The annual mean temperature, annual precipitation and precipitation of the warmest quarter were used as representatives of climatic variables and altitude, because other variables were highly correlated (R2 > 0.8) for analyzed locations. The number of latent factors used was K = 3, based on phylogenetic relationships among individuals and population genetic structure (see Results). Five independent simulations were conducted with 50,000 iterations after burn-in periods of 25,000. The z-scores calculated from multiple runs for an environmental variable were combined and re-adjusted to p values. Association between allele frequencies for each SNP and the environmental variable were evaluated by adjusted p values (q values) obtained using the Benjamini–Hochberg procedure (Benjamini and Hochberg 1995) and considered to be outliers if q < 0.00001. Mean allele frequencies and 95% confidence intervals for each candidate SNP and population were calculated by bootstrap resampling, based on 1000 replicates using R.

Results

Phylogenetic patterns and population genetic structure

The RAD-seq library yielded 197,461,828 reads with 1,645,515 on average per sample of Rhododendron tsusiophyllum (from 120 individuals), and 9,274,411 reads with 515,245 on average per sample for the outgroup species, R. tschonoskii (from 18 individuals) (Table S1). We detected 4407 SNPs, with a genotyping rate of 91.50% excluding the outgroup (set 1), and 3284 SNPs with a genotyping rate of 92.17% including the outgroup (set 2). Average coverage per locus (across all individuals) ranged from 4.23 to 387.89 (35.44 on average) excluding the outgroup, and from 4.90 to 357.84 (38.82 on average) including the outgroup. The number of polymorphic SNPs and expected heterozygosity within each population of R. tsusiophyllum (set 1) ranged from 389 to 1967 and from 0.025 to 0.115, respectively (Table 1). These variables were low for the two populations of the Izu Islands. The analysis of phylogenetic relationships by neighbor-net identified three groups, corresponding to populations of R. tsusiophyllum in Chichibu, Hakone and the Izu Islands, and individuals from each population formed a distinct sub-group (Supplementary Figs. 1 and 2). The grouping was the same according to neighbor-nets with and without the outgroup. The neighbor-net including the outgroup indicated that the outgroup was rooted between individuals in Chichibu and the other two regions (Hakone and the Izu Islands). The genetic structure of R. tsusiophyllum inferred by sNMF corroborated these findings, indicating that three was the optimal number of clusters (K) based on the cross-entropy criterion. With K = 3, each cluster dominated in one of the three regions, in accordance with the neighbor-net results (Fig. 2). When the number of clusters was increased to K = 5, the second most optimal number, populations of Chichibu and the Izu Islands separated into different clusters, respectively. Two runs of coalescent-based phylogenetic relationships among populations by SNAPP converged and combined data also identified the three regions (Fig. 3). The consensus tree showed high support for all nodes with the highest posterior probability (PP = 1.0), suggesting lineage sorting among three regions. The topology indicated that the outgroup was rooted between populations in Chichibu and the other two regions, corroborating the relationships identified by the neighbor-net analysis. High relative migration rates among regions were assessed between Chichibu and Hakone, and low relative migration rates were assessed between Chichibu and the Izu Islands (Supplementary Fig. 3). Asymmetric migrations were supported only the pair between Hakone and the Izu Islands.

Fig. 2: Individual-based genetic structure inferred by sNMF.
figure 2

Distribution of ancestry coefficients with numbers of clusters (K) ranging from 2 to 5.

Full size image
Fig. 3: Coalescence-based phylogenetic relationships among populations of Rhododendron tsusiophyllum including an outgroup species estimated using SNAPP.
figure 3

Thin gray lines indicate individual trees and black bold lines indicate the consensus tree. Posterior probabilities for all of nodes are 1.00.

Full size image

Population size and its historical changes in the three regions

The number of SNPs detected in Chichibu, Hakone and the Izu Islands were 4949, 7211 and 1403 SNPs, respectively, with a missing rate filter (-r 0.8); 3304, 4091 and 814 SNPs, respectively (-r 0.9); and 799, 822 and 227 SNPs, respectively (-r 1.0). Violin plots of expected heterozygosity distribution for each region under three mentioned missing rate filters showed similar pattern that loci with low heterozygosity were most frequent in Hakone despite large number of polymorphic SNPs (Supplementary Fig. 4). It was concordant with the pattern that loci with rare alleles (minor allele frequency < 0.05) were most frequent in Hakone with the folded SFS obtained for each region with no missing values (Supplementary Fig. 5). The patterns of effective population size changes estimated by stairway plots including or excluding singletons were similar within each region (Supplementary Fig. 6). Effective population sizes tended to be smaller in the Izu Islands than in the other two regions (Fig. 4). They also notably increased in Hakone and (to a lesser degree) Chichibu from 10,000 to 100,000 generations ago, while that of the Izu Islands slightly increased from 1000 to 5000 generations ago.

Fig. 4: Changes in effective population size in each region reconstructed by stairway plots based on site frequency spectrum.
figure 4

Lines show median and shaded areas show 95% highest posterior probability densities (95% HPDs). Numbers of generations are log-scaled.

Full size image

Habitat differentiation based on environmental factors and identification of outlier loci

Each environmental variable was extracted for 25 localities including analyzed populations. Three representative variables—altitude, annual mean temperature, and annual precipitation—ranged from 518 m to 2411 m, 0.8 °C to 14.2 °C, and 1678 mm to 2816 mm, respectively. The three regions were distinguished by the PCA based on seven environmental and geographical variables, especially localities of the Izu Islands, which were clearly distinguished from localities of Chichibu and Hakone along the first principal component (PC1, Fig. 5). The first axis explained 74.4% of the total variance, and was mainly related to the temperature and altitude variables.

Fig. 5: Differentiation of habitats of the analyzed populations and other observed individuals based on six climatic variables and altitude detected by principal component analysis (PCA).
figure 5

Gray areas show three regions defined from distribution of individuals (Fig. 1) and population genetic structure. The axis loading of each variable to the first two principal components (PC1 and PC2) are shown as a vector.

Full size image

Associations between allele frequency and PCA axes, as representatives of climatic condition were assessed by LFMM analysis. Histograms of adjusted p values were uniformly distributed, indicating that the number of latent factors (K = 3) adequately controlled effects of population structure. LFMM analysis identified 49 and 6 of candidate SNPs associated with PC1 and PC2 (q < 0.00001), and allele frequencies of candidates showed gradual changes with PCA axes (Supplementary Fig. 7).

Discussion

North to south divergence among the three regions

Analysis of phylogenetic relationships and genetic structure identified three genetic groups of Rhododendron tsusiophyllum corresponded to three geographic regions: Chichibu, Hakone and the Izu Islands (Fig. 3). Phylogenetic analysis also indicated that populations in Chichibu diverged from the other populations relatively soon after speciation between R. tsusiophyllum and R. tschonoskii (Yoichi et al. 2019). The distribution of R. tsusiophyllum clearly seems to have extended from cold northern to warm southern regions, considering marginal location of the Izu islands in the species range and the distribution and habitats of the related species. The results support the hypothesis that the distributions of Fossa Magna element plants, which are present on both the islands and the mainland, have been shaped by speciation from closely distributed relatives on the mainland and subsequent southern expansion into the islands (Bellemain and Ricklefs 2008; Takahashi 1971; Vaxevanidou et al. 2006).

We detected clear genetic divergence between the islands and mainland populations, as previously reported for other plants, such as between Lilium auratum var. auratum and L. auratum var. platyphyllum in the mainland and islands (Yamamoto et al. 2017). Further, the simple network obtained from neighbor-net analysis and robust relationship among regions identified by SNAPP support evidence of lineage sorting (Bryant et al. 2012; Foote and Morin 2016), in addition, few admixtures of ancestries detected by sNMF suggest there has been restricted migrations between regions (Fig. 2). Although migration rates were not estimated directly, relative migrations between Chichibu and Hakone were more frequent than between Hakone and the Izu Islands. Migrations between Hakone and the Izu Islands may be restricted by not only the geographic isolation but also environmental heterogeneity between the regions, as detected by environmental PCA (Fig. 5). Moreover, there is weaker genetic divergence between populations within the islands than between populations in Chichibu, and between populations of other plants found on the islands, such as Weigela coraeensis (Yamada and Maki 2012). The weak divergence has been attributed to the narrow distribution of R. tsusiophyllum on the Izu Islands, associated with altitudinal restrictions (few of the islands have altitudes exceeding 550 m) and other environmental limitations of the species in the region. The divergence of populations between Honshu and the Izu Islands seems to be easily explained by vicariance associated with the sea separating them. However it should be explained by montane isolation, because habitat of the species is restricted in upsides of mountains and there is also notable divergence of populations between Chichibu and Hakone within Honshu. The clear phylogenetic relationships and genetic structure also suggest that there has been restricted migrations among the three regions. Thus, the results clearly illustrate the promotion of genetic divergence by montane isolation, sometimes called ‘sky islands’ (DeChaine and Martin 2005; Warren et al. 2015) by lower altitude areas. Differences in population survival and other demographic changes associated with these geographical factors are discussed below.

Differences in population survival among the three regions

Demographic histories of the genetically diverged populations in the three identified regions seem to have substantially differed since the splits. The small Izu Islands also had smaller effective population size than two regions of the larger Honshu according to stairway plot (Fig. 4), as reported for birds in the Canary Islands (Stervander et al. 2015). This suggests that the low effective population size in the Izu Islands may have been maintained since the island colonization (Stuessy et al. 2014). In addition, the number of polymorphic SNPs and expected heterozygosity of the two populations of the Izu Islands are less than one-third that of the populations of Honshu (Table 1), the low genetic diversity, with individuals growing in isolated and restricted habitats, would also have experienced genetic drift. All these findings and theoretical considerations suggest that isolated and restricted habitats strongly influence to population survival on the small islands (Warren et al. 2015).

In contrast, effective population size in Hakone seems to have dramatically increased 10,000–100,000 generations ago. This may have been due to the numerous active or inactive volcanos (Table 1), as well as mountains uplifted by tectonic forces (Taira et al. 1998), providing more suitable habitats in this region. There are no volcanos in Chichibu, where R. tsusiophyllum is found at the highest altitude with similar climate to that of R. tschonoskii. Population size in the region was presumably related to the area of that climatic belt on high mountains. It is difficult to estimate absolute times of these increases, due to lack of knowledge of the species generation time. However, assuming a generation time of 10–20 years for the shrub species, they probably occurred hundreds of thousands to millions years ago, suggesting the population growth occurred during glacial-interglacial periods before the last glacial (Clark et al. 2009; Petit et al. 1999). Although mutation rates are also difficult to estimate for highly variable loci (Garot et al. 2019), both mutation rate and generation time are unlikely to have varied substantially among the regions. Thus, the inferred changes in population size in the regions can be compared (Koch et al. 2000; Ossowski et al. 2010). The differences in demographic histories strongly indicate that changes in suitable habitat area strongly influenced each population’s size under their isolated circumstances (Frankham 1997; Stuessy et al. 2014).

Possibility of local adaptation on isolated islands

We probed the possible adaptive evolution of populations in the Izu Islands by a combination of genome-wide scans and environmental analyses. PCA clearly revealed differences in habitats of the species in the three regions, especially between the Izu Islands and two regions of Honshu (Fig. 5). Essentially, R. tsusiophyllum has a wider altitudinal range (500–2400 m) and annual mean temperature (0.8–14.2 °C) than most plants globally (Chapin et al. 2002). Rhododendron tsusiophyllum seems to have shifted to a warmer niche on the Izu Islands from cold niche in Chichibu. In conclusion, some of the clear genetic divergence between the three regions presumably resulted from neutral processes associated with geographical isolation. However, adaptive evolution to habitats with different environmental characteristics (Izuno et al. 2017), particularly in temperature, also seems to have contributed to the genetic divergence between populations of Honshu and the Izu Islands. Volcanos may also be important environments for adaptive evolution, however volcanos are also located in Hakone with frequent eruptions, thus volcanos do not just characterize habitat of the Izu Islands. Rather, R. tsusiophyllum seems prefer simply wind-swept open sites in upsides of mountains formed by volcanic eruption or tectonic uplift.

Candidate SNPs with climatic-related changes in allele frequency, and thus possible association with adaptive evolution, were detected, especially we found 49 candidates against PC1 axis value. However this detection seems to contain false positives, because evolutionary process from Chichibu to the Izu Islands, based on all loci, and temperature gradient explained by PC1 are associated (Frichot et al. 2015; Supplementary Fig. 7). Caveats should be mentioned regarding these candidates’ potential links with adaptive evolution. Ideally the functions of candidates should be known or determined, but the shortness of our ddRAD-seq reads reduced the confidence of homology searches, in addition, no good reference data are available for the study species. Candidates have been detected by FST-based methods as outliers in other studies (Funk et al. 2016; Izuno et al. 2017), but these methods could not be applied here because of the variation in expected heterozygosity among populations (0.025–0.115), which affects estimation of FST (Hedrick 2005). Especially, fixation of allele, which frequently founded in the two populations of the Izu Islands, seems difficult to detect true candidates. Missing genotypes, which are frequent in ddRAD-seq datasets, also restricted calculation of allele frequencies. Thus, there is clearly scope for improving evaluation of the possible associations by addressing these restrictions (Rellstab et al. 2015).

Previous studies suggest that founder effects during initial colonization and subsequent genetic drift significantly hinder adaptive evolution of island populations (Hamabata et al. 2019). However, strong selective pressure due to environmental differences from source populations also play significant roles in adaptive evolution of populations under drift (Barton 1996), as shown by many examples of adaptive radiation in islands (Funk et al. 2016; Velo-Antón et al. 2012), and our analysis supports the opinion.