Genetic Diversity and Structure of Apomictic and Sexually Reproducing Lindera Species (Lauraceae) in Japan

Abstract

Research Highlights: genetic diversity in populations were compared among related shrub species with different reproductive systems. Background and Objectives: Lindera species are dioecious trees or shrubs that produce seeds by mating of males and females. To evaluate the importance of genetic diversity for the persistence of natural populations, we compared genetic information among four Lindera species in Japan. Three are dioecious shrubs (Lindera praecox, Lindera umbellata, and Lindera obtusiloba) that produce seeds by sexual reproduction. The remaining species, Lindera glauca, reproduces by apomixis; only female plants are found in Japan. Materials and Methods: all four species were sampled across a wide geographic area, from Tohoku to Kyushu, Japan. Single nucleotide polymorphisms (SNPs) were detected by multiplexed ISSR genotyping by sequencing (MIG-seq) and the resulting genetic diversity parameters were compared among populations. Results: in all sexually reproducing species, the values of observed heterozygosity were close to the expected ones and the inbreeding coefficients were nearly 0. These results were supposed to be caused by their obligate outcrossing. The genetic difference increased, in ascending order, between a mother plant and its seeds, within populations, and across geographic space. We observed a substantial geographic component in the genetic structure of these species. For L. glauca, the genetic difference between a mother and its seeds, within populations, and across space were not significantly different from what would be expected from PCR errors. Genetic diversity within and among populations of L. glauca was extremely low. Conclusions: apomixis has the advantage of being able to found populations from a single individual, without mating, which may outweigh the disadvantages associated with the extremely low genetic diversity of L. glauca. This may explain why this species is so widely distributed in Japan. Provided that the current genotypes remain suited to environmental conditions, L. glauca may not be constrained by its limited genetic diversity.

1. Introduction

Genetic diversity allows for environmental adaptation [1,2]. In plants, genetic diversity depends in part on the reproductive system of a given species; outcrossing species tend to be more genetically diverse [3]. Genetic diversity has been found to be higher in dioecious species, which are obligate outcrossers, compared to related cosexual species [4]. However, dioecious plants may have two reproductive disadvantages: pollination limitation because two individuals of opposite sexes are required to reproduce sexually and low colonization ability, because only half of the population contributes to seed dispersal. The seed-shadow handicap, related to reduced dispersal around female plants, typically results in a clumped distribution at the local scale [5], and potentially also at greater geographic scales. These advantages and disadvantages are reversed in apomictic plants.

Apomixis, the asexual formation of seeds, occurs in over 400 angiosperm genera [6]. Apomixis is considered to be advantageous for population growth and expansion, because populations can be founded by a single individual and there are no pollination limitations [7]. However, because apomictic seeds are produced without meiosis or fertilization, each seed is a genetic replica of the maternal plant unless somatic mutation occurs. Therefore, apomictic plants may have reduced disease resistance and adaptability to environmental changes.

Lindera species (Lauraceae) are dioecious trees or shrubs that produce seeds by mating of males and females. There are over 100 species of Lindera in the world, seven of which are native to Japan [8]. Of Japanese Lindera species, Lindera glauca Bl. occurs in Japan as strictly female plants that reproduce by apomixis. There are remaining questions as to how populations of this species are maintained in Japan. If genetic diversity is exceptionally low, this may support the commonly accepted theory that apomixis is advantageous in the short term due to enhanced colonization ability. However, if L. glauca has a similar level of genetic diversity to sexually reproducing species, this would suggest that the trade-off for enhanced colonization is too high, and that there is high pressure to maintain genetic diversity even in apomictic species. Previous research on Taraxacum officinale [9] and Miconia albicans [7] suggested that different genotypes were maintained in apomictic populations through somatic mutation or meiosis.

The purpose of this study was to evaluate the importance of genetic diversity for the persistence of wild plant populations. For the purpose, we compared the genetic diversity and geographic genetic structure of L. glauca to three sexually reproducing congeners, Lindera praecox Bl., Lindera umbellata Thunb., and Lindera obtusiloba Bl., to determine the extent of variation in terms of genotype coincidence between mother plants and seeds, within populations and across geographic space.

