Microsatellite analysis

The DNA of organelles with maternal inheritance has been used for phylogeography studies in many plant and animal species because of the clear genetic divergence between populations (McCauley, 1995). However, in this study, we used nuclear DNA markers to examine the historical genetic structure because organelle DNAs such as plastid and mitochondrial DNA show paternal inheritance in this species (Ohba et al., 1971; Kondo et al., 1998). We determined the genotypes of all samples using 20 microsatellite markers, including eight genomic (Moriguchi et al., 2003; Tani et al., 2004) and 12 EST (expressed sequence tag) microsatellite loci (Ueno et al., 2012). The genomic microsatellite loci were also used in previous studies (Takahashi et al., 2005; Kimura et al., 2013) and the EST microsatellite loci have been newly included in this analysis. Polymerase chain reaction (PCR) amplification was performed in 6 μL reaction volumes containing about 5 ng of genomic DNA solution, 3 μL of 2× Qiagen Multiplex Master Mix (Qiagen, Hilden, Germany) and each primer at a concentration of 2 μm (the forward primer was labelled with dye in each primer pair). Amplification was performed using a GeneAmp PCR System Model 9700 thermal cycler (PE Applied Biosystems, Warrington, UK). PCR conditions were as follows: denaturation at 95 °C for 15 min followed by 35 cycles consisting of 94 °C for 30 s, annealing at 57 °C for 90 s and 72 °C for 60 s, with a final extension step at 60 °C for 30 min. The PCR products and DNA size marker (LIZ600; Life Technologies, Foster City, CA, USA) were separated by electrophoresis on an ABI 3100 Avant Genetic Analyser (Applied Biosystems), and their fragment sizes were determined using GENEMAPPER version 3.7 analysis software (Applied Biosystems).

Genetic diversity and characteristics

To assess genetic diversity across all populations, the total number of alleles detected (TA), gene diversity in the total population (H_T), average gene diversity within populations (H_S; Nei, 1987) and observed heterozygosity (H_O) were calculated across all populations at each locus and over all loci using FSTAT v. 2.9.3.2 (Goudet, 2002). The fixation index (F_IS; Weir and Cockerham, 1984) was also calculated using FSTAT v. 2.9.3.2 (Goudet, 2002) across all populations at each locus and across all loci to measure departure from Hardy–Weinberg equilibrium (HWE).

Similarly, to assess genetic diversity within each population, the number of alleles (A), allelic richness (Ar; El Mousadik and Petit, 1996), H_O and unbiased expected heterozygosity (Nei, 1987) were calculated at each locus and across all loci in each population using FSTAT v. 2.9.3.2 (Goudet, 2002). Ar was calculated for a standard sample size of nine (18 gene copies), which was the smallest sample size, with complete genotypes at all 20 loci, among populations. We also calculated the number of rare alleles (RA; defined as alleles with a frequency <1 % in the total population) and private alleles (PA; unique to one population) per individual in each population. The fixation index (F_IS; 1 – H_O/H_E) was calculated at each locus and over all loci in each population to measure departure from HWE. The deviations of F_IS from zero were tested in each population by false discovery rate test using the ‘stats’ package in R 2.15.0 (R Development Core Team, 2012).

We also evaluated whether the populations studied had experienced a reduction in effective population size through founding events or population bottlenecks using Wilcoxon's signed-rank test (one-tailed) with the null hypothesis, H_E < H_EQ, and the alternative hypothesis, H_E > H_EQ, under the assumption of mutation–drift equilibrium in the infinite allele model (IAM) and the stepwise mutation model (SMM) with 1000 simulation iterations using BOTTLENECK v. 1.2.02 (Piry et al., 1999).

Bayesian genetic clustering

We performed Bayesian cluster analysis using STRUCTURE 2.3.1 (Pritchard et al., 2000; Falush et al., 2003a) to assign individuals to clusters based on their multilocus genotypes. We first estimated the optimal number of genetic clusters (K) using 761 individuals from 37 natural populations of C. japonica. Using an admixture model with correlated allele frequencies and location information (LOCPRIOR model), 20 independent runs were conducted for each K (putative cluster number, from 1 to 7) with 100 000 iterations after a 10 000-step burn-in. All other parameters were set at their default values. The optimal number of clusters, K, was identified using lnPr (X|K) (Pritchard et al., 2000) or delta K; the rate of change of lnPr (X|K) between successive K values (Evanno et al., 2005) was maximal. The genetic relationships among the genetic clusters were evaluated based on net nucleotide distance calculated in STRUCTURE, and a Neighbor–Joining (NJ) tree based on the net nucleotide distances was generated using POPULATIONS v. 1.2.30 (Langella, 1999). The net nucleotide distance is the average probability that the two alleles in a pair from two different genetic clusters are different (Falush et al., 2003b).

