Genetic variation and structure of Ubame oak, Quercus phillyraeoides, in Japan revealed by chloroplast DNA and nuclear microsatellite markers

ABSTRACT

Genetic variation and structure of Ubame oak (Quercus phillyraeoides A. Gray), a shrub tree on the Pacific coast of Japan, were examined to elucidate historical population dynamics using five chloroplast DNA (cpDNA) marker regions and 11 nuclear microsatellite loci. Three cpDNA haplotypes (A, B and C) were identified in a screen of 41 populations across the entire distribution range in Japan. Haplotype A was the most prevalent and was found in the entire range. Haplotype B was locally restricted to the Kii Peninsula, while haplotype C was restricted to the southwestern part of Japan. These haplotypes corresponded with the nuclear genetic constitution revealed by microsatellite markers. Two genetically differentiated major groups were identified by STRUCTURE analysis applied to 536 individuals from 28 populations, and they mostly corresponded with the two major cpDNA haplotypes, A and C. These populations were further divided into three geographically identified groups: group 1 in the area including Kanto and Tokai regions, the Kii Peninsula and the Muroto-misaki Cape on the Pacific coast; group 2 in the Bungo-suido Channel area; and group 3 in southern Kyushu and Okinawa. Populations in the Seto Inland Sea were divided into two groups: one was included in group 1 and appears to have originated from last glacial maximum (LGM) refugia located in the Kii Peninsula, while the other was included in group 2 and appears to have originated from LGM refugia located in southern Kyushu. These groups can be considered as conservation units for the preservation of unique seashore ecosystems, or as a seed source to foster coastal protection forests and next-generation production forests. Considerable care should be taken to protect isolated populations that may be specialized to unique local environments, such as those on the islands of Koshikijima and Izenajima.

INTRODUCTION

Population genetic analysis using molecular markers and assisted by recent advances in computationally intensive statistical methods has provided us with a deeper understanding of the phylogeographic pattern and demographic history of natural populations for a wide range of taxa. This has given rise to new research fields, such as molecular ecology, and has become the theoretical background of practical fields, such as conservation genetics and landscape genetics (Holderegger et al., 2006). Population genetic techniques have also been applied to forest tree populations in various climate classes in the world. Among them, microsatellite markers are the most commonly used technique because of their codominant nature and the large number of available loci with abundant variation (Balloux and Goudet, 2002). Microsatellites are essentially neutral and the allele frequency changes are based on purely stochastic processes (Li, 1997). Thus, the markers have great potential to detect genetic structures shaped by stochastic demographic processes, such as population expansion, colonization, isolation and reunion, that are associated with historical climate changes. This can often be useful for predicting potential natural vegetation (Tüxen, 1956), as well as setting outlines for seed and seedling transfer zones to define conservation units and management units (Manel et al., 2003).

The Japanese archipelago extends approximately 3,000 km from the north to the south in a narrow arc with a maximum width of 250 km. Present-day vegetation in Japan, like other temperate climate forests worldwide, has been largely influenced by past climate changes, especially the repeated ice ages in the late Quaternary period (the last one million years) (Hewitt, 2000). At the last glacial maximum (LGM) at about 20,000 years ago, annual temperatures decreased 6–8 ℃ and the sea level dropped 80–140 m below the present level in the temperate zones in East Asia (Yonekura et al., 2001). This has had a profound effect on present-day vegetation: vegetation shifted toward the south or lower altitudes and often became isolated into small refugia (Comes and Kadereit, 1998). In the course of Earth’s warming since the LGM, northward and higher altitudinal vegetation range shifts have occurred (Comes and Kadereit, 1998). These processes must have been accompanied by population expansion and founding events during a series of colonizations, consequently altering genetic variation and structure.

Forest in Japan is largely divided into two types, the temperate broad-leaved deciduous forests in the north and the warm-temperate broad-leaved evergreen forests in the south. The transition zone occurs in eastern Honshu, with an isotherm of 85 ℃ on the warmth index (Kira, 1948). The broad-leaved evergreen forests contain many plant elements common to eastern China and Taiwan, which constitute the Sino-Japan floristic region (Kitamura, 1957). These forests are sometimes called lucidophyllous forests – forests with shiny leaves (Kira, 1977). Ubame oak, Quercus phillyraeoides A. Gray, is a short shrubby evergreen oak that grows on the Pacific coasts of Honshu from the Boso Peninsula in the northeast to the south on Izena Island, Okinawa. It is also abundant along the coasts of the Seto Inland Sea in western Japan, which separates Kyushu and Shikoku from Honshu (Fig. 1) (Horikawa, 1972), as well as in southern and eastern China below 34 °N (Xie et al., 2011). It has small hard leaves with slight serration and grows on dry, rocky, windy hills facing oceans. Quercus phillyraeoides often forms pure stands in its most favorable habitats, but commonly grows with other shrubby oceanic tree species, such as Pittosporum tobira, Rhaphiolepis indica and Eurya emarginata, and forms one of the typographic forest communities of lucidophyllous forests called, in socio-botanical terms, “Pittosporeto-Quercetum phillyraeoidis Suz.-Tok. et Hatiya 1951” (Miyawaki et al., 1983). Quercus phillyraeoides is the only species that belongs to the Ilex group (evergreen oaks in the subgenus Quercus) in Japan, whereas in China, at least 21 species, mostly endemic, are known in this group (Yang et al., 2017). The Chinese oaks in Ilex are further divided into three morphology-based sections (i.e., Heterobalanus, Englariana and Echinoilepides) (Pu et al., 2002; Denk and Grimm, 2009). Quercus phillyraeoides belongs to the section Englariana with six other species. A molecular genetic study revealed that Q. phillyraeoides is the sister species of Q. acrodonta in the section Englariana (Yang et al., 2017). These findings strengthen the inference that Q. phillyraeoides originated in mainland China, probably on the Yunnan-Guizhou Plateau, during the Miocene and migrated to Japan afterward (Xie et al., 2011).

Fig. 1.

Collection sites of Quercus phillyraeoides and chloroplast DNA haplotype distribution. The Kii Peninsula area is enlarged in the inset on the lower right. Characteristic locations are shown as: a, the Boso Peninsula; b, the Izu Peninsula; c, the Kii Peninsula; d, the Kii-suido Channel; e, the Bungo-suido Channel; f, the Seto Inland Sea; and g, Bisan-seto.

The wood of Q. phillyraeoides is one of the hardest amongst Japanese trees, and is used for the production of a high-quality charcoal called Bincho-tan. Annual production of Bincho-tan charcoal has been maintained at 3,000 tonnes over the last two decades, mostly shared by Wakayama and Kochi prefectures (Forest and Forestry Statistics Handbook, 2016). Although there is not a large supply of wood for making Bincho-tan charcoal, it is most commonly harvested from natural stands and often depletes available trees. Raising next-generation forests for a sustainable supply has been a serious issue for local people involved in the industry (Consortium of Japan Bincho-tan meeting, 2011). Quercus phillyraeoides is also commonly used for gardens and hedges, and as roadside trees in western Japan because of its luxuriant foliage and vigorous branching, even in poor soils. Large parts of original Japanese coastal areas have been utilized for civilian purposes, such as living, industrial, and agricultural lands that include large areas of artificial forests to protect coastlines, which have generally been made on flat ground including sandy beaches. Many areas of coastal vegetation on rocky cliffs, however, still remain because of limited human access and are expected to have retained their original composition of natural vegetation.