2. Materials and Methods

2.1. Focal Species

All four species are dioecious, deciduous forest understory shrubs. Lindera praecox is endemic to Japan, while the remaining three species are found on the Asian Continent and the Japanese Archipelago. Pliocene fossils of L. praecox, L. umbellata, and L. obtusiloba have been found in Japan [10], and Pleistocene fossils of L. glauca have been found [11]. The species flower in spring, from late February to early May, with each species flowering for 2–4 weeks [8,12]. The flowers, which are borne in pseudo-umbel-shaped inflorescences, are small (4–9 mm in diameter), dish-shaped, yellow or green in color, and unisexual, with sterile rudiments of the other sex [12].

In the past, the seeds of L. praecox were used to press oil for lanterns [13]. Lindera umbellata extracts had the effects of prevention of onset and progression of nephropathy [14]. Lindera obtusiloba has been widely used in a traditional medicine as the anti-allergic inflammatory effect [15]. The fruits of L. glauca contain some aromatic essential oils, alkaloids, terpenoids, and flavonoids, which are widely employed in the essence, medicine, and chemical industries [16].

The sex ratio of populations of L. praecox, L. umbellata, and L. obtusiloba is approximately 1:1 in forests in western Japan [12]. Both male and female flowers attract insect pollinators, primarily Diptera and Coleoptera [8]. Regarding L. glauca, only female plants are found in Japan, but both males and females are found on the Asian Continent [17,18]. The sex ratio of L. glauca in the wild in mainland China is estimated to be predominantly female [19]. Dupont [17] conducted bagging experiments to clarify whether or not pollination was necessary for seed set of Lindera species. She covered flower buds on branches of female plants of L. umbellata, L. obtusiloba, and L. glauca with nylon meshes to keep pollen out. The results for L. umbellata and L. obtusiloba showed that pollination was essential for their seed set. On the other hand, females of L. glauca produced seeds without pollination or male flowers, indicating apomixis. For some woody plants, clonal growth arises by wide-ranging roots that produce new shoots, termed stolon. However, the four Lindera species do not reproduce laterally by stolons. Hence, stems derived from an individual can easily identified.

2.2. Field Collections and DNA Extraction

We collected leaves from 105 individuals of L. glauca in natural vegetations at 29 locations, representing the entirety of its distribution in Japan (Supplementary Material Table S1). We tried to collect up to five individuals from one location, but in some locations, we could only find one. As a result, leaves were taken from one to five individuals per location. From six locations among the 29 locations, i.e., Takaragaike, Kasugayama, Myokensan, Shigisan, Okayama-shi, and Umamiyama, five seeds from each selected mature individual were taken (Supplementary Material Table S2). For L. praecox, we collected leaves from 45 individuals (one to five individuals per location) in natural vegetations at 11 locations (Table S1). From two locations, among the 11 locations, i.e., Mitarai Valley and Tsuno, six to twelve seeds were collected from each of six mature individuals (Table S2). For L. umbellata, we collected leaves from 27 individuals at six locations on Honshu Island (two–five individuals per location) (Table S1). From a mature individual among six locations, i.e., Mitarai Valley, five seeds were collected (Table S2). For L. obtusiloba, we collected leaves from 39 individuals at eight locations on Honshu, Shikoku, and Tsushima Islands (four–five individuals per location) (Table S1). From two individuals at Mitarai Valley, 11 and 12 seeds were collected (Table S2). We recorded the geographic location, i.e., latitude and longitude, and diameter at breast height (DBH; 1.3 m) of each sampled individual (Table S1).

DNA was extracted from leaves and seed cotyledons for each species using a DNA extraction kit (DNeasy Plant Mini Kit; QIAGEN, Tokyo, Japan) or a DNA extraction buffer (UniversAll™ Extraction Buffer II; Yeastern Biotech, New Taipei City, Taiwan).

2.3. Multiplexed ISSR Genotyping by Sequencing (MIG-Seq) Analysis