CLUMPP1.1.1 (Jakobsson and Rosenberg, 2007) was used to calculate the average membership coefficient (Q value) for each individual by aligning and converging the results of the 20 runs described above. We also calculated the average membership coefficient (Q value) for each individual excluding the Yakushima Island population (K = 3) for ABC analyses (see below) because the genetic cluster of Yakushima Island was particularly different from the other gene pools. In addition, principal co-ordinate analysis (PCoA) we conducted based on genetic distance in GeneAlEx 6·41 (Peakall and Smouse, 2006) to visualize the relationships among populations.

Approximate Bayesian computation (DIY ABC)

In the absence of a fossil record, the time scales which estimated solely on genetic data may not be as precise as geological time scales, but are good enough to make it possible to generate the temporal hierarchies required for understanding interspecific gene flow and demographic dynamics cause by distributional changes (Wakeley, 2008). For probabilistic analysis among alternative hypotheses for the history of gene pool divergence, we used the approximate Bayesian computation procedure (Beaumont et al., 2002) in DIY ABC v. 1.0.4.39 (Cornuet et al., 2008, 2010). We tested the following six competing scenarios (for details, see Fig. 2 and Supplementary Data Fig. S1).

Open in a separate window

Fig. 2.

Distribution of cluster memberships at individual and population levels in the 37 Cryptomeria japonica populations obtained by STRUCTURE analysis (Pritchard et al., 2000). Number of cluster (K) is 4. Each horizontal bar in the histograms represents the proportion of cluster memberships for each individual. The pie charts superimposed on the maps show the average proportions of cluster memberships among the individuals sampled in each population.

Genetic structure

The population genetic structure was explored by the STRUCTURE analysis. The values of log-likelihood of the data, ln P(X|K), reached a plateau when the numbers of clusters were K > 4 (Supplementary Data Fig. S4). Furthermore, the values of delta K, based on the rate of change of ln P(X|K) between successive K values, were highest when K was 2. Therefore, we examined the proportional membership of each cluster of each individual and population at K = 4 and K = 2. When we set K to 4, the distribution of cluster membership exhibited not only the pattern of population divergence between the Japan Sea and Pacific Ocean sides, but also genetic sub-structures among populations along each side (Fig. 2). On the Japan Sea side, the membership proportion of gene pool B1 was higher in the populations from northern Tohoku district (population nos 1–7) to Sado Island in Niigata Prefecture (nos 11, 12), whereas that of gene pool B2 was higher in the populations (nos 17, 26–33) from the Chubu to Chugoku districts (Fig. 2). Similarly, among the populations along the Pacific Ocean side, the populations (nos 10, 18–25) from Tohoku to Shikoku district had higher proportions of gene pool B3, whereas the Yakushima Island populations (nos 34–37) in Kyushu district had higher proportions of gene pool B4. High admixture of several gene pools was found in Kyoto, Toyama Prefecture and southern Tohoku district (nos 8, 9, 14–16, 22). The F-values of gene pools B1, B2, B3 and B4 were 0·05937, 0·04393, 0·02779 and 0·05026, respectively (Fig. 2). The NJ tree using net nucleotide distance showed that gene pool B4 was markedly divergent from the other three gene pools (Fig. 2). The distribution of membership of the two clusters (A1 and A2) was found to be geographically structured between the populations along the Japan Sea and Pacific Ocean sides (Supplementary Data Fig. S2). The proportion of gene pool A1 (blue) was higher in the populations of northern Tohoku district, whereas that of gene pool A2 (red) was higher in populations of Yakushima Island in the Kyushu district. The F-values of gene pools A1 and A2 were 0·02810 and 0·03779, respectively (Fig. S2).

This result of population genetic structure by STRUCTURE analysis corresponded to that of PCoA analysis (Fig. 3). The first two axes accounted for 65·6 % of the total variation. The first axis reflected Yakushima Island, other Pacific Ocean side populations and Japan Sea side populations. The second axis reflected northern Tohoku populations and other populations on the Japan Sea side (Fig. 3).

Open in a separate window

Fig. 3.

Genetic similarity among 37 populations, as depicted from principal co-ordinate analysis. Pie charts represent the proportion of cluster memberships in each population using STRUCTURE analysis with K = 4 (see Fig. 2). Gene pool B1 = blue, gene pool B2 = sky blue, gene pool B3 = orange, gene pool B4 = red.

ABC analysis

