Keiko Kitamura
Kanji Namikawa
Yoshiaki Tsuda
Makoto Kobayashi4
Tetsuya Matsui- 1Hokkaido Research Centre, Forestry and Forest Products Research Institute, Sapporo, Japan
- 2Biological Laboratory, Hokkaido University of Education, Sapporo, Japan
- 3Sugadaira Montane Research Center, University of Tsukuba, Ueda, Japan
- 4Department of Education and Culture, Echigo-Matsunoyama Museum of Natural Science, Tokamachi, Japan
- 5Center of Biodiversity and Climate Change, Forestry and Forest Products Research Institute, Tsukuba, Japan
- 6Faculty of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Japan
Siebold’s beech, Fagus crenata, is widely distributed across the Japanese Archipelago and islands in Japan Sea. Similar to the northern limit of the geographical distribution of F. crenata on the mainland of Hokkaido, the northern limit of the distribution of F. crenata on islands in the Japan Sea is observed on Okushiri Island (ca 42°N). To understand the genetic relationships of F. crenata on Okushiri Island, we examined chloroplast (cp) DNA haplotypes and 11 nuclear microsatellite (SSR) loci among 1,838 individuals from 44 populations from Okushiri Island, mainland Hokkaido, and the northern part of the Tohoku region on Honshu Island. We identified 2 cpDNA haplotypes, which represent not only populations on the Japan Sea coast but also those on the Pacific coast and this suggested the Okushiri Island populations might not be formed by single colonization. Genetic diversity of the Okushiri Island populations of nuclear SSR was not lower than the mainland and the STRUCTURE analysis revealed the Okushiri Island individuals were admixed between Hokkaido and Tohoku clusters. Approximate Bayesian computation inferred that divergence between Tohoku and Hokkaido, and admixture between two populations which generated Okushiri populations occurred before the last glacial maximum (LGM), that is, 7,890 (95% hyper probability density (HPD): 3,420 – 9,910) and 3,870 (95% HPD: 431– 8,540) generations ago, respectively. These inferences were well supported by a geological history which suggested an isolation of Okushiri Island from Hokkaido started prior to the Middle Pleistocene. We discuss the possible persistence of F. crenata during the last glacial maximum on northern islands in the Japan Sea such as Okushiri Island.
Introduction
Plant populations on islands differ from mainland populations in terms of ecological and demographic processes such as colonization events, population persistence, and expansion because of their geographical isolation. From a genetic perspective, island populations are assumed to have less diversity than their mainland counterparts due to evolutionary processes such as founder effects, limited gene dispersal, and small population size (Barrett, 1996; Frankham, 1996; Frankham, 1997; Takayama et al., 2013; Stuessy et al., 2014; Takayama et al., 2015; Stuessy, 2020). For example, the initial colonizing event inevitably involves a genetic bottleneck which determines the loss of gene diversity in founder populations (Barrett, 1996). Therefore, newly founded populations carry less gene diversity than source populations; however, repeated gene flow from source populations can compensate for a lack of gene diversity and contribute to long-term population persistence and expansion (Alsos et al., 2015). Theoretically, long-distance gene flow is less frequent than short-distance gene flow, so that spatial isolation and small population size may lead to an erosion of genetic variation and increased interpopulation differentiation (Young et al., 1996). In general, geographical isolation is a major contributor to shallow gene diversity (DeJoode and Wendel, 1992; Chiang et al., 2006; Takayama et al., 2013; Alsos et al., 2015).
On the contrary, several studies of terrestrial plant species revealed that island populations have substantial genetic diversity in spite of their isolated geographic locations which might restrict frequent gene flow (Chiang et al., 2006; Rosas Escobar et al., 2011; Désamoré et al., 2012; García-Verdugo et al., 2013; García-Verdugo et al., 2015). This is thought to be due to the fact that the genetic diversity of island populations is not only related to founder effects but also by different geographic and biological conditions and contexts (Stuessy et al., 2014), such as adaptation and maladaptation to the new environment (St Clair and Howe, 2007; Kuparinen et al., 2010) and chronic changes of gene immigration (Elleouet and Aitken, 2019; Kitamura and Nakanishi, 2021).
Another interesting issue regarding island populations is the colonization history during the species range expansion (Taberlet et al., 1998; Magri et al., 2006; Alsos et al., 2015). Compared to mainland populations, the routes and frequency of dispersal to islands are restricted due to geographic isolation, and long-distance gene dispersal represents the initial stage in the colonization process.
Fagus crenata Blume (2n = 24) is monoecious and allogamous, with barochorous and synzoochorous by rodents and birds (Miguchi, 1996). The geographic range of F. crenata extends across the Japanese Archipelago, from southern Hokkaido to Kyushu (Peters, 1997), with sharp geographic clines in genetic diversity (Tomaru et al., 1997; Hiraoka and Tomaru, 2009) and morphology (Hagiwara, 1977; Ishii et al., 2018) observed. Our previous study of genetic structure at the northernmost distribution range on mainland Hokkaido discovered a decline in gene diversity toward the northern edge of populations (Kitamura et al., 2015). In addition to the mainland distribution, natural populations of F. crenata are observed on off-shore islands in the Japan Sea, with Okushiri Island (ca. 42°N) representing the northernmost distribution of F. crenata among such islands (Tatewaki, 1948). The island is reported to have been separated from the mainland in the Middle Pleistocene and never became reattached to the mainland thereafter (Ohshima, 1980). The northward expansion of F. crenata after the Last Glacial Maximum (LGM) started about 6,000 years BP on the mainland (Kito and Takimoto, 1999); however, the origins of the northernmost island population of F. crenata is still unknown. Long-distance vs. short-distance gene dispersals are driving forces to develop genetic diversity within populations. As tracers of gene flow, neutral molecular markers provide useful information, although care needs to be taken in separating historical connectivity from current gene flow. A combination of maternally inherited chloroplast (cp) genome for angiosperms and bi-parental inherited nuclear genome provides relevant information such as the colonization history because comparisons of uniparentally (cp) vs. biparentally (nuclear) inherited loci can reveal differences in dispersal between seeds and/or pollen (Petit and Excoffier, 2009). Moreover, recent advances in genetic-based population demographic inference make it possible to investigate a colonization history of species from different ancestral lineages, considering time scales (Tsuda et al., 2015).
The aims of this study are: 1) to evaluate the genetic structure of F. crenata on Okushiri Island by 44 populations from northern distribution using nuclear and cpDNA variation, 2) to infer the species’ population demography on the island, and 3) to discuss the colonization and northern expansion dynamics of F. crenata, considering a possibility of northern persistence of the species during the LGM.
Materials and methods
Study sites
Okushiri Island is located 16 km west of the mainland of Hokkaido in the island-mainland northernmost distribution area of F. crenata (Figure 1). This Island is a small island of 143 square km, where forest covers 80% of the island, yet the F. crenata forest covers approximately 60% of the forested area (Namikawa et al., 2021), and its stand volume is around 60 cubic meters per hectare (Tatewaki, 1948).
Figure 1 Locations of Fagus crenata populations examined in this study (Supplementary Table 1) Populations 18, 20, and 31 were referred in the Discussion.
In terms of the vegetation classification, the beech forests on Okushiri Island showed diverse floristic features from subarctic to warm temperate elements which is partly due to the warm current (Namikawa et al., 2021). The major vegetation type of beech forest on the island includes species that frequently appear in beech forests on the Japan Sea coast of Hokkaido and northern Tohoku. Another vegetation type found on Okushiri Island is associated with the southern part of the Pacific coast of Hokkaido, which has a short snow-cover period. Thus, despite the small size of the island, regional differences in species composition seen in mainland Hokkaido and Tohoku were condensed in the beech forests on Okushiri Island (Namikawa et al., 2021).
We selected 17 populations from across the F. crenata forests on the island, 16 populations from the mainland of Hokkaido, and 11 from the northern Tohoku region of Honshu as study sites, covering most of the species’ distribution in these areas (Figure 1, Supplementary Table 1). To cover most of the natural distribution, study populations on Okushiri Island were selected from both east- and west-facing slopes of the island, as well as coastal and inland habitats. Populations from mainland Hokkaido and the northern Tohoku region were selected across the contrasting climate conditions of the Pacific coast that shows the dry and low temperature winter, and the Japan Sea coast, which experiences heavy snow fall in winter.
DNA extraction, chloroplast DNA SNPs and microsatellite genotyping
Mature trees of diameter at breast height (DBH) > 20 cm (11 to 65 trees per population) were chosen at random, and leaves were collected from a total of 1,838 individuals from 44 populations (Supplementary Table 1). Collected leaves were kept in a cooling box with ice until they could be stored in the refrigerator or freezer in the laboratory. Total DNA was extracted with the DNeasy Plant Mini Kit (Qiagen K. K., Tokyo, Japan) from 50 mg of leaf sample that was ground into a powder using liquid nitrogen and a Multi-Beads Shocker cell disruptor (Yasui Kikai Co. Ltd., Osaka, Japan).
It has been revealed that there are two chloroplast (cp)DNA haplotypes in these regions, that is, A (GenBank: AB046492.1) and B (GenBank: AB046493.1) (Fujii et al., 2002). We determined the cpDNA haplotypes of each population by SNPs according to the method described by Katai et al. (2014); Takahashi et al. (2022). The first polymerase chain reaction (PCR) was performed to amplify trnK fragment by 0.2 µM final concentration of each forward (5’-TTATTCTTAGCGGATCGGTCCA-3’) and reverse (5’-CCGTGCTTGCATCTTTCATTG-3’) primer using Multiplex PCR Kit (Qiagen K. K., Tokyo, Japan) containing 10 to 30 ng of genomic DNA in a final volume of 6 µl. Amplification condition was 95°C for 15 min., 35 cycles of [94°C for 30 sec., 56°C for 90 sec., 72°C for 2 min.], 72°C for 10 min, then hold at 4°C. Amplicons were treated with 0.8 U of ExoI (Takara Bio Inc., Shiga, Japan) and 2 U of Shrimp Alkaline Phosphatase (SAP, Takara Bio Inc., Shiga, Japan) with ExoI buffer (Takara Bio Inc., Shiga, Japan) to remove unincorporated primers and dNTPs for 1 hour at 35°C followed by 15 min. at 75°C for inactivation of enzyme, then hold at 4°C. The single-base primer extension was performed using SNaPshot Multiplex Kit (Thermo Fisher Scientific K. K., Tokyo, Japan) and trnK1754 primer (5’- T14CTAGCATTTGACTCCGCACCACTGAAG -3’). The total volume of the reaction mix was 7 µl which contained 1 µl of SNaPshot Master Mix, 0.7 µl of 2 µM primer, 1 µl of BigDye Terminator 5 x Sequencing buffer (Thermo Fisher Scientific K. K., Tokyo, Japan), and 2 µl of the first PCR product. The extension was performed under condition of 25 cycles of [96°C for 10 sec., 50°C for 5 sec., 60°C for 30 sec.], then hold at 4°C. SNaPshot extension reaction was treated with 1 U of SAP (Takara Bio Inc., Shiga, Japan) to remove unincorporated primer and ddNTPs for 1 hour at 35°C followed by 15 min. at 75°C for inactivation of enzyme, then hold at 4°C. The fragments were detected with ABI PRISM 3130xl Genetic Analyzer (Thermo Fisher Scientific K. K., Tokyo, Japan) using POP-7 polymer and 36 cm capillary. An aliquot of SNaPshot reaction (1 µl) was mixed with 9.8 µl of Hi-Di Formamide (Thermo Fisher Scientific K. K., Tokyo, Japan) and 0.2 µl of GeneScan 120 LIZ dye Size Standard (Thermo Fisher Scientific K. K., Tokyo, Japan). SNPs were determined by GENESCAN for Windows (Thermo Fisher Scientific K. K., Tokyo, Japan). We first analyzed four individuals from each population to determine cpDNA haplotypes. When four individuals from the same population showed the same haplotypes, the cpDNA haplotype of the population was determined to be monomorphic. When different haplotypes were detected from four individuals from the same population, we further analyzed 12 to 22 individuals from that population to determine the cpDNA haplotype ratio of the population. The cpDNA haplotype are fixed in most of the populations of this species, and we analyzed the limited number of individuals to detect at least 6% of polymorphism.
Eleven nuclear microsatellite (SSR) primer pairs were employed to identify genotypes of all individuals for the following loci: mfc2, mfc12 (Tanaka et al., 1999), FS1-03, FS4-46 (Pastorelli et al., 2003), sfc7, sfc18, sfc36, sfc378, sfc1063, sfc1105, and sfc1143 (Asuka et al., 2004), according to the method described in (Kitamura et al., 2015). PCR was performed with the Multiplex PCR Kit (Qiagen K. K., Tokyo, Japan) containing 10 to 30 ng of genomic DNA in a final volume of 10 µl. Amplification condition was 95°C for 15 min., 35 cycles of [94°C for 30 sec., 57°C for 90 sec., 72°C for 1 min.], 60°C for 30 min, then hold at 15°C. PCR products (1 µl) was mixed with 10 µl Hi-Di Formamide (Thermo Fisher Scientific K. K., Tokyo, Japan) and 0.15 µl of GeneScan 600 LIZ dye Size Standard v2.0 (Thermo Fisher Scientific K. K., Tokyo, Japan) for estimating DNA fragment sizes. The length of amplified fragments was analyzed using the ABI PRISM 3130xl Genetic Analyzer using POP-7 polymer, 36 cm capillary, and GENESCAN for Windows (Thermo Fisher Scientific K. K., Tokyo, Japan).
The number of individuals used for analyses are shown in Supplementary Table 1.
Data analyses
The genetic diversity parameters within each population were evaluated by determining the expected heterozygosity (HE) (Nei 1987), allelic richness (RS) (El Mousadik and Petit, 1996), and the fixation index (FIS) (Wright 1965) using the FSTAT 2.9.3 computer program (hereafter, FSTAT) (Goudet, 2001). Significant deviations in FIS values from 0 were estimated by 95% confidence intervals (CI) obtained through 999 permutations of bootstrapping using GENODIVE (Meirmans and Van Tienderen, 2004). The degree of genetic differentiation among the populations was evaluated by Nei’s G’ST. The number of alleles and the effective number of alleles within a region were calculated by GENODIVE. Frequency of null allele were calculated by CERVUS (Marshall et al., 1998).
To determine the coancestry composition among populations and admixtures, we employed multilocus model-based cluster analysis using STRUCTURE 2.3.4 (hereafter, STRUCTURE) (Pritchard et al., 2000). All of the runs consisted of 30,000 Markov chain Monte Carlo (MCMC) generations, after a burn-in period of 70,000 iterations. We confirmed that these settings for the length of the MCMC generations and burn-in period are sufficient to obtain valid genetic structures by comparing the results from initial runs with those involving longer MCMC generations. This analysis was based on the LOCPRIOR model described by Hubisz et al. (2009), an admixture model, and the correlated allele frequencies model (hereafter, the F-model) described by Falush et al. (2003). Ten runs were performed for each value of K, ranging from 1 to 10; as we confirmed that the genetic structure among population could be detected at K < 10 in the pilot run. CLUMPAK server (Kopelman et al., 2015) was employed to evaluate the probability of the data [LnP (D)] for each K and to calculate ΔK according to the method described by Evanno et al. (2005), and to evaluate multimodality among runs and major clustering patterns at each K.
In addition, principal coordinate analysis (PCoA) of 44 populations based on the allele frequencies of 11 SSR loci, as well as Mantel test for isolation by distance between the pairwise genetic distance and geographic distances were conducted by GENODIVE (Meirmans and Van Tienderen, 2004). Phylogenetic networks among populations on Okushiri Island was obtained from neighbor-net split network method based on pairwise genetic distance using SplitsTree4 (Huson and Bryant, 2006).
Inference of past population demographic history using the approximate bayesian computation
To infer past demographic history of colonization, we employed the Approximate Bayesian Computation (ABC) approach and used software DIYABC v2.1 (Cornuet et al., 2008; Cornuet et al., 2014). To make the scenarios in the ABC analysis simple, we focused on three representative populations: (1) the mainland of Hokkaido, (2) Okushiri Island, and (3) Tohoku region. These representative populations were made as follows. According to the results of the above-mentioned STRUCTURE analysis, three clusters were detected that corresponded to these three areas at K = 3. Then the 50 individuals which showed the highest ancestry within each cluster were selected, and named to Pop1 (the mainland of Hokkaido), Pop2 (Okushiri Island) and Pop3 (Tohoku region). These three representative populations, that is a total of 150 individuals, were examined six simple population demography scenarios (Figure 2).
Figure 2 Six simple population demography scenarios used for Approximate Bayesian Computation (ABC).
Since DIYABC requires a population which can be traced back to an ancestral population, we set Pop3 to be the ancestor because the Tohoku region showed higher genetic diversity than the mainland of Hokkaido (Takahashi et al., 1994; Hiraoka and Tomaru, 2009) (see also, Table 1). This treatment is reasonable according to Tsuda et al. (2015). In these scenarios, t1 to t6 represents the time scale, measured by generation time, N1 to N3 represents the effective population size of the corresponding populations (Pop1, 2 and 3) and ra to rd represents the ratio of admixture. Although the DIYABC basically does not take gene flow after the split into the scenario, a recent study by Chapuis et al. (2020) modified the DIYABC to allow symmetrical admixture for two population data and we extended this approach to make it more flexible for multiple population data, making it possible to put directional admixture (e.g., from Pop1 to 2, from Pop2 to 1) in the scenarios. Thus, this directional admixture can be analogues to gene flow. The inferred value of admixture parameter (r#) from Pop1 to Pop2 could be considered as the relative amount of gene flow at t# and 1-r# was the value of relative amount of gene flow within a parental population for the admixture (gene flow).
Table 1 Genetic parameters for each region.
Scenario 1: Hierarchical split model I: Pop2 was merged to Pop3 at t3 and Pop1 was merged to Pop3 at t6, assuming populations on the Okushiri Island were split from Tohoku region. The directional admixture from Pop2 to Pop3 (ra), Pop3 to Pop2 (rb), from Pop1 to Pop3 (rc), and from Pop3 to Pop1 (rd) were set at t1, t2, t4, and t5, respectively. As we did not set any orders of these time scale for the admixtures, the time scale order relied on the inference (e.g., t2 could be earlier than t1) and it was the same in the following Scenarios (2-6).
In Scenarios 3, 4, and 5, the unidirectional admixture was adopted whilst populations on Okushiri Island were assumed to be too small to influence mainland populations by emigrant genes.
Scenario 2: Hierarchical split model II: Pop2 was merged to Pop1 at t3 and Pop1 was merged to Pop3 at t6, assuming populations on Okushiri Island were split from the mainland of Hokkaido. Similar to Scenario 1, the directional admixtures between populations after the splitting were set.
Scenario 3: Simple split model with admixture: all 3 populations diverged at the same time at t3, allowing admixture from Pop1 to Pop2 (ra) and Pop3 to Pop2 (rb) at t1 and t2, respectively.
Scenario 4: Hierarchical split model III: The hierarchical split pattern was similar to Scenario 2. However, we assumed only one directional admixture from Pop1 to Pop2 (ra) at t1.
Scenario 5: Isolation with admixture model I: Pop1 was split from Pop3 at t4 and Pop2 was created by admixtures from Pop3 to Pop2 (rc), and from Pop1 to Pop2 (1-rc) at t3. The secondary directional admixture from Pop1 to Pop2 (ra) and from Pop3 to Pop2 were set at t1 and t2, respectively.
Scenario 6: Isolation with admixture model II: Although this scenario was similar to Scenario 5, we assumed all possible combinations of secondary directional admixture from Pop2 to Pop3 (ra), from Pop3 to Pop2 (rb), from Pop1 to Pop2 (rc), from Pop2 to Pop1 (rd), at t1, 2, 3 and 4, respectively.
The prior distributions of these parameters were shown in Supplementary Table 2. We employed the higher mutation rate from the generalized stepwise mutation model (GSM) (Estoup et al., 2002) and the lower rate for single nucleotide indels (SNI) as mutation models of SSRs. To summarize the observed and simulated data, the mean values for expected heterozygosity (HE), the number of alleles (A), and the variance of allele size were used as summary statistics for each of the three populations. HE, A, variance of allele size, classification index (analogues to genotype likelihood) and FST were the summary statistics for each of the population pairs. One million simulations were run for each scenario. After all the simulations had been run, the most likely scenario was determined by comparing the posterior probabilities using the logistic regression method. The goodness of fit of the scenario was assessed by the option ‘model checking’ with PCA in DIYABC, which measures the discrepancy between data sets simulated with the prior distributions of parameters and the observed data as well as data sets from the posterior predictive distribution in the scenario.
Results
cpDNA polymorphism
Two haplotypes were distinguished; namely, A (GenBank: AB046492.1) and B (GenBank: AB046493.1) (Fujii et al., 2002), among the populations studied (Supplementary Tables 1, 3, Figure 3). The population was determined to be fixed to either haplotype, when the first analysis of 4 individuals showed the same haplotypes. An additional 12 to 22 individuals were analyzed at populations with mixed haplotypes. Most of the populations were fixed to either haplotype; 16 populations to haplotype A and 23 populations to haplotype B. Most of the populations in Hokkaido were fixed to haplotype A, but most of those on Okushiri Island and in Tohoku were fixed to haplotype B. The other 5 populations included both haplotypes, that is, 1 population in Hokkaido, 3 populations on Okushiri Island, and 1 population in Tohoku. However, the ratio of haplotype A/B differed among populations (Supplementary Table 1, Figure 3). The lower representative haplotype in each region showed specific localization; namely, haplotype B in the southeastern part of Hokkaido, haplotype A in the southeastern part of Okushiri Island, and haplotype A in the northwestern part of Tohoku (Figure 3).
Figure 3 Pie diagrams of cpDNA haplotypes. Red and yellow represent haplotype A and B, respectively.
Nuclear gene diversity by SSR
The total heterozygosity (HT) for all 44 populations was 0.783, and the genetic differentiation for all populations (G’ST) was 0.026. The HE values of Hokkaido, Okushiri Island, and Tohoku were 0.776 (values for each population ranges from 0.711 to 0.792), 0.778 (from 0.726 to 0.792), and 0.797 (from 0.731 to 0.800), respectively (Table 1, Supplementary Table 1). The effective number of alleles of Hokkaido, Okushiri Island, and Tohoku were 5.639 (each population value ranges from 4.811 to 6.338), 5.332 (from 4.576 to 6.377), and 5.795 (from 4.848 to 6.809), respectively (Table 1, Supplementary Table 1). The G’ST for Hokkaido, Okushiri Island, and Tohoku were 0.016, 0.023, and 0.036, respectively (Table 1). FIS of Hokkaido and Okushiri Island did not show significant deviation from zero, while that of Tohoku was positively deviated from zero (Table 1). The positively significant FIS was observed at one population from Hokkaido and Okushiri Island, and 4 populations from Tohoku; while negatively significant FIS was observed at one population from Okushiri Island (Supplementary Table 1). There were no significant regional differences in genetic diversity parameters, such as HE, the effective number of alleles, G’ST, and FIS, among Hokkaido, Okushiri Island, and the Tohoku region (Table 1, Supplementary Table 1, Figure 4). The allelic richness ranged from 7.153 to 8.007, from 6.677 to 8.045, and from 7.056 to 8.325, for Hokkaido, Okushiri Island, and Tohoku, respectively (Supplementary Table 1). The number of alleles for Hokkaido, Okushiri Island, and Tohoku were 22.273 (each population value ranged from 7.182 to 13.364), 22.273 (from 7.818 to 13.636), and 23.455 (from 11.091 to 14.364), respectively (Table 1, Supplementary Table 1). High frequency of null allele (0.278) along with a significant deviation from Hardy-Weinberg equilibrium was estimated in mfc12 at population 8 (Supplementary Table 4), however, we believe that it has little influence for the analysis.
Figure 4 Genetic diversity parameters for 3 regions. (A) HE, (B) effective number of alleles, (C) G’ST, and (D) FIS.
Coancestry composition among population
Although the results from STRUCTURE analysis revealed that ΔK was the highest when K = 2 (Supplementary Figure 1), the Ln P(D) increased with K. However, variances of Ln P(D) values among runs were larger after K = 4. As regards ΔK, even though a hierarchical population scenario is not likely, the highest value can be detected at the uppermost hierarchical level of the population structure (Evanno et al., 2005). Similarly, whilst the highest ΔK is often the case of K = 2, further population structure could be revealed in K > 2 (Pritchard et al., 2000). In this study, the second highest ΔK was obtained at K = 4, whereas the bar plots of K = 3 and 4 show similar population structures. Therefore, K = 3 could be still a meaningful clustering with the variance of Ln P(D) being smaller than K > 4.
We compiled the individual coancestry to each population and drew pie diagrams on geographical map (Figures 5A, B). When K = 2 (Figure 5A), the F values (analogue to FST between assumed common ancestral population and cluster) were 0.0156 and 0.0275 for Cluster I and II, respectively. The Cluster I was dominant in Tohoku and Cluster II in both Hokkaido and Okushiri Island. Cluster II was more dominant in Okushiri Island than Hokkaido.
Figure 5 Pie diagrams of STRUCTURE results. The different colors indicate the different ancestral clusters. (A) K = 2, red and green indicate Cluster I and II, respectively. (B) K = 3, red, green, and blue indicate Cluster I, II, and III, respectively. (C) Relationships among the 3 clusters in K = 3 based on FST.
When K = 3 (Figure 5B), the F values were 0.0125, 0.0122 and 0.0141 for Cluster I, II, and III, respectively. Three different clusters dominated the three different regions; that is, Custer I, II, and III for Tohoku, Hokkaido, and Okushiri Island, respectively. Unrooted neighbor-joining cladograms among the 3 clusters based on FST genetic distance matrix showed a closer relationship between Cluster II and III than between Cluster II or III and Cluster I (Figure 5C).
Principle coordinate analysis, isolation-by-distance, and within-island phylogenetic network
The results of principle coordinate analysis (PCoA) based on the correlation matrix are shown in Figure 6. The first coordinate substantially split the populations in Tohoku from those in Hokkaido and on Okushiri Island. The second coordinate weakly divided populations on Okushiri Island from those in Hokkaido. Also, the second coordinate isolated one population (site 16) located at the tip of the southeastern peninsula of Hokkaido. The cpDNA haplotype differences were reflected in the PCoA diagram (Figure 6); populations fixed with haplotype A (red symbols in Figure 6) were plotted at the center and those with haplotype B (black) scattered on the periphery of the coordinate plane.
Figure 6 Principle coordinate analysis based on the correlation matrix. The first and second axes represent the first and second main components, respectively. Contribution rates of each axis are given in parentheses. Crosses, filled circles, and open circles indicate Hokkaido, Okushiri Island, and Tohoku, respectively. Red, black, and grey indicate cpDNA haplotype A, B, and mixed A and B populations, respectively.
Significant isolation-by-distance was detected between study sites (Mantel test by 1,000 permutations; R2 = 0.198, p = 0.001) (Figure 7).
Figure 7 Isolation-by-distance between populations. Genetic distance matrix is based on FST/(1 -FST). Mantel test by 1,000 permutation detected positive correlation between genetic and geographic distances (y = 0.015 + 9E-05x, R2 = 0.198, p = 0.001).
A phylogenetic network among populations within Okushiri Island was revealed by neighbor-net split tree (Figure 8). The tree distinguished two geographic localities within the island; southern and northwestern populations. The existence of cpDNA haplotype A on the island was associated with phylogenetic network.
Figure 8 Phylogenetic network among populations within Okushiri Island based on pairwise genetic distance. Population numbers are identical to Figure 1 and Supplementary Table 1. Two distinguished localities, south and northwest, are indicated by blue dotted lines. Population in red letters include cpDNA haplotype A.
Inference of past population demographic history using the approximate bayesian computation
In DIYABC, the highest posterior probability was found in Scenario 5 (Isolation with admixture model I) (Table 2, Supplementary Figure 2). The value (0.5388, 95% CI: 0.5305 - 0.5471) was much higher than other scenarios and the 95% CI of Scenario 5 was not overlapped with other scenarios. For Scenario 5, the median values of effective population size of N1 (Pop1), N2 (Pop2), and N3 (Pop 3) were 7,820 (95% hyper probability density (HPD): 4,600-9,640), 5,680 (95% HPD; 2,150-8,640) and 9,180 (95% HPD; 6,450-9,970), respectively (Supplementary Table 5). The median values of t1, t2, t3, and t4 were 1,780 (95%HPD: 65.8 -7,500), 4,380 (95%HPD: 998-8,730), 3,870 (95%HPD: 431-8,540) and 7,890 (95%HPD: 3,420-9,910) (Supplementary Table 5). The median values of ra, rb and rc were 0.615 (95%HPD: 0.0585-0.983), 0.615 (95%HPD: 0.0402-0.982) and 0.505 (95%HPD: 0.0272-0.975), respectively (Supplementary Table 5). The median values of the mean mutation rate of SSR and SNI were 8.68 ×10-4 (95% HPD; 4.98×10-4-9.99×10-4) and 8.69 ×10-7 (95% HPD; 1.65×10-8-8.50×10-6), respectively (Supplementary Table 5). The median value of mean P, the parameter of the geometric distribution to generate multiple stepwise mutations was 0.244 (95%HPD; 0.117-0.300) (Supplementary Table 5). Only one of the 27 summary statistics of the simulated data was significantly different from the observed data and the PCA yielded a large cloud of data from the prior and observed datasets, centered around a small cluster from the posterior predictive distribution (Supplementary Figure 3), suggesting a good fit of the scenario together with high posterior probability.
Table 2 Posterior probability of each scenario and its 95% confidence interval based on the logistic estimate by DIYABC.
Discussion
Comparable amount of genetic diversity on Okushiri Island to Hokkaido and Tohoku region
The calculated values of HT of nuclear SSRs of this study (0.711-0.800, mean value: 0.736) were lower than that for the whole geographical range of F. crenata (0.793-0.873, mean value: 0.839) (Hiraoka and Tomaru, 2009). This is because the northern F. crenata populations show lower genetic diversity than the geographically central populations (Hiraoka and Tomaru, 2009; Kitamura et al., 2015). Our previous study on the northernmost F. crenata populations on mainland Hokkaido revealed a decline in genetic diversity at the northern edge of the species distribution (Kitamura et al., 2015). Although Okushiri Island represents the northernmost distribution of F. crenata on islands in the Japan Sea, our data did not show any decline in gene diversity (Table 1; Figure 4). The northernmost populations on mainland Hokkaido were relatively new founder populations at the northern front of its expansion (Kitamura et al., 2015; Kitamura and Nakanishi, 2021). On the other hand, the high genetic diversity on Okushiri Island indicated that the F. crenata populations on that island were not new founders, but might be populations that have persisted for a long time. This was also supported by the fact that F. crenata on Okushiri Island were well differentiated with a similar level of GST as Hokkaido and the Tohoku region (Figure 4). Thus, although the island is surrounded by ocean, which is a strong practical barrier against gene flow, the effective population size of F. crenata remained large enough to compensate for the isolation from the mainland population.
In general, the loss of genetic diversity in island populations can be explained by founder effect during population foundation, limited gene flow and genetic drift due to geographical isolation and small population size (Barrett, 1996; Frankham, 1997). However, our nuclear SSR results revealed that populations on Okushiri Island showed comparable amount of genetic diversity when we focused on Hokkaido and the Tohoku region (Figures 4A–D). As a result, the common assumption that island populations have lower genetic diversity than mainland populations (Barrett, 1996; Frankham, 1997) was not applicable to the F. crenata populations on Okushiri Island. García-Verdugo et al. (2015) have reported similar results that the nuclear markers do not tend to show lower genetic variation in island populations than in mainland ones.
Relevant amount of genetic diversity on Okushiri Island is also due to the fact that this island is covered with deciduous-broadleaved forest including a substantial population of F. crenata (Tatewaki, 1948; Namikawa et al., 2021). Outcrossing by wind-pollination of the species might compensate for the topographic gene flow barriers such as mountain ridges and valleys (Austerlitz et al., 2004; Hu et al., 2008). When genetic connectivity was secured, a large population size might prevent the loss of gene diversity on the island. Also, it is probable that ancestral alleles were retained within the population (Su et al., 2018) as the species has long life span of 250 years on average in Hokkaido (Kitamura et al., 2007) and the oldest tree on Okushiri Island was confirmed to be 374 years by annual ring counts (T. Matsui, M. Kobayashi, and K. Kitamura, unpublished data). Another supportive evidence of persistence of F. crenata on Okushiri Island is that the number of private alleles of SSR on Okushiri Island was 16 which was comparative to 14 in Hokkaido. Therefore, the levels of genetic diversity among F. crenata on Okushiri Island might be attributable to the large population size, which allows a significant degree of gene flow, and long demographic history, rather than new founder populations from a few individuals. Indeed, although the ABC-based demographic inference showed that the effective population size of Okushiri Island population was smaller than Hokkaido and Tohoku region, the inferred median value was 5,680 (Supplementary Table 5), which could be a criterion at least to maintain genetic diversity.
CpDNA genetic differentiation among F. crenata populations in Okushiri Island, Hokkaido, and Tohoku region
Our results from cpDNA haplotypes clearly separated Hokkaido from the Tohoku region (Supplementary Table 1, Figure 3). A primary difference between F. crenata in Hokkaido and the Tohoku region was supported by phytosociology, which categorized the different vegetation types in F. crenata forests in Hokkaido and the Japan Sea side of Honshu (Hukusima et al., 1984). Furthermore, the vegetational component of several F. crenata populations on the Pacific Ocean side of Hokkaido and the northwestern peninsula of Tohoku were classified into the same category as that on Okushiri Island (Namikawa et al., 2021), which was compatible with the existence of the same haplotype B on Okushiri Island and the Pacific Ocean side of southern Hokkaido (Figure 3). The areas are also characterized by short snow-cover period (Namikawa et al., 2021). These vegetational data might afford circumstantial evidence to support the existence of relict F. crenata populations on Okushiri Island. In addition, Takahashi et al. (2022) revealed that the leaf area of F. crenata with haplotype B was significantly smaller than that of haplotype A individuals, which might be adapted to the dry winter environment, namely short snow-cover period.
Our cpDNA haplotype data suggested that a primary differentiation between Hokkaido and the Tohoku region. Fagus crenata on Okushiri Island was dominated by the Tohoku lineage yet showed partial establishment of the Hokkaido lineage in the southeastern populations. In addition, nuclear SSR showed coancestry lineages of Okushiri Island belong to Hokkaido lineages at K = 2 of the STRUCTURE analysis. The pollen immigration occurred at the northernmost F. crenata population at least 12 km in distance (Kitamura and Nakanishi, 2021). This suggested that gene flow by pollen might occur between Okushiri Island and Hokkaido, where the nearest distance between the two regions is 16 km.
Organelle DNA, such as cpDNA, is maternally inherited and dispersed only by seeds (Mogensen, 1996). Thus, the spatial structure of cpDNA haplotypes might indicate the maternal genetic lineages. Previous studies of organelle DNA diversity of F. crenata revealed single or closely related lineages in the northern distribution range (Koike et al., 1998; Tomaru et al., 1998; Fujii et al., 2002; Okaura and Harada, 2002). We focused on local populations in the northern distribution range in this study and observed two haplotypes which can be distinguished by a single SNP difference. These different haplotypes appear to indicate the existence of at least two ancestral origins in the northern distribution range (Hewitt, 1999). Haplotype B was shared among northwestern Okushiri Island and the majority of the Tohoku region, while haplotype A was shared among southeastern Okushiri Island and Hokkaido. Hence, the northwestern populations on Okushiri Island might be more closely related to those in the Tohoku region in this study, while those in southeastern Okushiri Island are closely related to the populations in Hokkaido. This spatial structuring of cpDNA haplotypes might reflect the different expansion routes by seed into Okushiri Island. Like the former studies which characterized refugia of this species by cpDNA haplotypes (Okaura and Harada, 2002; Katai et al., 2011; Katai et al., 2017), different haplotypes on the island might indicate the existence of refugia.
Ooi (2016) resolved the species distribution history by an extensive pollen record. According to this, the most recent increase of Fagus pollen in Hokkaido observed after period 5 ka (thousand years ago), with its dominance expanding northward. This is in concordance of the species increase along the Japan Sea coast of Tohoku reflecting the increase of snowfall in this area. However, during the late-glacial (from 14 ka to 12 ka), Fagus pollen occurred in southern Hokkaido (Ooi, 2016), an area that corresponded to populations 5, and 11 to 16, which shared haplotype B. Moreover, Hoshino (2001) revealed that F. crenata have occurred continuously on Okushiri Island from 6,650 ± 120 years BP to the present. Being so, the species might have existed prior to the most recent northward increase along the Japan Sea coast. Given that the haplotype B is of an earlier prosperous lineage, this lineage, namely relict of F. crenata, has been sustained within Okushiri Island. Then, the other lineage, which should be represented by haplotype A, increased dominance from the relict source on the island (Ooi, 2016) or immigrated from the Japan Sea coast of Hokkaido. Possible seed dispersal agents of F. crenata were spotted nutcrackers (Nucifraga caryocatactes) (Vander Wall and Balda, 1977; Kobayashi and Watanabe, 2003; Nishi and Bekku, 2015) and jays (Garrulus glandarius brandtii) (Johnson and Adkisson, 1985). Nishi and Bekku (2015) estimated that the spotted nutcracker flew 10.5 km for hoarding seeds, and Kitamura and Nakanishi (2021) revealed seed migration of F. crenata at least 12 km in distance in Hokkaido. Also, European nutcrackers carry seeds as far as 22 km (Vander Wall and Balda, 1977). Thus, a 16 km-long seed dispersal between Hokkaido and Okushiri Island likely occurs. In addition, the cpDNA haplotype mixture such as populations 18, 20, and 31, might have resulted from recent seed dispersal on the island, given that the haplotype A lineage increased dominance more recently than the haplotype B lineage.
Demographic history and northern persistence of F. crenata predated the LGM
The result of the demographic inference suggested that the scenario of the isolation with admixture and secondary admixture was the most likely to explain the genetic structure of F. crenata populations examined in this study (Figure 5). The STRUCTURE analysis by nuclear SSR genotype revealed different coancestries between Hokkaido (including Okushiri Island) and the Tohoku region when K = 2 (Figure 5A), and a distinct difference between Hokkaido and Okushiri Island when K = 3 (Figure 5B). Moreover, PCoA separated the Tohoku region at the first coordinate and secondary separation between Hokkaido and Okushiri Island (Figure 6). Although a clear admixture pattern was not detected in Okushiri Island populations at K = 2, this might be due to recent secondary gene flow from both Hokkaido and Tohoku populations.
The highest probability was obtained from Scenario 5 (Table 2, Figure 2, Supplementary Figure 2). Although time scales for the demographic events were inferred in this study, transformation of generation time to year is still challenging in tree species as they are long-lived with overlapping of generations (Tsuda et al., 2015). Therefore, generation time assumption causes uncertainty in the inference. However, ecological information is still useful to consider generation time, and a feasible substitution for the generation time can be the earliest reproductive age of the species in the northern range. There are several observations for the earliest reproductive age of F. crenata of the northern distribution. In the Kuromatsunai Depression, naturally dispersed F. crenata in the open environment have flowered at 20 and 21 years of age (K. Kitamura and H. Saito, personal observations 2013). An ornamental F. crenata tree at Kuromatsunai Depression flowered at 18 years after the juvenile plantation (H. Saito, personal observations, 2013). Provenance tests of F. crenata at two nurseries of Tokyo University (Chichibu, Saitama, Japan and Furano, Hokkaido, Japan) revealed that the northern originated F. crenata flowered 15 and 17 years after juvenile plantation at the earliest (S. Goto and M. Takahashi, personal observations 2013). On the other hand, the age of reproductive trees at the northernmost population was estimated to be 80 years (Kitamura and Nakanishi, 2021). From these facts, we assumed two different generation time; 20 years at the earliest and 100 years at the fully matured status for F. crenata in this study (Table 3). Thus, the time scales of t1 (admixture from Pop1 to Pop2), t2 (admixture from Pop3 to Pop2), t3 (admixture to generate Pop2) and t4 (splitting time of Pop1 from Pop3) were inferred 35,600 (95% HPD: 1,316 - 150,000), 87,600 (95% HPD: 19,960 - 174,600), 77,400 (95% HPD: 8,620-170,800) and 157,800 (95% HPD: 68,400-198,200) years BP, respectively, under the assumption of generation time of 20 years. When we assumed 100 years for the generation time, they are 178,000 (95% HPD: 6,580 - 750,000), 438,000 (95% HPD: 99,800 - 873,000), 387,000 (95% HPD: 43,100-854,000) and 789,000 (95% HPD: 342,000-991,000) years BP, respectively. The results suggested the divergence between Pop1 (Hokkaido) and Pop3 (Tohoku) predated the LGM even under the assumption of a short generation time of 20 years and even when we considered the 95% HPDs. The divergence time was much earlier when we assumed 100 years for the generation time. Moreover, the demographic inference suggested that the Okushiri Island population was generated by admixture between Hokkaido and Tohoku populations 77,400 (95% HPD: 8,620-170,800) and 387,000 (95% HPD: 43,100-854,000) years BP, under the assumption of the generation time of 20 and 100 years, respectively. The inferred rc value (0.505) suggested that both parental populations equally contributed to the admixture to Okushiri Island population. Notably, asymmetrical secondary admixtures to Okushiri Island population from both parental populations was inferred in this study. However, as the median value of t2 (4,380), time of the secondary admixture from Tohoku to Okushiri Island populations, was slightly larger than t3 (3,870) and their 95% HPDs were overlapped, it is difficult to discuss the time scale of this secondary admixture in detail. On the other hand, the secondary asymmetrical admixture from Hokkaido to Okushiri Island populations was inferred to be 35,600 (95% HPD: 1,316 - 150,000) years BP, under the assumption of the generation time of 20 years, suggesting recent gene flow to Okushiri Island population from Hokkaido, probably in relation to the demographic episode during the last ice age or the LGM. Although it is difficult to estimate locations where ancestors of Hokkaido and Tohoku populations distributed in the past and the secondary admixture could not be a single event nor at a specific time but rather by continuous gene flow over time, these finding shed a new light on past population demography of the northern populations of F. crenata, and these results were different from the previous studies which discussed post-LGM northward expansion to Hokkaido from Tohoku region (Tsukada, 1982; Tomaru et al., 1997; Fujii et al., 2002).
Table 3 Demographic parameters of Scenario 5 obtained by DIYABC.
Recent studies of distribution shift of species based on palaeoecological and/or genetic data revealed the existence of cryptic refugia and northern persistence of forest tree species in more northerly latitudes than initially expectations which assumed southern refugia in lower altitude, not only in Europe and North America but also Japan (Petit et al., 2002; McLachlan and Clark, 2004; Tsuda and Ide, 2005; Petit et al., 2008; Provan and Bennett, 2008; Liepelt et al., 2009; Parducci et al., 2012; Tsuda et al., 2015; Tsuda et al., 2016; Tsuda et al., 2017). Bradshaw et al. (2010) reviewed long-term distribution dynamics of 3 Fagus species; namely, F. sylvatica, F. grandifolia, and F. crenata, and indicated the recent northward spread from several outlying founder populations as well as the location of glacial refugia (Kito and Takimoto, 1999; McLachlan et al., 2005; Magri et al., 2006; Magri, 2008). This hypothesis of northern persistence of cold-tolerant tree species was well supported by species distribution model (SDM) (Svenning et al., 2008) and a combined analysis of the SDM and genetic based-demographic inference (Tsuda et al., 2015). Indeed, a number of palynological studies suggest that there must have been cryptic refugia on the Japan Sea side of Hokkaido (Nakamura and Tsukada, 1960; Nakamura, 1968; Uemura and Takeda, 1987; Kito, 2015; Ooi, 2016). Actually, Tsuda et al. (2015) showed a high possibility of persistence of Betula maximowicziana, a birch species that often distributes in cool temperate forests with F. crenata, during the LGM mainly in the western part of Hokkaido as well as the northern part of the Tohoku region. The distribution of F. crenata in the Pleistocene interglacial periods expanded further north on the Japan Sea side and east to the Pacific Ocean side of Hokkaido (Yano, 1972). Palynological studies revealed that F. crenata existed in Hokkaido (Sakaguchi, 1989; Kito and Takimoto, 1999; Kito and Ohkuro, 2012) and Okushiri Island (Yamanoi, 1992; Hoshino, 2001) before LGM. Therefore, F. crenata might have been distributed further north and east than the present geographic margin in Hokkaido before the straits divided the Japanese Archipelago (Ito, 1987; Ooi, 2016). Also, phytogeographical evidence of the disjunctive distribution of plants on Okushiri Island and the eastern part of Hokkaido supported the notion that temperate species were distributed further north and east than the present geographic margin (Tatewaki, 1954). Moreover, the Okushiri Strait separated Okushiri Island from Hokkaido prior to the Middle Pleistocene, and then the Tsugaru Strait divided Hokkaido from Honshu in the Riss-Würm interglacial period (Ohshima, 1980). These straits could be indeed physical barriers to prevent long-distance gene flow between regions. Thus, it appears that Okushiri Island, besides the Japan Sea coast of Hokkaido, was a possible persistence of F. crenata in the northern range. The similar amount of genetic diversity in the northern populations from Tohoku to Hokkaido regions might be mirrored by long-term persistence of populations as well as no experience of severe bottlenecks by the post-LGM colonization. Indeed, we did not detect any recent bottlenecks in the examined populations in the pilot run of the ABC inference in this study (data not shown), while newly founded northern marginal populations showed the effects of genetic drift (Kobayashi et al., 2013; Kitamura et al., 2015; Kitamura and Nakanishi, 2021). Moreover, F. crenata on Okushiri Island showed a comparative number of private alleles of SSR to Hokkaido, which indicated the existence of relict on this island.
Finally, the results of this study coincided with palaeoecological and vegetation studies as well as geology of the Okushiri Island, and supported the existence of possible persistence in the northern distribution of F. crenata with secondary admixture from the both parental populations in Hokkaido and Tohoku region. This means Okushiri Island included cryptic refugia of F. crenata which persisted for a longer period than expected from the post-LGM expansion.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
Conceptualization, KK and YT. Methodology, KK and YT. Formal analysis, KK and YT. Investigation KK, KN, TM, and MK. Writing – original draft preparation, KK. Writing- review and editing, KK, KN, YT, TM, and MK. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by JSPS KAKENHI (JP17K07852 and JP20K06152) and Core-to-Core Program (Asia-Africa Science Platforms: JPJSCCB20220007) from the Japan Society for the Promotion of Science and the 27th Pro Natura Fund Grant Program from the Pro Natura Foundation Japan.
Acknowledgments
We thank M. Ooue, M. Takahashi, H. Saitou, T. Nagamitsu, A. Nakanishi, and A. Takazawa for their assistance in the field and laboratory work.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.990927/full#supplementary-material
References
Alsos, I. G., Ehrich, D., Eidesen, P. B., Solstad, H., Westergaard, K. B., Schönswetter, P., et al. (2015). Long-distance plant dispersal to north Atlantic islands: Colonization routes and founder effect. AoB Plants 7 (1), 1–19. doi: 10.1093/aobpla/plv036
Asuka, Y., Tani, N., Tsumura, Y., Tomaru, N. (2004). Development and characterization of microsatellite markers for Fagus crenata blume. Mol. Ecol. Notes 4 (1), 101–103. doi: 10.1046/j.1471-8286.2003.00583.x
Austerlitz, F., Dick, C. W., Dutech, C., Klein, E. K., Oddou-Muratorio, S., Smouse, P. E., et al. (2004). Using genetic markers to estimate the pollen dispersal curve. Mol. Ecol. 13 (4), 937–954. doi: 10.1111/j.1365-294X.2004.02100.x
Barrett, S. C. H. (1996). The reproductive biology and genetics of island plants. Philos. Trans. R. Soc. London B 351, 725–733. doi: 10.1098/rstb.1996.0067
Bradshaw, R. H. W., Kito, N., Giesecke, T. (2010). Factors influencing the Holocene history of Fagus. For. Ecol. Manage. 259 (11), 2204–2212. doi: 10.1016/j.foreco.2009.11.035
Chiang, Y. C., Hung, K. H., Schaal, B. A., Ge, X. J., Hsu, T. W., Chiang, T. Y. (2006). Contrasting phylogeographical patterns between mainland and island taxa of the Pinus luchuensis complex. Mol. Ecol. 15 (3), 765–779. doi: 10.1111/j.1365-294X.2005.02833.x
Cornuet, J. M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., et al. (2014). DIYABC v2.0: A software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30 (8), 1187–1189. doi: 10.1093/bioinformatics/btt763
Cornuet, J. M., Santos, F., Beaumont, M. A., Robert, C. P., Marin, J. M., Balding, D. J., et al. (2008). Inferring population history with DIY ABC: A user-friendly approach to approximate Bayesian computation. Bioinformatics 24 (23), 2713–2719. doi: 10.1093/bioinformatics/btn514
DeJoode, D. R., Wendel, J. F. (1992). Genetic diversity and origin of the Hawaiian islands cotton, Gossypium tomentosum. Am. J. Bot. 79 (11), 1311–1311. doi: 10.2307/2445059
Désamoré, A., Laenen, B., González-Mancebo, J. M., Jaén Molina, R., Bystriakova, N., Martinez-Klimova, E., et al. (2012). Inverted patterns of genetic diversity in continental and island populations of the heather. Erica scoparia s.l. J. Biogeography 39 (3), 574–584. doi: 10.1111/j.1365-2699.2011.02622.x
Elleouet, J. S., Aitken, S. N. (2019). Long-distance pollen dispersal during recent colonization favors a rapid but partial recovery of genetic diversity in Picea sitchensis. N. Phytol. 222 (2), 1088–1100. doi: 10.1111/nph.15615
El Mousadik, A., Petit, R. J. (1996). High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) skeels] endemic to Morocco. Theor. Appl. Genet. 92 (7), 832–839. doi: 10.1007/BF00221895
Estoup, A., Jarne, P., Cornuet, J.-M. (2002). Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. Mol. Ecol. 11, 1591–1604. doi: 10.1046/j.1365-294X.2002.01576.x
Evanno, G., Regnaut, S., Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14 (8), 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x
Falush, D., Stephens, M., Pritchard, J. K. (2003). Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics 164, 1567–1587. doi: 10.1093/genetics/164.4.1567
Frankham, R. (1996). Relationship of genetic variation to population size in wildlife. Conserv. Biol. 10 (6), 1500–1508. doi: 10.1046/j.1523-1739.1996.10061500.x
Frankham, R. (1997). Do island populations have less genetic variation than mainland populations? Heredity 78 (3), 311–327. doi: 10.1038/hdy.1997.46
Fujii, N., Tomaru, N., Okuyama, K., Koike, T., Mikami, T., Ueda, K. (2002). Chloroplast DNA phylogeography of Fagus crenata (Fagaceae) in Japan. Plant Systematics Evol. 232 (1-2), 21–33. doi: 10.1007/s006060200024
García-Verdugo, C., Calleja, J. A., Vargas, P., Silva, L., Moreira, O., Pulido, F. (2013). Polyploidy and microsatellite variation in the relict tree Prunus lusitanica l.: How effective are refugia in preserving genotypic diversity of clonal taxa? Mol. Ecol. 22 (6), 1546–1557. doi: 10.1111/mec.12194
García-Verdugo, C., Sajeva, M., La Mantia, T., Harrouni, C., Msanda, F., Caujapé-Castells, J. (2015). Do island plant populations really have lower genetic variation than mainland populations? Effects of selection and distribution range on genetic diversity estimates. Mol. Ecol. 24 (4), 726–741. doi: 10.1111/mec.13060
Goudet, J. (2001). FSTAT, a program to estimate and test gene diversities and fixation indices (Lausanne, Switzerland: University of Lausanne).
Hagiwara, S. (1977). Buna ni mirareru youmenseki no kurain ni tsuite [Clines of leaf area of beech]. Soc. Study Species Biol. / Shuseibutsugaku Kenkyu 1, 39–51.
Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. Biol. J. Linn. Soc. 68 (1-2), 87–112. doi: 10.1006/bijl.1999.0332
Hiraoka, K., Tomaru, N. (2009). Genetic divergence in nuclear genomes between populations of Fagus crenata along the Japan Sea and Pacific sides of Japan. J. Plant Res. 122 (3), 269–282. doi: 10.1007/s10265-009-0217-9
Hoshino, F. (2001). Okushirito chuobu no kahun bunseki [Fossil pollen analysis in central Okushiri Island]. Hoppo Sanso 18, 53–58.
Hubisz, M. J., Falush, D., Stephens, M., Pritchard, J. K. (2009). Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 9 (5), 1322–1332. doi: 10.1111/j.1755-0998.2009.02591.x
Hukusima, T., Nashimoto, M., Watanabe, I. (1984). Phytosociological studies on the beech forest in Hokkaido, Japan. Tech. Bull. Faculty Horticulture Chiba Univ. 33, 117–131.
Huson, D. H., Bryant, D. (2006). Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23 (2), 254–267. doi: 10.1093/molbev/msj030
Hu, L.-J., Uchiyama, K., Shen, H.-L., Saito, Y., Tsuda, Y., Ide, Y. (2008). Nuclear DNA microsatellites reveal genetic variation but a lack of phylogeographical structure in an endangered species, Fraxinus mandshurica, across north-east China. Ann. Bot. 102 (2), 195–205. doi: 10.1093/aob/mcn074
Ishii, H. R., Horikawa, S.i., Noguchi, Y., Azuma, W. (2018). Variation of intra-crown leaf plasticity of Fagus crenata across its geographical range in Japan. For. Ecol. Manage. 429, 437–448. doi: 10.1016/j.foreco.2018.07.016
Johnson, W. C., Adkisson, C. S. (1985). Dispersal of beech nuts by blue jays in fragmented landscape. Am. Midland Nat. 113 (2), 319–324. doi: 10.2307/2425577
Katai, H., Takahashi, M., Hiraoka, K., Yamada, S., Tomaru, N. (2017). Indigenous genetic lineages of Fagus crenata found in the Izu Peninsula suggest that there was one of refugia for the species during the last glacial maximum. J. For. Res. 18 (5), 418–429. doi: 10.1007/s10310-012-0368-8
Katai, H., Takahashi, M., Hiraoka, K., Yamada, S., Yamamoto, S., Kato, K., et al. (2011). Inference of genetic lineages of Fagus crenata populations in Shizuoka Prefecture based on chloroplast DNA and nuclear microsatellite variations. J. Japanese For. Soc. 93 (2), 73–78. doi: 10.4005/jjfs.93.73
Katai, H., Yamada, S., Hiraoka, K., Hoshikawa, T., Tomaru, N., Takahashi, M. (2014). Genetic lineage and diversity of Fagus crenata trees planted in Shizuoka Prefecture. For. Genet. Tree Breed. 3, 101–110.
Kitamura, K., Kobayashi, M., Kawahara, T. (2007). Age structure of wind-felled canopy trees for Siebold's beech (Fagus crenata) in the northernmost population in Karibayama, Hokkaido. J. For. Res. 12 (6), 467–472. doi: 10.1007/s10310-007-0026-8
Kitamura, K., Matsui, T., Kobayashi, M., Saitou, H., Namikawa, K., Tsuda, Y. (2015). Decline in gene diversity and strong genetic drift in the northward-expanding marginal populations of Fagus crenata. Tree Genet. Genomes 11. doi: 10.1007/s11295-015-0857-y. Article Number 36.
Kitamura, K., Nakanishi, A. (2021). Recovery process of genetic diversity through seed and pollen immigration at the northernmost leading-edge population of Fagus crenata. Plant Species Biol. 36 (3), 489–502. doi: 10.1111/1442-1984.12332
Kito, N. (2015). Expansion of beech forest in Tohoku and Hokkaido since Last Glacial. Japanese J. For. Environ. 57 (2), 69–74. doi: 10.18922/jjfe.57.2_69
Kito, N., Ohkuro, Y. (2012). Vegetation response to climatic oscillations during the last glacial-interglacial transition in northern Japan. Quaternary Int. 254, 118–128. doi: 10.1016/j.quaint.2011.06.050
Kito, N., Takimoto, F. (1999). Population growth and migration rate of Fagus crenata during the Holocene in southwestern Hokkaido, Japan. Quaternary Res. 38 (4), 297–311.
Kobayashi, M., Kitamura, K., Matsui, T., Kawano, S. (2013). Genetic characteristics reflecting the population size and disturbance regime of Siebold's beech (Fagus crenata blume) populations at the northernmost distribution. Silvae Genetica 62 (1-2), 1–7. doi: 10.1515/sg-2013-0001
Kobayashi, M., Watanabe, S. (2003). Stand structure of the northern bound population of Fagus crenata, located at Tsubamenosawa, Hokkaido, Japan. Bull. Geo-environmental Sci. 5, 1–23.
Koike, T., Kato, S., Shimamoto, Y., Kitamura, K., Kawano, S., Ueda, K., et al. (1998). Mitochondrial DNA variation follows a geographic pattern in Japanese beech species. Botanica Acta 111 (1), 87–91. doi: 10.1111/j.1438-8677.1998.tb00682.x
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., Mayrose, I. (2015). CLUMPAK: A program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15 (5), 1179–1191. doi: 10.1111/1755-0998.12387
Kuparinen, A., Savolainen, O., Schurr, F. M. (2010). Increased mortality can promote evolutionary adaptation of forest trees to climate change. For. Ecol. Manage. 259 (5), 1003–1008. doi: 10.1016/j.foreco.2009.12.006
Liepelt, S., Cheddadi, R., de Beaulieu, J. L., Fady, B., Gömöry, D., Hussendörfer, E., et al. (2009). Postglacial range expansion and its genetic imprints in Abies alba (Mill.) - a synthesis from palaeobotanic and genetic data. Rev. Palaeobotany Palynology 153 (1-2), 139–149. doi: 10.1016/j.revpalbo.2008.07.007
Magri, D. (2008). Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica). J. Biogeography 35 (3), 450–463. doi: 10.1111/j.1365-2699.2007.01803.x
Magri, D., Vendramin, G. G., Comps, B., Dupanloup, I., Geburek, T., Gömöry, D., et al. (2006). A new scenario for the Quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences. N Phytol. 171 (1), 199–221. doi: 10.1111/j.1469-8137.2006.01740.x
Marshall, T. C., Slate, J., Kruuk, L. E. B., Pemberton, J. M. (1998). Statistical confidence for likelihood-based paternity inference in natural populations. Mol. Ecol. 7 (5), 639–655. doi: 10.1046/j.1365-294x.1998.00374.x
McLachlan, J. S., Clark, J. S. (2004). Reconstructing historical ranges with fossil data at continental scales. For. Ecol. Manage. 197 (1-3), 139–147. doi: 10.1016/j.foreco.2004.05.026
McLachlan, J. S., Clark, J. S., Manos, P. S. (2005). Molecular indicators of tree migration capacity under rapid climate change. Ecology 86 (8), 2088–2098. doi: 10.1890/04-1036
Meirmans, P. G., Van Tienderen, P. H. (2004). GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Mol. Ecol. Notes 4 (4), 792–794. doi: 10.1111/j.1471-8286.2004.00770.x
Miguchi, H. (1996). Dynamics of beech forest from the view point of rodents ecology - ecological interactions of the regeneration characteristics of Fagus crenata and rodents. Japanese J. Ecol. 46, 185–189.
Mogensen, H. L. (1996). The hows and whys of cytoplasmic inheritance in seed plants. Am. J. Bot. 83 (3), 383–404. doi: 10.1002/j.1537-2197.1996.tb12718.x
Nakamura, J. (1968). Palynological aspects of the Quaternary in Hokkaido. V. pollen succession and climatic change since the upper Pleistocene. Research Reports of Kochi University. Natural Sci. 17 (3), 39–51.
Nakamura, J., Tsukada, M. (1960). Palynological aspects of the Quaternary in Hokkaido. I. the Oshima Peninsula. Research Reports of Kochi University. Natural Sci. 9, 117–138.
Namikawa, K., Kitamura, K., Matsui, T., Ishikawa, Y. (2021). The comparative study on the species composition of Fagus crenata Blume forests in Okushiri Island of southwestern Hokkaido and its surrounding area. Vegetation Sci. 38, 175–190.
Nishi, N., Bekku, Y. S. (2015). The spotted nutcracker hoarding seeds of Japanese white pine in Mt. Fuji where Japanese stone pine is not distributed. Strix 31, 113–123.
Ohshima, K. (1980). Recording the Late-Quaternary sea-level change on the topographic feature of the straits of the Japanese islands. Quaternary Res. Tokyo 19 (1), 23–37. doi: 10.4116/jaqua.19.23
Okaura, T., Harada, K. (2002). Phylogeographical structure revealed by chloroplast DNA variation in Japanese beech (Fagus crenata Blume). Heredity 88, 322–329. doi: 10.1038/sj.hdy.6800048
Ooi, N. (2016). Vegetation history of Japan since the last glacial based on palynological data. Japanease J. Histrical Bot. 25 (1-2), 1–101.
Parducci, L., Jørgensen, T., Tollefsrud, M. M., Elverland, E., Alm, T., Fontana, S. L., et al. (2012). Glacial survival of boreal trees in northern Scandinavia. Science 335 (6072), 1083–1086. doi: 10.1126/science.1216043
Pastorelli, R., Smulders, M.J.M., Van’T Westende, W.P.C., Vosman, B., Giannini, R., Vettori, C., et al (2003). Characterization of microsatellite markers in Fagus sylvatica L. and Fagus orientalis Lipsky. Mol. Ecol. Notes 3 (1), 76–78. doi: 10.1046/j.1471-8286.2003.00355.x
Petit, R. J., Brewer, S., Bordács, S., Burg, K., Cheddadi, R., Coart, E., et al. (2002). Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. For. Ecol. Manage. 156 (1-3), 49–74. doi: 10.1016/S0378-1127(01)00634-X
Petit, R. J., Excoffier, L. (2009). Gene flow and species delimitation. Trends Ecol. Evol. 24 (7), 386–393. doi: 10.1016/j.tree.2009.02.011
Petit, R. J., Hu, F. S., Dick, C. W. (2008). Forests of the past: A window to future changes. Science 320 (5882), 1450–1452. doi: 10.1126/science.1155457
Pritchard, J. K., Stephens, M., Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi: 10.1093/genetics/155.2.945
Provan, J., Bennett, K. D. (2008). Phylogeographic insights into cryptic glacial refugia. Trends Ecol. Evol. 23 (10), 564–571. doi: 10.1016/j.tree.2008.06.010
Rosas Escobar, P., Gernandt, D. S., Piñero, D., Garcillán, P. P. (2011). Plastid DNA diversity is higher in the island endemic guadalupe cypress than in the continental tecate cypress. PloS One 6 (1), e16133–e16133. doi: 10.1371/journal.pone.0016133
Sakaguchi, Y. (1989). Some pollen records from Hokkaido and Sakhalin. Bull. Department Geography Univ. Tokyo 21, 1–17.
St Clair, J. B., Howe, G. T. (2007). Genetic maladaptation of coastal Douglas-fir seedlings to future climates. Global Change Biol. 13 (7), 1441–1454. doi: 10.1111/j.1365-2486.2007.01385.x
Stuessy, T. F. (2020). The importance of historical ecology for interpreting evolutionary processes in plants of oceanic islands. J. Systematics Evol. 58 (6), 751–766. doi: 10.1111/jse.12673
Stuessy, T. F., Takayama, K., López-Sepúlveda, P., Crawford, D. J. (2014). Interpretation of patterns of genetic variation in endemic plant species of oceanic islands. Botanical J. Linn. Soc. 174 (3), 276–288. doi: 10.1111/boj.12088
Su, J., Yan, Y., Song, J., Li, J., Mao, J., Wang, N., et al. (2018). Recent fragmentation may not alter genetic patterns in endangered long-lived species: Evidence from Taxus cuspidata. Front. Plant Sci. 871. doi: 10.3389/fpls.2018.01571
Svenning, J. C., Normand, S., Kageyama, M. (2008). Glacial refugia of temperate trees in Europe: Insights from species distribution modelling. J. Ecol. 96 (6), 1117–1127. doi: 10.1111/j.1365-2745.2008.01422.x
Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., Cosson, J. F. (1998). Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7 (4), 453–464. doi: 10.1046/j.1365-294x.1998.00289.x
Takahashi, M., Goto, S., Fukuda, Y., Watanabe, A. (2022). Utility of chloroplast DNA haplotype data for ecological restoration using Fagus crenata seedlings in case of incomplete seed source information availability. Ecol. Res. doi: 10.1111/1440-1703.12351
Takahashi, M., Tsumura, Y., Nakamura, T., Uchida, K., Ohba, K. (1994). Allozyme variation of Fagus crenata in northeastern Japan. Can. J. For. Res. 24 (5), 1071–1074. doi: 10.1139/x94-141
Takayama, K., Lopez-Sepulveda, P., Greimler, J., Crawford, D. J., Penailillo, P., Baeza, M., et al. (2015). Relationships and genetic consequences of contrasting modes of speciation among endemic species of Robinsonia (Asteraceae, senecioneae) of the Juan Fernandez Archipelago, Chile, based on AFLPs and SSRs. N. Phytol. 205 (1), 415–428. doi: 10.1111/nph.13000
Takayama, K., Sun, B. Y., Stuessy, T. F. (2013). Anagenetic speciation in Ullung Island, Korea: Genetic diversity and structure in the island endemic species, Acer takesimense (Sapindaceae). J. Plant Res. 126 (3), 323–333. doi: 10.1007/s10265-012-0529-z
Tanaka, K., Nakamura, T., Tsumura, Y. (1999). Development and polymorphism of microsatellite markers for Fagus crenata and the closely related species, F. Japonica. Theor. Appl. Genet. 99 (1-2), 11–15. doi: 10.1007/s001220051203
Tatewaki, M. (1948). Buna no hokugenkai [Northern limit of Fagus crenata]. Ecol. Rev. / Seitaigaku Kenkyu 11, 46–51.
Tatewaki, M. (1954). Disjunctive distribution of the flowering plants in Hokkaido, Japan. Bull. Soc. Plant Ecol. 3 (4), 250–170.
Tomaru, N., Mitsutsuji, T., Takahashi, M., Tsumura, Y., Uchida, K., Ohba, K. (1997). Genetic diversity in Fagus crenata (Japanese beech): Influence of the distributional shift during the late-Quaternary. Heredity 78 (3), 241–251. doi: 10.1038/hdy.1997.38
Tomaru, N., Takahashi, M., Tsumura, Y., Takahashi, M., Ohba, K. (1998). Intraspecific variation and phylogeographic patterns of Fagus crenata (Fagaceae ) mitochondrial DNA. Am. J. Bot. 85 (5), 629–636. doi: 10.2307/2446531
Tsuda, Y., Chen, J., Stocks, M., Källman, T., Sønstebø, J. H., Parducci, L., et al. (2016). The extent and meaning of hybridization and introgression between Siberian spruce (Picea obovata) and Norway spruce (Picea abies): Cryptic refugia as stepping stones to the west? Mol. Ecol. 25 (12), 2773–2789. doi: 10.1111/mec.13654
Tsuda, Y., Ide, Y. (2005). Wide-range analysis of genetic structure of Betula maximowicziana, a long-lived pioneer tree species and noble hardwood in the cool temperate zone of Japan. Mol. Ecol. 14 (13), 3929–3941. doi: 10.1111/j.1365-294X.2005.02715.x
Tsuda, Y., Nakao, K., Ide, Y., Tsumura, Y. (2015). The population demography of Betula maximowicziana, a cool-temperate tree species in Japan, in relation to the last glacial period: Its admixture-like genetic structure is the result of simple population splitting not admixing. Mol. Ecol. 24 (7), 1403–1418. doi: 10.1111/mec.13123
Tsuda, Y., Semerikov, V. L., Sebastiani, F., Vendramin, G. G., Lascoux, M. (2017). Multispecies genetic structure and hybridization in the Betula genus across Eurasia. Mol. Ecol. 26 (2), 589–605. doi: 10.1111/mec.13885
Tsukada, M. (1982). Late-quaternary shift of Fagus distribution. Botanical Magazine Tokyo 95 (2), 203–217. doi: 10.1007/BF02488586
Uemura, S., Takeda, Y. (1987). Phytogeographical study on the distribution of the main temperate plants composing the natural forests in Hokkaido, Japan. Papers Plant Ecol. Taxonomy to Memory Dr. Satoshi Nakanishi, 259–269.
Vander Wall, S. B., Balda, R. P. (1977). Coadaptations of the Clark's Nutcracker and the piñon pine for efficient seed harvest and dispersal. Ecol. Monogr. 47 (1), 89–111. doi: 10.2307/1942225
Yamanoi, T. (1992). Palyno-flora of middle Miocene sediments of Okushiri Island, southwest Hokkaido. Japanese J. Palynology 38 (2), 106–115.
Yano, M. (1972). On the plant remains from the fossil elephant bearing bed in Tokachi Plain, Hokkaido. Assoc. Geological Collaboration Japan /Chikyu Kagaku 26 (1), 12–19.
Keywords: the northernmost geographic range, chloroplast DNA, nuclear SSR, refugia, structure
Citation: Kitamura K, Namikawa K, Tsuda Y, Kobayashi M and Matsui T (2022) Possible northern persistence of Siebold’s beech, Fagus crenata, at its northernmost distribution limit on an island in Japan Sea: Okushiri Island, Hokkaido. Front. Plant Sci. 13:990927. doi: 10.3389/fpls.2022.990927
Received: 11 July 2022; Accepted: 16 November 2022;
Published: 15 December 2022.
Edited by:
Fang Du, Beijing Forestry University, ChinaReviewed by:
Filipa Monteiro, University of Lisbon, PortugalWataru Ishizuka, Hokkaido Research Organization, Japan
Copyright © 2022 Kitamura, Namikawa, Tsuda, Kobayashi and Matsui. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Keiko Kitamura, a2l0YW1xQGZmcHJpLmFmZnJjLmdvLmpw
†These authors have contributed equally to this work