DNA sequences were subjected to MIG-seq analysis [20] and single nucleotide polymorphisms (SNPs) were identified. In MIG-seq analysis, the region between single sequence repeats (SSRs), i.e., the inter-SSR (ISSR), is amplified by polymerase chain reaction (PCR). Approximately 1.5 ng/µL of DNA was used for the first PCR, as template DNA. The concentration of the DNA solution was measured using a Quantus^TM Fluorometer (Promega, Madison, WI, USA), and the concentration of the template DNA solution was adjusted by adding sterile water.

For the first PCR, a mixture of primers was prepared from eight types of forward and reverse primers intended to amplify the ISSR region. The fragments were amplified with a Multiplex PCR Assay Kit Ver. 2 (TaKaRa Bio Inc., Shiga, Japan) using 7 μL reaction volumes containing 1.0 μL of template DNA, 1.4 μL of PCR primers (mix of primers, each at a concentration of 2.0 μM), 3.5 μL of 2 × Multiplex PCR Buffer, 0.035 μL of Multiplex PCR Enzyme Mix and 1.065 μL of water. The PCR cycling conditions used for the Verti 96-Well Thermal Cycler (Applied Biosystems, Foster City, CA, USA) were as follows: 94 °C for 1 min followed by 27 cycles of 94 °C for 30 s, 38 °C for 1 min, 72 °C for 1 min, and finally 72 °C for 10 min.

The second PCR used the products of the first PCR diluted to a concentration of 1:50 by adding sterile water. The second PCR was performed to add complementary sequences as an index to both ends of the first PCR products to identify individuals. The fragments were amplified with PrimeSTAR GXL DNA polymerase (TaKaRa Bio Inc., Shiga, Japan) using 12 μL reaction volumes containing 2.5 μL of diluted first PCR products as template DNA, 2.4 μL of second PCR primers (mix of primers, each at a concentration of 2.0 μM), 2.4 μL of 5 × PrimeSTAR Buffer, 0.96 μL of dNTP mixture, 0.24 μL of PrimeSTAR GXL polymerase, and 3.5 μL of water. The PCR cycling conditions were as follows: 12 cycles of 98 °C for 10 s, 54 °C for 15 s, and 68 °C for 1 min. After purification of the second PCR products with AMPure XP magnetic beads, only 300–800 bp fragments with a size suitable for sequencing were selected by the Blue Pippin DNA Size Selection System (Sage Science, Beverly, MA, USA). Tape Station 4200 (Agilent Technologies, Santa Clara, CA, USA) and Qubit fluorometer (Life Technology, Carlsbad, CA, USA) were used to measure the frequency distribution of fragment lengths and the amount of fragments. Prepared to a concentration of 12 pM, the DNA library was used for sequencing on an Illumina MiSeq sequencing platform using a MiSeq Reagent Kit v3 (Illumina, San Diego, CA, USA). The sequencing results were input into Stacks v.1.47 [21] to identify loci and detect SNPs. Following Catchen et al. [22], we used “ustacks” for building loci, “cstacks” for creating the loci catalogue for all samples, and “sstacks” for checking against the catalog. Finally, the “population” module was used to estimate a “population genetics statistic” from a population of individual samples.

2.4. Data Analyses

Only loci with more than 50% of the population and less than 60% heterozygosity were extracted. For each species, the number of alleles (Na), number of effective alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), and an inbreeding coefficient (F_IS) were estimated using GenAlEx v.6.5 [23]. Genotypes, including those from seed samples, were compared to identify clones and to assess genotype diversity for each species.