A comparison of posterior probabilities of the six scenarios using local linear regression indicated that Scenario 6 had the highest support, with a probability of 0·811 (Table 2). Scenario 6 had the smallest number of marginal summary statistics (Supplementary Data Table S4). ABC analysis provided good support for Scenario 6, with the same time divergence of four gene pools and no bottleneck (Supplementary Data Fig. S1, the type I error rate = 0·22).

The median value of the divergence time (t2) was estimated as 507 generations (Supplementary Data Figs S1 and S5; 95 % CI = 161–1170). For Scenario 6, the median values of effective population size of N_B1, N_B2, N_B3, N_B4 and N_A before population divergence were 4520, 8040, 8840, 5330 and 2630, respectively (Table 3). Posterior distributions are shown in full in Table 3.

Prediction of divergence time using ABC analysis

The ABC analysis using 14 SSR loci supported Scenario 6, in which four gene pools separated at the same time, rather than gradual divergence scenarios. The time of divergence (t2) was predicted as 507 generations. These results suggested that four ancestral populations were present before approx. 500 generations. This divergence time is dependent on the generation time of this species, but the generation time of woody species can be much longer than their maturation age (Petit and Hampe, 2006). As C. japonica begins to produce cones at 20 years of age in artificial open habitats (e.g. Endo et al., 2012), we therefore assumed a generation time for this species of approx. 150 years under natural conditions. If we assume the generation time to be 150 ± 50 years, the divergence time for the four gene pools could be 76 050 ± 25 305 years bp. Although the estimated divergence time also changes by the mutation rate, the posterior distribution of model parameters under the most likely scenario was used to make inferences about the timing of divergence before the LGM. The four gene pools of C. japonica probably indicated refugia populations at the LGM.

Climatic conditions controlling C. japonica distribution

The distribution of the potential habitats for C. japonica was predicted on the basis of the GAM using the present distribution data of C. japonica and a current climatic database. The prediction accuracy of the potential habitats was evaluated as ‘good’ based on the AUC. Species distribution modelling showed that the ecological niche of C. japonica is constrained primarily by winter precipitation (Bio19, 50·2 %). The distribution probability was low in Kanto and Hokkaido districts (Fig. 4). In these areas, climatic variables also showed low coefficient values (mean temperature of summer and winter precipitation in Kanto, mean daily minimum temperature of winter and summer precipitation in Hokkaido; Supplementary Data Fig. S6). Tsukada (1982) noted that the probable glacial refugia in Cryptomeria can be characterized by climatic variables: annual precipitation of at least 1200 mm, annual temperature >5·0 °C, minimum mean January temperature of –7·0 to –8·2 °C and mean July temperature of 20·0–25·0 °C. This climatic regime is almost equivalent to that at the northern limit of planted Cryptomeria forest in Hokkaido today (Tsukada, 1980). Our results were consistent with the report by Tsukada (1982). In addition, the current distribution of the predicted potential habitats covered the natural distribution map of C. japonica by Hayashi (1960). In this model, past human activity, such as logging of natural populations, was not considered; therefore, the predicted current potential habitat area would be wider, especially in the lowlands, if we took past human activity into consideration, but no contradiction occurred with the current positions of natural populations.

Comparing pollen data with predicted potential habitats for the LGM, the MIROC simulation was suggested to be more realistic than CCSM in the case of Japan (Supplementary Data Fig.S3). Calculating the distances between C. japonica pollen data and the projected potential habitat during the LGM showed that the projected potential habitat on the MIROC is reasonably accurate. Potential glacial refugia of C. japonica were predicted not only in western Japan, e.g. around Oki Island, Kii Peninsula and Yakushima Island, but also in lowlands of Tohoku district, such as Oga Peninsula. These results suggested the presence of putative refugia in the coastal areas of Tohoku district.

Many suitable habitats for this species were predicted in Kyushu district at the LGM; however, these natural populations were thought to have been destroyed by volcanic eruptions. For example, a large-scale eruption of Mt. Aso occurred 90 000 years ago and the Kikai caldera erupted violently in the southern Kyushu district 7300 years ago (Machida and Arai, 1992), causing catastrophic damage in these areas where plant life was thought to have been virtually annihilated. Although the small old forests of C. japonica distributed on Kyushu Island, this population might be established by the immigration from Chugoku populations, such as Azouji and others, after the extinction of natural forest on Kyushu Island (Takahara, 1998). The supporting result was obtained by the genetic analysis using SNP markers (Tsumura et al., 2012).

Evidence of cryptic refugia of C. japonica in northern Japan