The purpose of our study was to examine the genetic variation and structure of Q. phillyraeoides as one of the key species of unique seaside ecosystems, and to elucidate how the species has established itself along the coastal areas of southwestern Japan. This will enable us to outline genetically differentiated groups across the entire range of habitats and will provide important information for designating seed and seedling transfer zones, which may be the first step in the conservation and management of economically and environmentally important trees such as Q. phillyraeoides. We examined genetic variation using nuclear microsatellite markers and chloroplast DNA (cpDNA) sequences in non-coding and coding regions. Microsatellites are biparentally transmitted and, with a higher mutation rate, tend to show changes in genetic structure that occur in a relatively short time. As the chloroplast is maternally inherited in angiosperm plant species, the genetic structure revealed by cpDNA variation will reflect the expansion process mediated by seed migration. By combining the results obtained from markers with different modes of inheritance, we expect to obtain a comprehensive view of the demographic history of Q. phillyraeoides in Japan.

MATERIALS AND METHODS

Plant materials and DNA extraction

Tree leaves were collected from 743 individuals in 41 populations across the entire distribution range of Q. phillyraeoides in Japan, including 19 populations used in a previous study (Liu et al., 2013). Collection sites are listed in Table 1 and shown geographically in Fig. 1. Collected leaves were kept dry in plastic bags with silica gel, and 50 mg of leaf tissue was crushed into powder using a TissueLyser (Qiagen, Hilden, Germany). DNA was extracted using the DNeasy Plant Mini Kit (Qiagen). The quality and concentration of DNA was confirmed by electrophoresis on a 0.8% agarose mini-gel along with a λ-HindIII digest as a DNA size standard. Extracted DNA was stored at −30 ℃ until used.

Table 1. Sampling sites and cpDNA haplotype frequency in Quercus phillyraeoides

Population			Coordinates		Sample sizea	cpDNA haplotype
Name	Code	Area	Latitude	Longitude	n	A	B	C
Iwai	IW	Kanto	35.1019	139.8319	32 (5)	5	0	0
Taraimisaki	TR	Tokai	34.6196	138.8991	10 (2)	2	0	0
Nakagi	NK	Tokai	34.6197	138.8164	29 (2)	2	0	0
Takatoriyama	TK	Tokai	34.7077	138.7576	4 (2)	2	0	0
Toi	TI	Tokai	34.9220	138.8145	7 (2)	2	0	0
Ogasayama	OG	Tokai	34.7455	137.9711	26	5	0	0
Hazu	HZ	Tokai	34.6975	136.9702	31	5	0	0
Toba	TB	Kii Peninsula	34.4547	136.8711	20 (2)	0	13	0
Kumano	KM	Kii Peninsula	33.8947	136.1273	9	8	0	0
Kiho	KH	Kii Peninsula	33.7690	136.0278	8	7	0	0
Kushimoto	KS	Kii Peninsula	33.4361	135.7622	20 (5)	17	3	0
Susami	SS	Kii Peninsula	33.5205	135.5406	6 (2)	3	3	0
Shirahama	SR	Kii Peninsula	33.6621	135.3427	23 (1)	20	0	0
Inami	IN	Kii Peninsula	33.7813	135.2390	5 (1)	4	0	0
Mihama	MH	Kii Peninsula	33.8867	135.0656	16 (2)	16	0	0
Minabe	MN	Kii Peninsula	33.8327	135.3454	5 (1)	4	0	0
Akizugawa	AK	Kii Peninsula	33.8031	135.4274	19 (1)	19	0	0
Nakahechi	NH	Kii Peninsula	33.7590	135.4983	4	4	0	0
Muroto-misaki	MR	Shikoku Pacific	33.2900	134.1519	19 (5)	19	0	0
Suma	SM	Seto Inland Sea	34.6419	135.0919	24	5	0	0
Korahachiman	KR	Seto Inland Sea	34.7230	134.2386	20	0	0	5
Tsuyama	TY	Seto Inland Sea	34.9436	134.0125	40	5	0	0
Syodosima	SY	Seto Inland Sea	34.5097	134.3522	42	5	0	0
Suga	SG	Seto Inland Sea	34.3455	133.9036	25 (2)	0	0	5
Toomiyama	TM	Seto Inland Sea	34.3883	134.1425	20	5	0	0
Yoshima	YS	Seto Inland Sea	34.3858	133.8194	19	4	0	2
Nagahama	NG	Seto Inland Sea	34.3333	132.9772	20 (1)	5	0	0
Iguchijima	IG	Seto Inland Sea	34.2880	133.0744	7	5	0	0
Oomishima	OM	Seto Inland Sea	34.2241	132.9877	9	6	0	0
Ooshima	OS	Seto Inland Sea	34.1172	133.0207	7	6	0	0
Nogutsunajima	NT	Seto Inland Sea	33.9763	132.6889	7	0	0	6
Sadamisaki	SD	Seto Inland Sea	33.3763	132.1161	20 (2)	0	0	5
Toshima	TS	Seto Inland Sea	33.2063	132.3691	20	0	0	5
Ashizuri-misaki	AS	Shikoku Pacific	32.7816	132.9550	20 (2)	5	0	0
Usuki	US	Kyushu	33.1394	131.8108	20	0	0	5
Satamisaki	ST	Kyushu	31.0461	130.6775	19	0	0	5
Yakushima	YK	Kyushu	30.2779	130.4152	2	1	0	1
Noma	NM	Kyushu	31.3994	130.2277	20	3	0	2
Koshikijima	KK	Kyushu	31.3994	130.2280	45 (8)	8	0	0
Ryugadake	RY	Kyushu	32.4063	130.3808	20 (5)	8	0	0
Izenajima	IZ	Okinawa	26.9063	127.9402	24 (5)	0	0	5
				Sum	743 (58)	215	19	46

^a   Numbers in parentheses indicate the number of samples in which all five cpDNA regions were sequenced.

Chloroplast DNA sequencing

Five cpDNA regions were amplified by PCR using universal primers: trnL-trnF spacer (Taberlet et al., 1991), atpB-rbcL spacer (Terachi, 1993), matK gene (Demesure et al., 1996; Kamiya et al., 2002), trnH-psbA spacer (Demesure et al., 1996; Okaura et al., 2007) and trnQ-trnS spacer (Kanno et al., 2004). Each reaction mixture contained 2x GoTaq Hot Start Colorless Master Mix (Promega, Madison, WI, USA) with 400 μM dNTPs, 4 mM MgCl₂, 10 ng of DNA, and 5 pmol of primers in a total volume of 25 μl. PCR conditions included initial denaturation at 94 ℃ for 2 min, followed by 35 cycles of denaturation at 94 ℃ for 1 min, annealing at 56 ℃ (51 ℃ for trnT-trnL spacer and 60 ℃ for trnQ-trnS spacer) for 1 min, extension at 72 ℃ for 2 min, and a final extension at 72 ℃ for 7 min using a GeneAmp 9700 Thermal Cycler (Applied Biosystems, Foster City, CA, USA) or a Veriti 96-Well Thermal Cycler (Applied Biosystems). PCR products were subsequently held at 4 ℃ until electrophoresis on a 0.8% agarose mini-gel to determine the concentration of the amplified DNA. Successfully amplified PCR products were purified using ExoSAP-IT (Affymetrix, Santa Clara, CA, USA).

Sequencing reactions were performed directly on the PCR products using the BigDye Terminator (version 1.1) Cycle Sequencing Kit (Applied Biosystems). Template DNA (50 ng) was subjected to labeling in a reaction cycle consisting of preheating for 2 min at 96 ℃ and 25 cycles of denaturation at 90 ℃ for 30 s, annealing at 50 ℃ for 30 s, and extension at 60 ℃ for 4 min. Labeled products were purified by ethanol precipitation and then applied to an ABI PRISM 310 Genetic Analyzer (Applied Biosystems). ATGC v.6.0 (GENETYX, Tokyo, Japan) was used to assemble the sequence data. Newly sequenced data were deposited in DDBJ (accession nos. LC275066 – LC275099 and LC275287 – LC275426 for matK, LC275100 – LC275138, LC275602 – LC275765 and LC368465 – LC368474 for trnQ-trnS, LC275244 – LC275286 for atpB-rbcL, LC275427 – LC275577 for psbA-trnH, and LC275578 – LC275601 for trnL-trnF).