We calculated the proportion of different alleles on each SNP loci as follows: if individuals shared an identical genotype (e.g., A/T and A/T), the proportion of different alleles was 0; if one allele differed between individuals (e.g., A/T and A/G), it was 0.5; and if the genotypes differed completely (e.g., A/T and G/C), it was 1. Once all loci had been assessed, the mean proportion of different alleles was defined as the genetic difference between individuals. To investigate the extent to which genetic variation occurs in apomictic seed production relative to sexual reproduction, the genetic difference between mother plants and their seeds was determined. We then determined the genetic difference between grown-up individuals to assess genotype diversity within a population by pairings limited to individuals growing in the same location. For the assessment the genotype diversity across the entire study area, the genetic difference was calculated for all pairs of grown-up individuals after combining all individuals sampled from all locations into a single pool. We estimated the possibility that pseudo-variant loci occurred originating from PCR errors using the genotype data from pairs of replicates of the same individual. We performed PCR amplification and genotyping twice for seven individuals of L. glauca and six individuals of L. praecox, L. umbellata, and L. obtusiloba. We defined the pseudo-variant rate as the genetic difference between replicates of the same individual.

To examine the distribution of genetic variation within and among studied locations and to explore genetic relationships among individuals and geographic genetic structure, pairwise individual-by-individual genetic distance matrixes were calculated for codominant data [24] using GenAlEx v.6.5. To deal with missing data, we selected the option "interpolate missing". Based on the matrixes, we conducted analysis of molecular variance (AMOVA) using the R package “pegas” v.0.14 [25] and a principal coordinate analysis (PCoA) to understand the genetic relationships among individuals and among locations in detail using the R package “vegan” v.2.5 [26]. To evaluate the relationship between geographic distance and genetic divergence, we calculated a fixation index (F_ST) value between populations, and assessed the correlation between F_ST and the geographic distance between populations using the R package “ape” v.5.4 [27].

3. Results

3.1. Genetic Diversity Parameters

For L. glauca, we obtained 246 SNP loci from 147 samples, including 105 individuals, replicates of seven individuals and 35 seeds, taken from 29 populations. The mean Na and Ne were 0.903 and 0.894, respectively. He varied from 0.000 to 0.018, whereas Ho varied from 0.000 to 0.027. F_IS was negative excluding Hagi (HG) (range: −0.760 to 0.429; mean = −0.441) (

Table 1. Genetic diversity parameters of the 29 populations of Lindera glauca (Lauraceae) in Japan based on 246 single nucleotide polymorphism (SNP) loci. Letters in parentheses are abbreviations of sampling location.

Sampling Location	N	Na	Ne	He	Ho	FIS
Sendai, Miyagi (SD)	1	0.687	0.687	0.004	0.008	-
Hachioji, Tokyo (HT)	2	0.796	0.792	0.012	0.021	−0.760
Sagamihara, Kanagawa (SH)	2	0.936	0.927	0.012	0.019	−0.538
Kawasaki, Kanagawa (KW)	3	0.909	0.896	0.01	0.014	−0.401
Nagano, Nagano (NG)	5	1.023	1.005	0.014	0.021	−0.504
Kasugai, Aichi (I)	5	1.023	1.009	0.012	0.018	−0.461
Aioiyama, Aichi (A)	5	1.034	1.021	0.014	0.02	−0.381
Hino, Shiga (HN)	3	0.002	0.003	0.01	0.018	−0.749
Takaragaike, Kyoto (T)	5	1.023	1.017	0.013	0.018	−0.328
Kasugayama, Nara (K)	5	1.019	1.008	0.011	0.018	−0.598
Myokensan, Osaka (M)	5	1.038	1.029	0.018	0.027	−0.499
Shigisan, Osaka (S)	5	1.026	1.020	0.011	0.016	−0.434
Takedao, Hyogo (TK)	4	0.996	0.981	0.014	0.02	−0.496
Okayamashi, Okayama (O)	5	1.000	0.995	0.008	0.012	−0.578
Hiroshimashi, Hiroshima (H)	4	0.909	0.900	0.01	0.016	−0.665
Masuda, Shimane (MS)	5	1.038	1.022	0.016	0.019	−0.222
Tsuwano, Shimane (TW)	2	0.977	0.972	0.006	0.009	−0.527
Izumo, Shimane (IZ)	5	1.038	1.023	0.013	0.017	−0.263
Yamaguchishi, Yamaguchi (YG)	5	1.023	1.012	0.007	0.01	−0.360
Hagi, Yamaguchi (HG)	2	0.962	0.961	0.003	0.002	0.429
Bizan, Tokushima (B)	3	1.019	1.015	0.012	0.019	−0.565
Iinoyama, Kagawa (IN)	5	0.985	0.961	0.014	0.019	−0.391
Shiunzan, Kagawa (SU)	5	1.049	1.021	0.018	0.024	−0.306
Tsuno, Kochi (TN)	3	0.887	0.881	0.009	0.015	−0.618
Sameura, Kochi (SM)	5	1.038	1.018	0.011	0.016	−0.377
Tsushima, Nagasaki (TS)	1	0.596	0.596	0.008	0.015	-
Umamiyama, Fukuoka (UM)	3	0.898	0.886	0.011	0.016	−0.435
Sasga, Saga (SA)	1	0.811	0.811	0.008	0.015	-
Sannotake, Kumamoto (SN)	1	0.457	0.457	0	0	-
Regional mean		0.903	0.894	0.011	0.016	−0.441

N, number of individuals; Na, number of observed alleles; Ne, effective number of alleles; He, expected heterozygosity; Ho, observed heterozygosity; F_IS: inbreeding coefficient.

).