Bayesian clustering analysis with the highest log-likelihood values also divided into four gene pools: northern Tohoku district (B1); from Chubu to Chugoku districts (B2); from Tohoku to Shikoku district on the Pacific Ocean side (B3); and the Yakushima Island population (B4). The results of PCoA analysis also supported these gene pools. Axis 1 (37·7 %) represents genetic differences among Yakushima Island populations, the other Pacific Ocean side populations and Japan Sea side populations, and axis 2 (27·9 %) represents genetic difference between the northern Tohoku and the other Japan Sea side populations (Fig. 3).

Relatively high levels of within-population gene diversity were detected in the populations of Tohoku district (H_E: average 0·604, range 0·581–0·642) in comparison with the other districts (H_E: average 0·614, range 0·570–0·643). In addition, no evidence of a recent bottleneck was observed in Tohoku district (except the Oga population). These results suggest the potential of northern cryptic refugia and/or the potential of admixture events from several refugia between populations in Tohoku district. The results of STRUCTURE analysis in K = 4 indicated the relatively purity of gene pool B1 in northern Tohoku populations (Fig. 2). We also found that some populations in northern Tohoku possessed private alleles and/or rare alleles. This genetic structure suggested the presence of cryptic refugia in northern Japan at the LGM. Recent phylogeographic studies on temperate trees and animals in Asia suggested that small populations may have survived in the northernmost part or high elevation mountain area and separated from the southern large glacial refugia (Tsuda and Ide, 2005; Tsumura et al., 2007; Ohnishi et al., 2009; Opgenoorth et al., 2010). Tsukada (1982) predicted probable refugia at Toyama Bay using environmental data, such as annual precipitation and temperature, and fossil pollen records. In this case, the C. japonica forest had to expand at a rate of 130 m per year from Toyama Bay (near population no. 14) to Akita Prefecture (near population no. 6) after the LGM to establish the current distribution (Tsukada 1980). In contrast, a few studies found a low frequency (approx. 10 %) of fossil pollen of C. japonica 12 000 years ago at Oga Peninsula (Yoshida and Takeuti, 2009). Moreover the fossil pollen of C. japonica was found in northern Tohoku district at least 5000–6000 years ago (e.g. Shimokita Peninsula; Kawamura, 1979; Yamanaka et al., 1990). These results correspond to the consideration of Takahara (1998) that scattered and limited populations could have survived in geographically isolated pockets of favourable habitat within a forest mosaic dominated by boreal taxa in northern Japan during the LGM. Our SDM also predicted that this species could have been distributed in northern Tohoku district (the coastal side from population no. 9 and Oga Peninsula) during the LGM, although population sizes in these areas were limited. Tsukada (1982) suggested that small isolated refugial populations could hardly persist in the harsh sub-arctic climate, such as the northern region of Japan, during the full glacial period. However, the effective population size of gene pool B1 was sufficient to maintain populations and genetic diversity (median 4650, 95 % CI 1640–8630), although this number was the smallest compared with other gene pools.

Frequent layering from trunks and/or branches occurs in this species, especially in snowy regions (Kimura et al., 2013). As there are large amounts of snowfall over the mountainous regions on the Japan Sea side, the main regeneration system for C. japonica in these areas is generally accepted to be layering, which occurs when a flexible branch or stem is bent to the ground by snow pressure (Shimizu et al., 2002; Hirayama and Sakimoto, 2008). Several studies using genetic markers have shown a higher frequency of layering in regions with heavy snowfall (Taira et al., 1997; Moriguchi et al., 2001; Hirayama and Sakimoto, 2008; Kimura et al., 2013). The reproductive strategy of C. japonica was found to be strongly affected by genetic factors and snow depth (Kimura et al., 2013). These results suggest that clonal reproduction allows persistence under harsh conditions, such as in the LGM, where establishment of small individuals, e.g. seedlings, is disturbed by heavy snow. Clonal reproduction as a reproductive strategy in C. japonica contributes to maintenance of genetic variation by increasing the life span of rare genotypes with unique alleles and by preserving the few seedlings that are able to establish under conditions of frequent disturbance. Such a reproductive strategy may explain the persistence of the small refugial populations in northern Tohoku district.

Potential of anagenetic divergence on Yakushima Island