Microsatellite genotyping

Thirty-two candidate primer pairs developed for other Quercus species were examined for amplification of the loci in Q. phillyraeoides (Supplementary Table S1). We selected the following eleven polymorphic loci: bcqm42, bcqm76 and bcqm94 (Mishima et al., 2006); QM58TGT and QM69-2M1 (Isagi and Suhandono, 1997); QpZAG16 and QpZAG36 (Steinkellner et al., 1997); QrZAG7, QrZAG11 and QrZAG112 (Kampfer et al., 1998); and MSQ13 (Dow et al., 1995). Others were not successfully amplified or did not show polymorphism and were not used. The primer pairs for amplifying the loci bcqm42, bcqm76, bcqm94, QM58TGT and QM69-2M1 were redesigned and customized for Q. phillyraeoides. The primer sequences used for this study are listed in Supplementary Table S2. PCR amplification was performed using the Type-it Microsatellite PCR Kit (Qiagen K. K., Tokyo). Fluorescent dye-labeled primers were prepared by custom order from Thermo Fisher Scientific Japan (Yokohama, Japan). The 5’-end of the forward primers of bcqm76-Qp, QM69-2M1-Qp, QrZAG112 and QpZAG16 was labeled with FAM, bcqm42-Qp, bcqm94-Qp and QpZAG36 was labeled with HEX, QM58TGT-Qp, MSQ13 and QrZAG11 was labeled with NED, and QrZAG7 was labeled with PET. Two to four different dye-labeled primer pairs were combined to make primer mix sets containing 0.2 μM of each in a total volume of 50 μl. The multiple PCR reaction mixture contained 3 μl of RNase-free water, 5 μl of 2x Type-it Multiplex PCR Master Mix, 1 μl of the primer mix set, and 1 μl of template DNA (200 ng for each PCR reaction). The PCR conditions included an initial denaturation at 95 ℃ for 5 min followed by 28–31 cycles of denaturation at 95 ℃ for 30 s, annealing at 49−60 ℃ for 90 s, and extension at 72 ℃ for 30 s with a final extension at 60 ℃ for 30 min. PCR was conducted with a GeneAMP PCR System 9700 Thermal Cycler (Applied Biosystems) or a Veriti 96-Well Thermal Cycler (Applied Biosystems). The concentration of the amplified DNA was determined by electrophoresis on a 0.8% agarose mini-gel. PCR products were diluted 200-fold and applied to an ABI PRISM 310 Genetic Analyzer. Fragment size was determined using GeneScan (Applied Biosystems) and genotyped using Gene Mapper v.3.0 (Applied Biosystems).

Analysis of cpDNA sequence variation

Nucleotide diversity (π) (Nei and Li, 1979) and scaled mutation rate (θ_w) (Watterson, 1975) were calculated for the five cpDNA regions using DnaSP v.5.10 (Librado and Rozas, 2009). Gaps were excluded from the calculation. Tajima’s D neutrality test (Tajima, 1989) was performed using DnaSP software. Statistical significance was determined by coalescent simulation implemented in the software. Sequences were manually assembled into haplotypes based primarily on nucleotide substitutions. We also corrected and redefined the haplotypes previously reported by Liu et al. (2013) because two haplotypes were misidentified by reading errors (i.e., haplotype B in Muroto-misaki and haplotype E in Ryugadake). To clarify the genetic relationships among the haplotypes, a haplotype network analysis was conducted using TCS v.1.21 (Clement et al., 2000) including the corresponding sequences of Q. mongolica var. crispula (Okaura et al., 2007) and Q. variabilis (Chen et al., 2012) as outgroups. Haplotype diversity and genetic differentiation were also calculated using DnaSP.

Analysis of genetic variation in microsatellite markers

Microsatellite variation was analyzed for 28 populations sampled across the distribution range of Q. phillyraeoides. The total number of alleles (Na), the number of private alleles (AP), and the observed and expected heterozygosities within populations (H_O and H_E) for each locus, each population, and overall populations were calculated using GenAlEx v.6.4.1 (Peakall and Smouse, 2006). Inbreeding coefficients (F) for each population and their statistical significance were determined using FSTAT v.2.9.3.2 software (Goudet, 2002). The allelic richness (AR) (El Mousadik and Petit, 1996; Petit et al., 1998) was also determined using FSTAT. Fixation indices in the total populations (F_IT), among populations (F_ST) and within populations (F_IS) (Weir and Cockerham, 1984), as well as G_ST and G’_ST as defined by Nei (1987), were also determined for each locus and overall loci using FSTAT. AMOVA (analysis of molecular variance) (Excoffier et al., 1992) was performed using GenAlEx v.6.4.1 to partition genetic variance into components attributable to the variation among groups, among populations, within populations and within individuals. Pairwise genetic distance matrices for the 28 populations were made for several distance measures using GenAlEx and were used for PCoA (principal coordinate analysis) (Gower, 1966) to visualize the pattern of genetic relationships among populations. A phylogenetic tree was also constructed based on Nei’s standardized genetic distance (D_A) (Nei et al., 1983) using POPTREE2 (Takezaki et al., 2010). The relationship between genetic distance and geographic distance was tested under the isolation by distance (IBD) hypothesis (Rousset, 1997) using the Mantel test. The possibility of a recent bottleneck was examined using BOTTLENECK v.1.2 (Piry et al., 1999). When rapid population size reduction (bottleneck) events occur, rare alleles are easily lost (Nei et al., 1975); this causes an excess of heterozygosity computed from a sample of genes (H₀) compared to expected heterozygosity at mutation-drift equilibrium (H₁), and hence provides criteria for testing a recent bottleneck (Cornuet and Luikart, 1996). Bottleneck significance was determined by the Wilcoxon test under the two-phase model (TPM) in which the default setting [30% SSM (stepping stone model) and 70% IAM (infinite allele model)] was used.