For L. praecox, we obtained 355 SNP loci from 114 samples, including 45 individuals, replicates of six individuals and 63 seeds, taken from 11 populations. The mean Na and Ne were 1.025 and 0.968, respectively. He ranged from 0.008 to 0.075, whereas Ho ranged from 0.016 to 0.070. F_IS was positive in all populations excluding Hiroshima (H) (range: −0.223 to 0.151; mean = 0.057) (

Table 2. Genetic diversity parameters of the 11 populations of Lindera praecox (Lauraceae) in Japan based on 355 SNP loci. Letters in parentheses are abbreviations of sampling location.

Sampling Location	N	Na	Ne	He	Ho	FIS
Sendai, Miyagi (SD)	1	0.402	0.402	0.008	0.016	-
Tsukubasan, Ibaraki (TB)	5	1.136	1.064	0.063	0.056	0.108
Hachioji, Tokyo (HT)	1	0.786	0.786	0.022	0.043	-
Nagano, Nagano (NG)	5	1.125	1.059	0.059	0.054	0.094
Hino, Shiga (HN)	5	1.146	1.073	0.070	0.070	0.010
Mitarai, Nara (MT)	5	1.164	1.083	0.061	0.057	0.072
Kongosan, Osaka (KG)	5	1.195	1.114	0.075	0.068	0.094
Takedao, Hyogo (TK)	5	1.160	1.072	0.062	0.052	0.151
Hiroshimashi, Hiroshima (H)	3	0.996	0.960	0.043	0.052	−0.223
Tsuno, Kochi (TN)	5	1.094	1.031	0.055	0.047	0.140
Umamiyama, Fukuoka (UM)	5	1.070	1.009	0.054	0.050	0.071
Regional mean		1.025	0.968	0.052	0.051	0.057

N, number of individuals; Na, number of observed alleles; Ne, effective number of alleles; He, expected heterozygosity; Ho, observed heterozygosity; F_IS: inbreeding coefficient.

).

For L. umbellata, we obtained 280 SNP loci from 38 samples, including 27 individuals, replicates of six individuals and five seeds, taken from six populations. The mean Na and Ne were 1.210 and 1.115, respectively. He ranged from 0.056 to 0.096, whereas Ho ranged from 0.075 to 0.114. F_IS was negative in all populations excluding Hino (HN) and Myokensan (M) (range: −0.356 to 0.020; mean = −0.142) (

Table 3. Genetic diversity parameters of the six populations of Lindera umbellata (Lauraceae) in Japan based on 280 SNP loci. Letters in parentheses are abbreviations of sampling location.

Sampling Location	N	Na	Ne	He	Ho	FIS
Sendai, Miyagi (SD)	5	1.120	1.102	0.071	0.089	−0.252
Tsukubasan, Ibaraki (TB)	5	1.234	1.142	0.096	0.114	−0.183
Hino, Shiga (HN)	5	1.268	1.133	0.090	0.088	0.020
Mitarai, Nara (MT)	5	1.204	1.109	0.080	0.087	−0.089
Myokensan, Osaka (M)	5	1.257	1.138	0.084	0.084	0.006
Takedao, Hyogo (TK)	2	1.102	1.069	0.056	0.075	−0.356
Regional mean		1.210	1.115	0.079	0.089	−0.142

N, number of individuals; Na, number of observed alleles; Ne, effective number of alleles; He, expected heterozygosity; Ho, observed heterozygosity; F_IS: inbreeding coefficient.

).