The results of STRUCTURE analysis in K = 4 also indicated the specificity of gene pool B4 in Yakushima Island populations (Fig. 2). Abundant Cryptomeria pollen occurrence (>40 %) was recorded during the Pliocene on southern Kyushu Island, although sediments derived from the full- and late-glacial periods have yielded neither Cryptomeria pollen nor its microscopic remains, suggesting that it became extinct in Kyushu Island before the LGM (Tsukada, 1980, 1982a). The main island (Honshu), Shikoku and Kyushu of the Japanese Archipelago were connected and formed one large island during the LGM. However, the Kikai caldera near Yakushima Island erupted violently 7300 years ago (Machida and Arai, 1992), causing catastrophic damage to the forests in Kyushu district and Yakushima Island (Kimura et al., 1996). Fossil pollen data indicated that Yakushima Island was one of the main refugia at the LGM (Tsukada, 1982a). Therefore, Yakushima Island may have been isolated from the other C. japonica populations for a long time.

Anagenesis (also known as phyletic speciation) has attracted attention as a speciation process, especially on oceanic islands, in which the initial founder populations simply diverge over time without further specific differentiation (Stuessy et al., 1990, 2006; Whittaker et al., 2008; Takayama et al., 2012; Zou et al., 2013). The clear genetic distinctions and similar genetic diversity were found between common species distributed on the Asian continent and endemic species distributed on islands as evidence of anagenetic speciation (Takayama et al., 2012; Zou et al., 2013). Bayesian clustering analyses showed a relatively high F-value in gene pool B4, suggesting an episode of genetic drift and/or selection during colonization. No clear evidence of a bottleneck was detected based on allelic frequency distribution and excess of observed heterozygosity; indeed, genetic diversity was high compared with other populations. Both PA and RA were also found in Yakushima Island populations. The old stumps have remained in extremely fresh condition in spite of the high levels of precipitation and humidity on the island due to their high resin contents (Toda and Sato, 1969; Takahashi et al., 2008). This result seems to show the difference in wood traits, although the effect of environmental condition was not avoided. Needle shape is also different compared with those of trees from other regions (M. K. Kimura, unpubl. data). The needles of individuals from Yakushima Island are short, opened and hard regardless of planted environment. These results probably suggest an anagenetic process on Yakushima Island. In fact, many endemic species on Yakushima Island are found in several taxa (Yahara et al., 1987).

Admixture-like structure of C. japonica

High admixture of several gene pools was found in Kyoto, Toyama Prefecture and southern Tohoku district (population nos 8, 9, 14–16, 22). This result suggests that multiple refugia and/or admixture events, with an ambiguous genetic cline along latitudes, would be found in C. japonica. We found potential cryptic refugia in Tohoku district (B1), further north than previous pollen records (Tsukada, 1982a). On the main Island of Japan (Honshu), the other two clear gene pools (B2 and B3) were also found from Chubu to Chugoku district on the Japan Sea side and from Tohoku to Shikoku district on the Pacific Ocean side, reflecting a refugial area (Wakasa Bay, Kii Peninsula and Izu Peninsula) based on fossil pollen analysis (Tsukada, 1982a). Species distribution modelling also predicted that these areas would be suitable with high probability for C. japonica at the LGM. It has numerous geographical barriers (e.g. high mountain ranges, valleys and seaways), which are expected to have influenced the distribution patterns and past migration of tree species in Japan. For example, the high mountain range divides the main island of Honshu into the Pacific Ocean and the Japan Sea sides. However, it is thought that secondary contacts could easily have occurred in some areas of lower mountains, because the surroundings of the locations of populations showing admixture have no geographical barrier such as higher mountains.

Evidence of multiple refugia, such as high genetic diversity in marginal populations or contact populations, has been reported in many species (e.g. Hewitt, 2000; Petit et al., 2003). The existence of multiple refugia may have played important roles in maintaining the genetic diversity in these areas. Admixture of genetically differentiated populations would result in high genetic diversity in the newly established populations comparable with other populations. In admixed zones, elevated genetic diversity has been reported in several tree species, e.g. Fagus sylvatica, Fraxinus excelsior, Pinus resinosa and Kalopanax septemlobus (Comps et al., 2001; Walter and Epperson, 2001; Petit et al., 2003; Heuertz et al., 2004; Sakaguchi et al., 2011).

The authors are grateful to Dr T. Nagamitsu for valuable discussions. We also thank Y. Taguchi, M. Koshiba and Y. Komatsu for their help in DNA analyses. This study was partly supported by a grant from the Forest Agency of Ministry of Agriculture, Forestry and Fishery, and by the Programme for the Promotion of Basic and Applied Research for Innovations in Bio-oriented Industry. M.K.K., K.U., Y.M. and Y.T. conceived and designed the experiments. M.K.K. and L.S.J.M. performed molecular analyses and statistical analyses. K.N. estimated current and past distributions by species distribution modelling. M.K.K. and Y.T. wrote the manuscript; the other authors provided editorial advice.