The model-based clustering method implemented in STRUCTURE v.2.3 (http://pritch.bsd.uchicago.edu/) was applied to infer population structure. The program identifies distinct genetic populations (clusters) and assigns individuals to the clusters. To estimate the optimum number of clusters K, ΔK, the rate of change in the log probability of data [LnP(D)] between successive K values, was calculated according to the method of Evanno et al. (2005) using Structure Harvester (Earl and vonHoldt, 2012). A computational geometrical approach implemented in BARRIER was conducted to determine the locations of genetic barriers, which prevent gene flow, and to construct boundaries of genetic change that occurred at an abrupt rate (Manni et al., 2004). The significance of barriers was evaluated by bootstrap analysis implemented in the program.

RESULTS

CpDNA variation

A total of 58 individual samples from 21 populations covering the entire distribution range of Q. phillyraeoides were sequenced for the five cpDNA regions, trnL-trnF spacer (444 bp), atpB-rbcL spacer (765 bp), matK gene (1,454 bp), trnH-psbA spacer (940 bp) and trnQ-trnS spacer (572 bp), covering a total aligned length of 4,175 bp (Table 1). A total of four nucleotide substitutions (two in the matK gene and two in the trnQ-trnS spacer) and two mononucleotide-repeat length variations (one in the trnH-psbA spacer and one in the trnQ-trnS spacer) were found (Table 2). No variation in the trnL-trnF and atpB-rbcL spacers was detected. π and θ_w across the five regions were 0.00016 and 0.00023, respectively. Tajima’s D was −0.6837 and was not significantly different from zero.

Table 2. Polymorphic sites in three cpDNA regions in Quercus phillyraeoides

Haplotype	Frequency	Polymorphic sites
		matK gene		trnH-psbA spacer	trnQ-trnS spacer
		177	1098	298	260	274	895
A	76.8	T	T	A11	G	A15	A
B	6.8	T	C	A11	G	A15	A
C	16.4	C	T	A9	T	A16	G

The polymorphic site positions are based on the standard sequences from the 5’ end of accessions AB060064 (matK gene), AB650455 (trnH-psbA spacer) and AB095323 (trnQ-trnS spacer).

Three haplotypes (A, B and C) were identified from the sequence variation detected in the five cpDNA regions (Table 2). The results showed that the matK gene sequence is diagnostic for defining haplotypes. Next, whole or partial sequences of the matK gene were sequenced for a total of 280 individual samples from all 41 populations to determine the haplotype. The haplotype frequency in each population is summarized in Table 1. Haplotype A was the most common (76.8%) and prevalent on the Pacific coast from the northeast population in the Boso Peninsula (IW) to Shikoku, and was also found in some populations in the Seto Inland Sea and Kyushu. Haplotype B was locally restricted (6.8%) and only found in three populations on the Kii Peninsula (TB, KS and SS). Haplotype C was the second most common (16.4%), but was found only in the western part of the distribution, namely the Seto Inland Sea, Kyushu and Okinawa. The global haplotype diversity (H_d) was 0.303 ± 0.00671 and the coefficient of genetic differentiation (G_ST) was 0.894. Geographic haplotype distribution is shown in Fig. 1. A haplotype network of the five cpDNA regions was reconstructed, including corresponding Qm1 sequences in Q. mongolica var. crispula (Okaura et al., 2007) and H4, H5 and H11 in Q. variabilis (Chen et al., 2012), and showed that haplotype A was ancestral to haplotypes B and C. Haplotype B was derived from A with only one nucleotide substitution; haplotype C was derived from haplotype A with three nucleotide substitutions and two mononucleotide-repeat length changes (Fig. 2).

Fig. 2.

The chloroplast DNA haplotype network constructed using TCS 1.21 (Clement et al., 2000). Haplotype Qm1 is from Quercus mongolica var. crispula (Okaura et al., 2007). Haplotypes H4, H5 and H11 are from Q. variabilis (Chen et al., 2012). Small circles indicate putative intervening haplotypes.

Microsatellite variation

Microsatellite genotypes were determined for 536 individual samples from 28 populations. In total, 104 alleles were detected for the 11 microsatellite loci (Table 3). Population genetic parameters were estimated using GenAlEx v6.4.1 (Peakall and Smouse, 2006) for 24 populations, in which IG, OM and OO populations were combined because of their close geographic proximity. NT and YK populations were eliminated from the calculation because of small sample sizes. The population genetic parameters estimated in each locus and in each population are listed in Table 3 and Table 4, respectively. None of the inbreeding coefficients were significantly larger than zero (P < 0.05) after the application of Bonferroni correction for type I errors (Rice, 1989), although some populations (OG, HZ and RY) were very close to the significance threshold (Table 4). The Weir and Cockerham (1984) estimator of F_IS for overall loci was significantly larger than zero at the 5%, but not at the 1%, level (Table 5).

Table 3. Summary of genetic variation and genetic differentiation estimated for each locus

Locus name	Na	HO	HS	HT	GIS	GST	G’ST	P
bcqm42	8	0.426	0.439	0.490	0.030	0.106	0.110	0.0001
bcqm76	6	0.521	0.682	0.739	0.237	0.077	0.080	0.0001
bcqm94	5	0.222	0.225	0.240	0.014	0.062	0.064	0.0001
QM58TG	8	0.658	0.628	0.676	−0.048	0.072	0.075	0.0001
QM69-2M1	13	0.501	0.555	0.611	0.097	0.092	0.096	0.0001
QpZAG16	14	0.481	0.679	0.776	0.291	0.125	0.130	0.0001
QpZAG36	7	0.516	0.483	0.543	−0.068	0.111	0.116	0.0001
QrZAG7	18	0.800	0.818	0.890	0.022	0.081	0.084	0.0001
QrZAG11	23	0.738	0.819	0.886	0.098	0.076	0.079	0.0001
QrZAG112	10	0.582	0.574	0.655	−0.014	0.124	0.129	0.0001
MSQ13	3	0.112	0.143	0.151	0.219	0.053	0.056	0.0001
Overall	104	0.505	0.550	0.605	0.081	0.092	0.096	0.0001

Na, observed number of alleles; H_O, observed heterozygosity; H_S, average heterozygosity; H_T, total heterozygosity; G_IS, inbreeding coefficient; G_ST, coefficient of genetic differentiation; G’_ST, coefficient of unbiased genetic differentiation; P, probability of G_ST being larger than zero. G_IS, G_ST and G’_ST are from Nei (1987). Tests are based on 10,000 randomizations using the exact G-test (Goudet et al., 1996).

Table 4. Summary of genetic variation estimated using the 11 microsatellite markers in each population and the result of Bottleneck test

Population	n	Na	AP	AR	HO	HE	uHE	F	Bottleneck testa
									A	B	C
IW	32	39	0	37.0	0.580	0.533	0.541	−0.077	0.999	0.002	0.003
NK	23	52	0	48.0	0.530	0.537	0.549	−0.013	0.959	0.051	0.102
OG	26	58	0	52.4	0.507	0.584	0.596	0.135	0.840	0.183	0.365
HZ	31	63	0	55.6	0.496	0.570	0.579	0.078	0.740	0.289	0.577
TB	17	59	2	58.1	0.545	0.586	0.604	0.072	0.650	0.382	0.765
KS	20	61	0	56.8	0.509	0.560	0.574	0.111	0.517	0.517	1.000
SR	23	56	0	51.5	0.569	0.571	0.583	0.002	0.861	0.160	0.320
MR	19	49	0	46.9	0.435	0.475	0.488	0.098	0.423	0.615	0.846
KR	16	45	0	45.0	0.494	0.499	0.515	0.016	0.875	0.150	0.301
SY	36	65	1	52.3	0.528	0.572	0.580	0.130	0.650	0.382	0.765
SG	25	48	1	43.0	0.462	0.500	0.510	0.064	0.711	0.319	0.638
TM	20	49	1	47.0	0.550	0.563	0.577	0.013	0.949	0.062	0.123
YS	19	49	0	47.5	0.526	0.540	0.555	0.027	0.935	0.080	0.160
NG	20	48	3	46.6	0.577	0.600	0.616	0.055	0.999	0.001	0.002
IG+OM+OO	22	60	3	56.5	0.558	0.567	0.580	0.024	0.880	0.139	0.278
SD	20	54	0	51.5	0.495	0.553	0.567	0.090	0.920	0.097	0.193
TS	19	52	2	49.8	0.512	0.527	0.542	0.018	0.754	0.278	0.557
AS	19	58	1	55.7	0.574	0.584	0.600	−0.015	0.966	0.042	0.083
US	19	50	2	47.9	0.488	0.501	0.515	0.006	0.551	0.483	0.966
ST	19	49	1	47.2	0.478	0.503	0.517	0.043	0.633	0.410	0.820
NM	19	59	1	56.5	0.464	0.550	0.564	0.131	0.483	0.551	0.966
KK	21	42	1	40.6	0.481	0.447	0.458	−0.101	0.674	0.367	0.734
RY	18	55	0	53.5	0.434	0.553	0.569	0.232	0.618	0.416	0.831
IZ	24	35	1	33.7	0.330	0.377	0.385	0.078	0.722	0.313	0.625
Mean	21.9	4.754	0.83	49.2	0.505	0.535	0.548	0.052
SE	0.29	0.167	0.19	1.26	0.015	0.014	0.014	0.014

n, sample size; Na, observed number of alleles; AP, number of private alleles; AR, allelic richness; H_O, observed heterozygosity; H_E, expected heterozygosity; uH_E, unbiased heterozygosity; F, inbreeding coefficient.

^a   Bottleneck was examined by Wilcoxon test: A, probability from one-tailed test for heterozygosity deficiency; B, probability from one-tailed test for heterozygosity excess; C, probability from two-tailed test for heterozygosity deficiency or excess.

Table 5. F-statistics estimated for overall loci based on the genetic variation revealed by microsatellite markers

	FIT	FST	FIS
Estimate	0.169	0.097	0.079
95% CI	0.105–0.238	0.084–0.112	0.011–0.152
99% CI	0.089–0.26	0.081–0.117	−0.006–0.174

F-statistics were estimated after Weir and Cockerham (1984). Confidence intervals (CIs) were determined by bootstrapping.

The expected heterozygosity, which Nei (1987) called gene diversity, as well as allelic richness, are commonly used to evaluate genetic diversity in a population as they are altered in response to changes in population size. Relative values of these two estimates across the 24 populations are illustrated in Fig. 3. Both values decreased in the southern KK and IZ populations. In the northeastern populations, both values were relatively high in HZ, TB and KS, but a sharp decline of allelic richness occurred toward the northeastern population of IW. Lowered allelic richness relative to heterozygosity was also observed in some of the populations in the Seto Inland Sea, such as TM, YS and NG. To infer a recent population size change in Q. phillyraeoides populations, BOTTLENECK software was applied for the microsatellite data. Under the TPM, three populations showed significant excess in observed heterozygosity (H₀) compared to expected heterozygosity (H₁) (IW, P < 0.05; NG, P < 0.01 and AS, P < 0.05) by the Wilcoxon test (Table 4). The NK population showed a nearly significant heterozygosity excess (P < 0.051). None of the populations showed significant deficiency in observed heterozygosity (Table 4).

Fig. 3.

The geographic distribution of two population genetics estimates, expected heterozygosity (H_E) and allelic richness (AR). Standard deviation of each estimate is used as the graphed unit and scaled on the y-axis, adjusting the grand means to zero.

Genetic structure revealed by microsatellite markers

The IBD hypothesis was examined to determine the correlation between genetic distance and geographic distance (Fig. 4). A significant correlation (R² = 0.3079, P < 0.001) was observed between pairwise geographic distance and genetic differentiation expressed as F_ST/(1 − F_ST) (Rousset, 1997) by the Mantel test. Global genetic differentiation was examined using F-statistics (Weir and Cockerham, 1984). Genetic differentiation among populations, F_ST, was 0.097 and significantly larger than zero (P < 0.01) (Table 5). The G_ST (Nei, 1987) was also large at 0.092 and significantly larger than zero (P < 0.001; Table 3). AMOVA showed that 10% of the total genetic variance was attributed to the variation among populations, 83% was attributed to the variation within individuals, and 7% to that among individuals. Cluster analysis was performed by constructing a NJ dendrogram based on D_A distance (Nei et al., 1983) (Fig. 5). Three groups were recognized: one includes populations on the Pacific coast in Kanto-Tokai, the Kii Peninsula and the Muroto-misaki Cape, another includes populations in the Bungo-suido Channel (the channel between Kyushu and Shikoku) area, and the third includes the populations in southern Kyushu and Okinawa. Populations in the Seto Inland Sea belonged either to the first or to the second groups. Close genetic relationships supported by high bootstrap value were observed between KS and RY, and between KR and SG, in the dendrogram. Similar grouping was revealed by PCoA based on a pairwise genetic distance matrix using Shannon’s index, in which some populations in the Seto Inland Sea (YS, OO and OM) were situated between the groups in the northeastern Pacific coast to the Muroto-misaki Cape and those in the Bungo-suido Channel area (Supplementary Fig. S1).

Fig. 4.

Correlation between pairwise geographic distance (GGD) and genetic differentiation shown by F_ST/(1-F_ST). Linear regression (R² = 0.30789) is highly significant (P < 0.001) by the Mantel test.

Fig. 5.

Dendrogram constructed using the NJ method based on Nei’s D_A distance in POPTREE2 (Takezaki et al., 2010). The underlined population codes show populations in the Seto Inland Sea. The dendrogram was drawn by rooting the midpoint. The branch lengths are proportional to the distance shown by the scale bar.

Population structure at an individual level was examined using the STRUCTURE program. The ΔK plot (Supplementary Fig. S2) suggests that the most likely number of genetic clusters was K = 2, and the second was K = 4. Bar plots of individuals with an estimated membership fraction for K = 2 and K = 4 are shown in Fig. 6. At K = 2, populations were divided into two groups according to the dominating clusters: one was dominated by cluster I and included populations in the northeastern Pacific coast from IW to MR in Shikoku, and the other was dominated by cluster II and included populations in the southwestern Pacific coast from AS in Shikoku to IZ in Okinawa (Fig. 6). In the Seto Inland Sea, populations dominated by clusters I and II were intermingled (Fig. 6). The populations YS in the Seto Inland Sea and NM in Kyushu were a mixture of the two clusters. A significant genetic differentiation between the two groups of populations dominated by cluster I and cluster II was revealed by AMOVA and accounted for 3% of the total genetic variance (Supplementary Table S3). At K = 4, clusters I and II were further divided into Ia and Ib, and IIa and IIb, respectively (Fig. 6). Three regional groups of populations were recognized based on the dominating clusters. The first one (group 1), which is the same as the group with cluster I at K = 2, included populations on the northeastern Pacific coast to the Muroto-misaki Cape and mostly admixtures of cluster Ia and Ib. The second one (group 2) is dominated by cluster IIa and included the populations in the Bungo-suido Channel area. The third one (group 3) is dominated by cluster IIb and included populations in southern Kyushu and Okinawa. These groups mostly correspond to the three groups revealed by the NJ dendrogram (Fig. 5), and also by the scatter plot of PCoA (Supplementary Fig. S1). A significant genetic differentiation among the three groups of populations was revealed by AMOVA, accounting for 6% of the total genetic variance (Supplementary Table S4).

Fig. 6.

Bar plots of STRUCTURE analysis results. Each bar shows an individual with an estimated membership fraction in each of the inferred number of clusters, K. Upper: bar plot for K = 2. Lower: bar plot for K = 4. Population codes are shown at the top. Letters at the bottom show the chloroplast DNA haplotypes observed in the population. At K = 2, the cluster of dark green is designated as cluster I, and that of pink as cluster II. The F_ST values for clusters I and II were 0.0709 ± 0.0004 (SE) and 0.0966 ± 0.0005 (SE), respectively. At K = 4, the clusters of dark green and pale green are designated as clusters Ia and Ib, respectively. Similarly, the clusters of pink and red are designated as cluster IIa and IIb, respectively. The F_ST values for clusters Ia, Ib, IIa and IIb were 0.1843 ± 0.0010 (SE), 0.1023 ± 0.0007 (SE), 0.1797 ± 0.0005 (SE) and 0.1910 ± 0.0005 (SE), respectively.

To infer genetic barriers on a geographic map, BARRIER analysis (Manni et al., 2004) was conducted on the microsatellite data (Supplementary Fig. S3). One barrier with a high bootstrap value appeared between IZ in Okinawa and the other populations. Secondary barriers were observed between island populations in the Seto Inland Sea. KK in southern Kyushu is also isolated from the others. No barriers were observed among the populations on the Pacific coast.

DISCUSSION

Genetic variation revealed by cpDNA sequences and microsatellite markers

The π value across the five cpDNA regions of Q. phillyraeoides [0.00016 ± 0.00005 (SE)] was relatively low compared to other oak species in Japan. The value for Japanese oak, Q. mongolica var. crispula, was 0.00055 ± 0.00001 (SE) in four cpDNA regions (Okaura et al., 2007), and that of Japanese beech, Fagus crenata, was 0.00220 ± 0.00026 (SE) for three cpDNA regions (Okaura and Harada, 2002). Although the examined loci were different, the genetic diversity of microsatellite markers was also lower in Q. phillyraeoides. Gene diversity (H_E) averaged over the 24 populations for 11 microsatellite markers in this study was 0.535 ± 0.0104 (SE); in contrast, gene diversity was 0.724 ± 0.0056 (SE) in Q. mongolica var. crispula (Ohsawa et al., 2011) and 0.839 ± 0.0050 (SE) in F. crenata (Hiraoka and Tomaru, 2009). The high genetic diversity in Japanese oak and Japanese beech reflects their endemicity to the Japanese archipelago, which has enabled them to accumulate more genetic variation. Like Q. phillyraeoides, the oriental oak Q. variabilis found in Japan originated in mainland China. The oriental oak populations in Japan were found to be fixed to one cpDNA haplotype, while 26 haplotypes were detected in mainland China (Chen et al., 2012). The migration process must have been accompanied by a series of founding events and caused considerable loss of genetic variation in Q. phillyraeoides, especially regarding cpDNA variation. The Japanese archipelago belonged to the same ancient landmass at the eastern margin of the Eurasian continent and became separated by the opening of the Japan Sea, commencing 15 Ma (million years ago) (Taira, 1990; Osozawa et al., 2012). However, the western part of Japan was still connected with mainland China until the opening of the East China Sea in the early Middle Pleistocene (Minato and Ijiri, 1976). The connected area could have been the main migration route of Q. phillyraeoides (Xie et al., 2011). Fossils of Stegodon orientalis, which represents the Middle Pleistocene fauna of Southern China, have been found in strata aged between 0.5 and 0.3 Ma in a broad area of Japan (Kamei et al., 1988). The latest time probably corresponds to the opening of the East China Sea up to the Korean Peninsula. The Ryukyu Arc, which extends from Kyushu to Taiwan in a span of 1,200 km, is another possible migration route of Q. phillyraeoides to Japan. This route, however, was blocked by the development of the Tokara gap north of the Amami islands and the Kerama gap south of the Okinawa islands around 1.5 Ma, associated with the backarc rifting of the Okinawa Trough (Osozawa et al., 2012). After the opening of the Japan Sea and the East China Sea, immigration of animals and plants from the Eurasian continent was limited to land bridges formed during the glacial ages in the late Pleistocene in the north through Sakhalin and in the south through the Korean Peninsula (Taira, 1990). A migration of Q. phillyraeoides via a land bridge from the Korean Peninsula is unlikely, however, because the oak occurs neither in the Korean Peninsula nor in northern Kyushu. These considerations suggest that migration of Q. phillyraeoides from mainland China came to an end when the East China Sea opened about 0.3 Ma.

Both allelic richness and heterozygosity in microsatellite markers were remarkably reduced in the southernmost island populations of IZ and KK (Fig. 3). There was no evidence of a recent bottleneck in these populations (Table 4). Together with the observed monotypic pattern of the STRUCTURE diagram, this suggests that these populations have been isolated for a long time, reaching a genetic equilibrium. The water depth between KK and the west coast of Kyushu is at most 100 m, and, hence, these regions must have been connected at the LGM approximately 20,000 years ago. Isolation could have occurred through separation from Kyushu by rising seawater levels. The IZ population could be a remnant population that migrated from mainland China across the land bridge of the Ryukyu Arc. Because we have no record of Q. phillyraeoides in the Ryukyu Islands except in Okinawa-jima, including Izena-jima (Hatsushima and Amano, 1994), further study to investigate the origin of IZ, including a possible ancient human transplantation, is needed.

Genetic structure revealed by cpDNA haplotype distribution and microsatellite STRUCTURE analysis

cpDNA haplotypes mostly co-occurred with certain clusters revealed by STRUCTURE analysis of microsatellite markers. Haplotypes A and C associated with clusters I and II, respectively, at K = 2 (Fig. 6). This is obvious in the mixed populations of YS and NM. In these populations, individuals dominated by cluster I had haplotype A, while individuals dominated by cluster II had haplotype C (data not shown). This implies that haplotype A is essentially fixed in cluster I, and haplotype C is fixed in cluster II. The larger F_ST value in cluster II than in cluster I revealed by STRUCTURE analysis (Fig. 6 legend) implies that cluster II experienced more extensive genetic drift. The differentiated pattern in the geographical distribution of cluster I with haplotype A and cluster II with haplotype C suggests that they independently migrated to Japan as two different lineages. Populations AS and KK are exceptional as they are dominated by cluster II but fixed with haplotype A. This probably occurred by organelle genome capture when the two types grew together in southern Kyushu. A recent study has shown that introgression is likely to take place preferentially from the resident species toward the invading species. It was also demonstrated by simulation that introgression is stronger for genes experiencing little gene flow, such as organelle genes (Currat et al., 2008). This suggests that cluster I with haplotype A was resident in this area and was invaded by cluster II with haplotype C. Haplotype B may have arisen from a de novo mutation in the group dominated by cluster I in the southern part of the Kii Peninsula.

At K = 4, STRUCTURE analysis revealed three regional groups with differentiated genetic backgrounds. Many populations in group 1 are admixtures of clusters Ia and Ib, indicating the occurrence of frequent gene flow between them as shown by the absence of barriers on the northeastern Pacific coast (Supplementary Fig. S3). The F_ST values revealed by the STRUCTURE analysis (Fig. 6 legend) are similar for clusters Ia, IIa and IIb, but much smaller in cluster Ib. This suggests that the population size of cluster Ib remained fairly large even during the last glacial age. Cluster Ib dominates the populations in the Kii Peninsula (TB, KS and SR) including MR in the Muroto-misaki Cape. Cluster Ib probably existed as a refugial population on the coast of the Kii Peninsula and expanded its distribution both eastward and westward after the LGM. Another refugium with cluster Ia probably existed in the Boso Peninsula and expanded westward. The existence of multiple refugia of lucidophyllous forests in the LGM has been suggested at the southern ends of major peninsulas or promontories on the Pacific side of Honshu, Shikoku and Kyushu (Kamei and Research Group for the Biogeography from Würm Glacial, 1981; Hattori, 2002; Aoki et al., 2004).

In the Seto Inland Sea, populations SY, TM, YS, NG, IG, OM and OO shared a common genetic background of cluster Ib, whereas KR, SG, NT, SD and TS had the genetic background of cluster IIa. At the LGM, the Seto Inland Sea basin was completely dried (Yonekura et al., 2001). The present ecological preference of Q. phillyraeoides suggests the absence of this species in the dried basin. Two waterways had formed on both sides of the Bisan-seto south of the Kojima Peninsula as a divide (Fig. 1, g). One flowed eastward and into the Pacific Ocean through the present-day Kii-suido Channel, while the other flowed westward and into the Pacific Ocean through the present-day Bungo-suido Channel (Yonekura et al., 2001). The populations in the Seto Inland Sea with a genetic background of cluster Ib might have been derived from refugia in the Kii Peninsula following the opening of the Kii-suido after the LGM as a result of rising sea levels. Similarly, the populations with a genetic background of cluster IIa are likely to have been derived from refugia located in southern Kyushu when the Bungo-suido Channel opened. These two channels were connected at the position of Bisan-seto 8,500 years ago (Ohta et al., 2004). The barriers observed between the coastal and island populations in this area (Supplementary Fig. S3) suggest that the populations derived via the Bungo-suido and Kii-suido Channels reached Bisan-seto around this time and have since been isolated by rising sea levels. Some populations with a genetic background of cluster Ib (NG, IG, OM, OO) and those with cluster IIa (KR and SG) might have migrated further to the west and east of Bisan-seto, respectively, after the two channels connected.

Both cluster IIa and IIb occurred in southern Kyushu (Fig. 6). Based on palynological studies, Tsukada (1983) suggested that a major refugium of lucidophyllous forest was located in the paleo-Yaku Peninsula at the LGM, which included the present-day Osumi Peninsula, Satsuma Peninsula, Yakushima, and Tanegashima. A refugium of cluster II was probably located in the paleo-Yaku Peninsula at the LGM and split into clusters IIa and IIb afterward. Mixing of cluster IIa and IIb was observed in NM in the Osumi Peninsula and YK in Yakushima (Fig. 6). This suggests that splitting occurred before Yakushima (and Tanegashima) became separated from the southern part of Kyushu.

A close phylogenetic similarity was observed between RY and KS in the NJ dendrogram (Fig. 5). The RY population was most likely transplanted by humans. An example of human transplantation of Q. phillyraeoides by fishermen in the Kii Peninsula in the 17^th century was recorded on the Goto islands in Nagasaki prefecture near the Amakusa islands, where the RY population is located (http://shinkamigoto.nagasaki-tabinet.com/guide/52263).

Implications of the study for conservation of Q.phillyraeoides

Along the Pacific coast, Q. phillyraeoides populations are divided into several genetically differentiated groups, and may have been derived from multiple independent refugia at the LGM. These regional groups with genetically differentiated backgrounds should be considered as conservation units. Two very specialized populations, KK and IZ, are likely to have considerable conservation value. KK is fairly monomorphic due to long-term isolation, and forms pure dense stands along a 4 km-long sand reef (Nagamenohama), constructing a very specific landscape. IZ has been retained in very small numbers and is seriously endangered (http://www.pref.okinawa.jp/site/kankyo/shizen/hogo/). Hence, this population is valuable both for conservation and for studying the migration process of Q. phillyraeoides from mainland China through the Ryukyu Arc. Transplantation or replantation of Q. phillyraeoides to foster coastal protection forests and next-generation charcoal production forests should be done based on regional groups as management units, which share a common genetic background considering their genetic variation and differentiation.

ACKNOWLEDGMENTS

We would like to thank the many people who assisted with sample collection and provided important local information, especially Dr. Masatoshi Hara, Chiba Prefecture Museum, Dr. Tatuo Sweda, Faculty of Agriculture, Ehime University, Mr. Hideyuki Katai, Shizuoka Prefecture, and Ms. Akiko Uchida, Oita Prefecture. This study was supported in part by a Grant-in-Aid from the Japan Society for the Promotion of Science, no. 22380086.

REFERENCES

•  Aoki,  K.,  Suzuki,  T.,  Hsu,  T. W., and  Murakami,  N. (2004) Phylogeography of the component species of broad-leaved evergreen forests in Japan, based on chloroplast DNA variation. J. Plant Res. 117, 77–94.
•  Balloux,  F., and  Goudet,  J. (2002) Statistical properties of population differentiation estimators under stepwise mutation in a finite island model. Mol. Ecol. 11, 771–783.
•  Chen,  D.,  Zhang,  X.,  Kang,  H.,  Sun,  X.,  Yin,  S.,  Du,  H.,  Yamanaka,  N.,  Gapare,  W.,  Wu,  H., and  Liu,  C. (2012) Phylogeography of Quercus variabilis based on chloroplast DNA sequence in East Asia: multiple glacial refugia and mainland-migrated island populations. PLoS One 7, e47268.
•  Clement,  M.,  Posada,  D., and  Creandll,  K. A. (2000) TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659.
•  Comes,  H. P., and  Kadereit,  J. W. (1998) The effect of Quaternary climatic changes on plant distribution and evolution. Trends Plant. Sci. 3, 432–438.
• Consortium of Japan Bincho-tan meeting (2011) Bulletin (in Japanese).
•  Cornuet,  J. M., and  Luikart,  G. (1996) Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144, 2001–2014.
•  Currat,  M.,  Ruedi,  M.,  Petiti,  R. J., and  Excoffier,  L. (2008) The hidden side of invasions: massive introgression by local genes. Evolution 62, 1908–1920.
•  Demesure,  B.,  Comps,  B., and  Petit,  R. J. (1996) Chloroplast DNA phylogeography of the common beech (Fagus sylvatica L.) in Europe. Evolution 50, 2515–2520.
•  Denk,  T., and  Grimm,  G. W. (2009) Significance of pollen characteristics for infrageneric classification and phylogeny in Quercus (Fagaceae). Int. J. Plant Sci. 170, 926–940.
•  Dow,  B. D.,  Ashley,  M. V., and  Howe,  H. F. (1995) Characterization of highly variable (GA/CT)_n microsatellites in the bur oak, Quercus macrocarpa. Theor. Appl. Genet. 91, 137–141.
•  Earl,  D. A., and  vonHoldt,  B. M. (2012) Structure Harvester: a website and program for visualizing structure output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361.
•  El Mousadik,  A., and  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, 832–839.
•  Evanno,  G.,  Regnaut,  S., and  Goudet,  J. (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620.
•  Excoffier,  L.,  Smouse,  P. E., and  Quattro,  J. M. (1992) Analysis of molecular variance Inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131, 479–491.
• Forest and Forestry Statistics Handbook (2016) Forestry Agency, the Ministry of Agriculture, Forestry and Fishery of Japan (in Japanese).
•  Goudet,  J. (2002) FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3.2). http://www2.unil.ch/popgen/softwares/fstat.html. Accessed 19 April 2017.
•  Goudet,  J.,  Raymond,  M.,  de Meeüs,  T., and  Rousset,  F. (1996) Testing differentiation in diploid populations. Genetics 144, 1933–1940.
•  Gower,  J. C. (1966) Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika 53, 325–338.
•  Hatsushima,  Y., and  Amano,  T. (1994) Flora of the Ryukyus, South of Amami Islands. The Biological Society of Okinawa, Nishihara, Okinawa (in Japanese).
•  Hattori,  T. (2002) Conservation of lucidophyllous forest from the viewpoint of phytogeography. In: Biology of Conservation and Restoration. (eds.: Yahara, T., Kawakubo, N.), pp. 203–222. The Society for the Study of Species Biology, Bun-ichi Sogo Shuppan, Tokyo (in Japanese).
•  Hewitt,  G. M. (2000) The genetic legacy of the Quaternary ice ages. Nature 405, 907–913.
•  Hiraoka,  K., and  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, 269–282.
•  Holderegger,  R.,  Kamm,  U., and  Gugerli,  F. (2006) Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc. Ecol. 21, 797–807.
•  Horikawa,  Y. (1972) Atlas of the Japanese Flora, Vol. 1. Gakken, Tokyo.
•  Isagi,  Y., and  Suhandono,  S. (1997) PCR primers amplifying microsatellite loci of Quercus myrsinifolia Blume and their conservation between oak species. Mol. Ecol. 6, 897–899.
•  Kamei,  S., and Research Group for the Biogeography from Würm Glacial (1981) Fauna and flora of the Japanese Islands in the last glacial time. Quat. Res. (Daiyonki-Kenkyu) 20, 191–205 (in Japanese with English abstract).
•  Kamei,  T.,  Kawamura,  Y., and  Taruno,  H. (1988) Mammalian stratigraphy of the Late Neogene and Quaternary in the Japanese Islands. Mem. Geol. Soc. Japan. 30, 181–204 (in Japanese with English abstract).
•  Kamiya,  K.,  Harada,  K.,  Mahani,  M. C., and  Latiff,  A. M. (2002) Genetic variation of Trigonobalanus verticillata, a primitive species of Fagaceae, in Malaysia revealed by chloroplast sequences and AFLP markers. Genes Genet. Syst. 77, 177–186.
•  Kampfer,  S.,  Lexer,  C.,  Glössl,  J., and  Steinkellner,  H. (1998) Characterization of (GA)_n microsatellite loci from Quercus robur. Hereditas 129, 183–186.
•  Kanno,  M.,  Yokoyama,  J.,  Suyama,  Y.,  Ohyama,  M.,  Itoh,  T., and  Suzuki,  M. (2004) Geographical distribution of two haplotypes of chloroplast DNA in four oak species (Quercus) in Japan. J. Plant. Res. 117, 311–317.
•  Kira,  T. (1948) On the altitudinal arrangement of climatic zones in Japan. Kanchi-Nogaku 2, 143–173 (in Japanese).
•  Kira,  T. (1977) A climatological interpretation of Japanese vegetation zones. In: Vegetation Science and Environmental Protection. (eds.: Miyawaki, A., Tüxen, R.), pp 21–30. Maruzen, Tokyo.
•  Kitamura,  S. (1957) Distribution of Plants, Genshoku Nihon Shokubutsu Zukan. Vol 1. Hoikusha, Osaka (in Japanese).
•  Li,  W.-H. (1997) Molecular Evolution. Sinauer Associates, Massachusetts, USA.
•  Librado,  P., and  Rozas,  J. (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452.
•  Liu,  H.-Z.,  Takeichi,  Y.,  Kamiya,  K., and  Harada,  K. (2013) Phylogeography of Quercus phillyraeoides (Fagaceae) in Japan as revealed by chloroplast DNA variation. J. For. Res. 18, 361–370.
•  Manel,  S.,  Schwartz,  M. K.,  Luikart,  G., and  Taberlet,  P. (2003) Landscape genetics: combining landscape ecology and population genetics. Trends Ecol. Evol. 18, 189–197.
•  Manni,  F.,  Guérard,  E., and  Heyer,  E. (2004) Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmonier’s algorithm. Hum. Biol. 76, 173–190.
•  Minato,  M., and  Ijiri,  S. (1976) The Japanese Archipelago, 3rd ed. Iwanami Shoten, Tokyo (in Japanese).
•  Mishima,  K.,  Watanabe,  A.,  Isoda,  K.,  Ubukata,  M., and  Takata,  K. (2006) Isolation and characterization of microsatellite loci from Quercus mongolica var. crispula. Mol. Ecol. Notes 6, 695–697.
•  Miyawaki,  A.,  Okuda,  S., and  Mochizuki,  R. (1983) Handbook of Japanese Vegetation, Revised Ed. Shibundo, Tokyo (in Japanese).
•  Nei,  M. (1987) Molecular Evolutionary Genetics. Columbia University Press, New York.
•  Nei,  M., and  Li,  W.-H. (1979) Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc. Natl. Acad. Sci. USA 76, 5269–5273.
•  Nei,  M.,  Maruyama,  T., and  Chakraborty,  R. (1975) The bottleneck effect and genetic variability in populations. Evolution 29, 1–10.
•  Nei,  M.,  Tajima,  F., and  Tateno,  Y. (1983) Accuracy of estimated phylogenetic trees from molecular data. II. Gene frequency data. J. Mol. Evol. 19, 153–170.
•  Ohsawa.  T.,  Tsuda,  Y.,  Saito,  Y., and  Ide,  Y. (2011) The genetic structure of Quercus crispula in northeastern Japan as revealed by nuclear simple sequence repeat loci. J. Plant Res. 124, 645–654.
•  Ohta,  Y.,  Naruse,  T.,  Tanaka,  S., and  Okada,  A. (2004) Regional Geomorphology of the Japanese Islands. Vol. 6 Geomorphology of Kinki, Chugoku, and Shikoku. University of Tokyo Press, Tokyo (in Japanese).
•  Okaura,  T., and  Harada,  K. (2002) Phylogeographical structure revealed by chloroplast DNA variation in Japanese Beech (Fagus crenata Blume). Heredity (Edinb) 88, 322–329.
•  Okaura,  T.,  Quang,  N. D.,  Ubukata,  M., and  Harada,  K. (2007) Phylogeographic structure and late Quaternary population history of the Japanese oak Quercus mongolica var. crispula and related species revealed by chloroplast DNA variation. Genes Genet. Syst. 82, 465–477.
•  Osozawa,  S.,  Shinjo,  R.,  Armid,  A.,  Watanabe,  Y.,  Horiguchi,  T., and  Wakabayashi,  J. (2012) Palaeogeographic reconstruction of the 1.55 Ma synchronous isolation of the Ryukyu Islands, Japan, and Taiwan and inflow of the Kuroshio warm current. Int. Geol. Rev. 54, 1369–1388.
•  Peakall,  R., and  Smouse,  P. E. (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol. Ecol. Notes 6, 288–295.
•  Petit,  R. J.,  El Mousadik,  A., and  Pons,  O. (1998) Identifying populations for conservation on the basis of genetic markers. Conserv. Biol. 12, 844–855.
•  Piry,  S.,  Luikart,  G., and  Cornuet,  J.-M. (1999) BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. J. Hered. 90, 502–503.
•  Pu,  C.-X.,  Zhou,  S.-K., and  Luo,  Y. (2002) A cladistics analysis of Quercus (Fagaceae) in China based on leaf epidermis and architecture. Acta Botanica Yunnanica 24, 689–698.
•  Rice,  W. R. (1989) Analyzing tables of statistical tests. Evolution 43, 223–225.
•  Rousset,  F. (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145, 1219–1228.
•  Steinkellner,  H.,  Fluch,  S.,  Turetschek,  E.,  Lexer,  C.,  Streiff,  R.,  Kremer,  A.,  Burg,  K., and  Glössl,  J. (1997) Identification and characterization of (GA/CT)_n-microsatellite loci from Quercus petraea. Plant Mol. Biol. 33, 1093–1096.
•  Taberlet,  P.,  Gielly,  L.,  Pautou,  G., and  Bouvet,  J. (1991) Universal primers for amplification of three non-coding regions of chloroplast DNA. Plant Mol. Biol. 17, 1105–1109.
•  Taira,  A. (1990) Birth of the Japanese Archipelago. Iwanami Shoten, Tokyo (in Japanese).
•  Tajima,  F. (1989) Statistical method for testing neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.
•  Takezaki,  N.,  Nei,  M., and  Tamura,  K. (2010) POPTREE2: Software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. Mol. Biol. Evol. 27, 747–752.
•  Terachi,  T. (1993) Structural alterations of chloroplast genome and their significance to higher plant evolution. Bulletin of the Institute for National Land Utilization Development, Kyoto Sangyo University 14, 138–148 (in Japanese).
•  Tsukada,  M. (1983) Vegetation and climate during the last glacial maximum in Japan. Quat. Res. 19, 212–235.
•  Tüxen,  R. (1956) Die heutige potentielle natürliche Vegetation als Gegenstand der Vegetatonskartierung. Angew. Pflanzensoziol. 13, 5–42 (in German).
•  Watterson,  G. A. (1975) On the number of segregating sites in genetical models without recombination. Theor. Popul. Biol. 7, 256–276.
•  Weir,  B. S., and  Cockerham,  C. C. (1984) Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370.
•  Xie,  C.-P.,  Yan,  F., and  Fang,  Y.-M. (2011) Geographical distribution of Quercus phillyraeoides A. Gray. Trop. Geogr. 31, 8–13 (in Chinese with English abstract).
•  Yang,  J. Vázquez,  L.,  Chen,  X.,  Li,  H.,  Zhang,  H.,  Liu,  Z., and  Zhao,  G. (2017) Development of chloroplast and nuclear DNA markers for Chinese oaks (Quercus subgenus Quercus) and assessment of their utility as DNA barcodes. Front. Plant Sci. 8, 816.
•  Yonekura,  N.,  Kaizuka,  S.,  Nogami,  M., and  Chinzei,  K. (2001) Regional Geomorphology of the Japanese Islands. Vol. 1 Introduction to Japanese Geomorphology. University of Tokyo Press, Tokyo (in Japanese).