For L. obtusiloba, 688 SNP loci were obtained from 68 samples, including 39 individuals, replicates of six individuals and 23 seeds, taken from eight populations. The mean Na and Ne were 1.227 and 1.124, respectively. He ranged from 0.071 to 0.097, whereas Ho ranged from 0.073 to 0.084. F_IS was close to zero in all population (range: −0.087 to 0.149; mean = 0.022) (

Table 4. Genetic diversity parameters of the eight populations of Lindera obtusiloba (Lauraceae) in Japan based on 688 SNP loci. Letters in parentheses are abbreviations of sampling location.

Sampling Location	N	Na	Ne	He	Ho	FIS
Nagano, Nagano (NG)	5	1.268	1.14	0.097	0.083	0.149
Hino, Shiga (HN)	5	1.233	1.133	0.085	0.083	0.027
Mitarai, Nara (MT)	5	1.203	1.108	0.078	0.084	−0.078
Myokensan, Osaka (M)	5	1.248	1.134	0.084	0.073	0.133
Kongosan, Osaka (KG)	5	1.264	1.145	0.093	0.084	0.095
Takedao, Hyogo (TK)	5	1.249	1.138	0.085	0.084	0.008
Iinoyama, Kagawa (IN)	4	1.177	1.103	0.072	0.079	−0.087
Tsushima, Nagasaki (TS)	5	1.172	1.095	0.071	0.076	−0.070
Regional mean		1.227	1.124	0.083	0.081	0.022

N, number of individuals; Na, number of observed alleles; Ne, effective number of alleles; He, expected heterozygosity; Ho, observed heterozygosity; F_IS: inbreeding coefficient.

).

3.2. Clone Identification and Genotype Diversity

For L. praecox, L. umbellata, and L. obtusiloba, pseudo-variant rates ranged from 0.005 to 0.013, 0.005−0.012, and 0.004−0.070, respectively (Figure 1a–c). Pseudo-variant rates in the seven pairs of replicated L. glauca samples ranged from 0.004 to 0.017 (Figure 1d).

For L. praecox, L. umbellata, and L. obtusiloba, the genetic difference between a mother plant and its seeds, the genetic difference within populations, and the genetic difference across the entire study area were significantly greater than would be expected from pseudo-variant rate (p < 0.05; Kruskal–Wallis test and Wilcoxon test with sequential Bonferroni correction; Figure 1). Additionally, the genetic differences were significantly different among these three scales (p < 0.05; Wilcoxon test with sequential Bonferroni correction; Figure 1a–c). For L. glauca, the genetic difference at these three scales was not significantly different than what would be expected from pseudo-variant rates (p > 0.05, Kruskal–Wallis test and Wilcoxon test with sequential Bonferroni correction; Figure 1d). However, the genetic difference among the three scales was significantly different from each other (p < 0.05; Wilcoxon test with sequential Bonferroni correction; Figure 1d).

3.3. Genetic Structure

Results of AMOVA indicated that genetic variation within populations was greater than among populations for L. praecox, at 56% and 44%, respectively (

Table 5. Analysis of Molecular Variance (AMOVA) results for the studied four Lindera populations.

Source of Variation		df	Sum of Squares	Variance Components	%
(A)	L. praecox
	Among population	10	427,635.500	8034.000	44
	Within population	34	352,842.900	10,378.000	56
	Total	44	780,478.400	18,412.000
(B)	L. umbellata
	Among population	5	57,977.820	2084.700	47
	Within population	21	48,930.350	2330.000	53
	Total	26	106,908.170	4414.700
(C)	L. obtusiloba
	Among population	7	477,950.100	9319.100	76
	Within population	31	709,215.400	2877.900	24
	Total	38	1,187,165.500	12,197.000
(D)	L. glauca
	Among population	28	1005.172	2.706	9
	Within population	76	2007.043	26.066	91
	Total	104	3012.215	28.772

). For L. umbellata, genetic variation was greater within populations (53%), followed by among population (47%). Genetic variation among populations of L. obtusiloba (76%) was substantially larger than within populations (24%). For L. glauca, geographic differences contributed to only 9% of the genetic variation within this species.

Pairwise Mantel tests for L. praecox (Figure 2a), F_ST and distance were positively correlated (p = 0.025). For L. umbellata (Figure 2b) indicated no statistically significant correlation between F_ST and the geographic distance between populations (p = 0.897). Lindera obtusiloba and L. glauca (Figure 2c,d) also showed F_ST and distance were positively correlated (p = 0.023 and 0.010, respectively).

The PCoA of L. praecox and L. obtusiloba appears to have genetically distinct variation within their range from east to west (Figure 3a,c), and geographic genetic differentiation of L. umbellata was clear but its pattern was complex (Figure 3b). The PCoA plot of L. glauca showed no clear differentiation in genetic composition throughout its geographic range in Japan (Figure 3d).

4. Discussion

A key advantage of MIG-seq is that it can provide putative neutral loci [28], which are useful for analyses of genetic structure. Therefore, this approach has been used recently in population genetics (e.g., [29,30]).

Obligate outcrossing species is expected to have the values of Ho close to those of He as reported by Ye et al. [31] for L. obtusiloba. For the genetic diversity of the dioecious tree species Phoebe zhennan (Lauraceae), Ho was higher than He (Ho = 0.243, He = 0.145), indicating low inbreeding [32]. Our results are similar to those of previous studies. This may be an advantage of dioecy, which avoids inbreeding. Our study results also suggest that dioecious plants have greater heterozygosity than cosexual plants, based on a comparison of our F_IS values with previously published estimates based on the genome-wide SNPs data (Table 2, Table 3 and Table 4), 0.035–0.091 (Rhododendron japonoheptamerum) [29]; 0.05–0.19 (Platycladus orientalis) [33]; 0.05 (Thuja koraiensis) [33]; 0.051–0.754 (Tectona grandis) [30]). However, the level of heterozygosity is influenced not only by reproductive system of the focal species, but also other factors, e.g., population size [34], age of individuals [35]. Thus, careful verification is needed.

The all-female populations of L. glauca assessed here showed higher heterozygosity than the other three Lindera species, on average, and exhibited negative F_IS values at all sites except for one site. These results are similar to a previous study that employed nuclear microsatellite analyses [18]. Zhu et al. [18] suggested that seven Japanese populations of L. glauca exhibited greater heterozygosity than what was observed across a large range of populations in China. For an apomictic lineage, the high and well-preserved heterozygosity in L. glauca populations in Japan may attribute to asexual (apomictic) reproduction as reported for Prunus avium [36], due to the absence of sexual recombination and/or the lack of accumulation of somatic mutations over time.

We found evidence of genetic differentiation among populations of L. praecox, L. umbellata, and L. obtusiloba as has often been reported for other plant species [37,38,39]. First, AMOVA results indicated that regional variation explained a large proportion of the total variation (Table 5). Second, the PCoA results indicated geographic divisions in populations. Finally, we found that populations that were further apart in space were more genetically different in L. praecox and L. obtusiloba (Figure 2a,c). Collectively, these results reflect geographically limited gene flow, a varied history of population expansion, or strong adaptation to local environmental conditions.

By contrast, populations of L. glauca had substantially lower genetic variation than these three species. Even in apomictic plants, genetic diversity can occur in the process of seed production by incomplete meiosis and/or recombination [40]. In the case of L. glauca in Japan, the comparison of multi-level genetic distance and PCR error rates indicated no significant differences for any level, which suggests that seeds were exact copies of their mothers, as the genetic difference between mothers and seeds were within the range of PCR errors (Figure 1d). A seed production process that does not produce mutations would lead genetically uniform populations in Japan.

Here, we compare our own results in Japan with available results from China. Lindera obtusiloba produces seeds by sexual mating both in China and Japan. Genetically distinct populations in different regions were observed in China [31]. Our results are similar to those of the previous study. On the other hand, the situation is different for L. glauca. Both females and males of Lindera glauca exist in China, and they reproduce sexually. A variety of genotypes existed in different regions and geographic genetic structure was observed in China [18]. While in Japan, seed production by apomixis resulted in an extremely uniform genetic structure. These results of comparisons suggest that reproductive system has a strong influence on genetic structure of plant populations.

We were surprised that virtually no genetic variation was observed in L. glauca across its entire geographic range in Japan. The genetic difference between individuals in the whole range of Japan was not significantly different from PCR errors (Figure 1d). In addition, we found significant but weak correlation between genetic and geographic distance (Figure 2d), and the PCoA results indicated unclear geographic distinctions among populations (Figure 3d). These are contrast to the previous studies which detected a clear geographic genetic structure of apomictic plants, Miconia albicans [7] and Taraxacum officinale [9]. These results suggested that a population of L. glauca in Japan might be a huge clone. Bricker et al. [41] reported that a dioecious aquatic plant, Thalassia testudinum, formed a “mega-clone” that spread over nearly 50 km by vegetative fragments. The European populations of Gagea spathacea (Liliaceae) were dominated by a single “mega-clone” spreading all over Central Europe [42]. Similarly, the L. glauca population in Japan could also be considered as a “mega-clone” that spread out over at least 1100 km in distance from Sendai to Kumamoto. There are two possible reasons for the formation of genetically uniform populations. First, Japanese populations of L. glauca might be founded by one asexually reproducing female as suggested by Dupont [17]. Second, population of L. glauca in Japan may still be young on the geological time scale. The population of L. glauca has been in Japan since prehistoric times because Pleistocene fossils of L. glauca have been found [11]. Genetic diversification via accumulated mutations has not yet occurred, although it may occur in the future. However, there were significant dissimilarities in the degree of the genetic difference between a mother plant and its seeds, between individuals in a population, and between individuals in the entire study area (Figure 1d). The result suggested that slight mutations may have occurred and accumulated.

Although L. glauca has been present in Japan since the Pleistocene (0.01–2.58 million years ago (Ma)) [11], the genetic diversity of this species has remained at a fairly low level until now. According to Xiong et al. [19], wild L. glauca in mainland China have a highly skewed sex ratio with females predominating, and their genetic diversity is thought to be maintained by sexual reproduction with a small number of male plants. It is not clear whether Japanese L. glauca can cross with the Chinese males. However, if they can, immigration of males may increase genetic diversity as well as the Chinese L. glauca.

Finally, the loss of genetic diversity relates to an increase in extinction risk [43]. Despite this, populations with extremely low genetic diversity have been reported worldwide [42,44]. In the apomictic L. glauca, the advantages of population establishment or expansion by a single individual, without mating, may outweigh the disadvantages associated with the extremely low genetic diversity. If environmental conditions remain suitable for the extant genotype, genetic diversity may not be required to maintain these populations, at least in the short term.

5. Conclusions

The three species that produce seeds by sexual reproduction, Lindera praecox, L. umbellata, and L. obtusiloba had lower values for inbreeding coefficient, indicating that there is indeed an effect of dioecy. Genetic variations were also found both at the local scale and the geographic scale. It is possible that they have higher genetic diversity than cosexual plants, and the effect of inbreeding depression may be suppressed. On the other hand, in L. glauca, the advantage of being able to found populations from single individual, without mating, may outweigh the disadvantages associated with extremely low genetic diversity. This may explain why this species is so widely distributed in Japan. Provided that the current genotypes remain suited to environmental conditions, L. glauca may not be constrained by its limited genetic diversity.