Population genetic structure and variability in Lindera glauca (Lauraceae) indicates low levels of genetic diversity and skewed sex ratios in natural populations in mainland China

Sample collection

During field expeditions carried out from 2013 to 2017, 300 individuals were sampled from 20 wild populations and two cultivated populations, representing nearly the entire natural distribution of L. glauca in mainland China (Table 1; Fig. 1A). Most L. glauca individuals are able to form clones via vegetative reproduction with stolons (Tsui, Xia & Li, 1982), as found in our field survey. In order to avoid the collection of several ramets from the same genet, a single sample was obtained from each cluster of shrubs in close proximity to a main tree, excluding the surrounding young branches growing on the ground. Each sample (individual) in same population was collected at least 10 m apart. In some smaller populations, fewer than 10 plants of putatively nonclonal origin were available. Overall sample sizes varied from 5 to 30 per wild population (Table 1). There was definite apomixis in the individuals of two cultivated populations (based on our survey for three years), and we added it to global analysis as a reference for wild populations. In the field, fresh leaves were immediately dried in silica gel after collection, and preserved at room temperature until DNA extraction.

DNA extraction and microsatellite genotyping

Genomic DNA was extracted from 100–150 mg of dried leaves per sample using a modified cetyltrimethylammonium bromide (CTAB) method (Doyle & Doyle, 1987). Microsatellite loci of all 300 individuals of L. glauca were screened, including 13 polymorphic nuclear microsatellite markers (EST-based microsatellites; hereafter nSSRs; Table S1) and five polymorphic chloroplast microsatellite markers (hereafter cpSSRs; Table S2). All nSSRs and cpSSRs were labeled with fluorescently labeled nucleotides (forward primer with M13F) and detected by capillary gel electrophoresis. Subsequent\ steps and the PCR assay were conducted according to Xiong et al. (2016). Genotyping was performed using an ABI 3730XL DNA Analyzer (Applied Biosystems, California, USA ) with a GeneScan 500 LIZ Size Standard, and alleles for each locus were manually scored using GeneMarker version 2.2.0 software (SoftGenetics, State College, PA, USA).

Data analyses

Raw data matrices containing information of alleles and haplotypes for 13 nSSR and 5 cpSSR loci were checked for scoring errors. All SSR analyses were conducted with 300 samples. Data editing and formatting were performed using GenAlEx v. 6.502 (Smouse, Whitehead & Peakall, 2015).

The related indexes of genetic diversity were calculated as followed: for the nSSR data set, genetic diversity indices, including the number of alleles (N_A), observed heterozygosity (Ho), expected heterozygosity (He), percentage of polymorphic loci (PPB), Wright’s inbreeding coefficient (F_IS), and Nei’s (Nei, 1978) genetic distances, were estimated using GenAlEx v. 6.502 (Smouse, Whitehead & Peakall, 2015) and POPGENE v. 1.32 (Yeh, Yang & Boyle, 1999). The online package GENEPOP v. 4.1.4 (Rousset, 2008) was used to perform exact Hardy–Weinberg equilibrium (HWE) tests and to test for the presence of private (null) alleles. The differentiation index F_ST was computed for pairs of populations using Arlequin v. 3.5.1.3 (Excoffier, Laval & Schneider, 2005). Allelic richness (A_R) was calculated using the software FSTAT v. 1.2 (Goudet, 1995). For the cpSSR data set, number of haplotypes (Nb), genetic diversity (Dv), haplotype richness (Hrs), the number of private alleles (Prv), and the polymorphism information content (PIC) per locus were estimated using HAPLOTYPE v. 1.05 (Eliades & Eliades, 2009).

The population genetic structure was analyzed as followed: for the nSSR data set, the genetic structure of the 22 populations (20 wild and 2 cultivated) was analyzed using the Bayesian clustering approach implemented in STRUCTURE v. 2.3.4 (Pritchard, Stephens & Donnelly, 2000), assuming an admixture model. In order to determine the most appropriate number of genetic clusters or groups (K value), K was set from 1 to 20, and the analysis was run with 20 iterations for each K with a burn-in of 1,000,000 generations followed by 50,000 generations for the Markov chain Monte Carlo (MCMC) simulation. The admixture level for each individual (Q) was also inferred. The program STRUCTURE HARVESTER v. 0.6.94 (Earl & Von Holdt, 2012) was used to estimate the number of population clusters based on the ΔK parameter according to Evanno, Regnaut & Goudet (2005). Based on the most appropriate number of clusters suggested by Bayesian clustering, analysis of molecular variance (AMOVA) was performed using Arlequin, with 10,000 iterations for the permutation test. A neighbor-joining (NJ) tree was generated using POPTREE2 (Takezaki, Nei & Tamura, 2010) based on pairwise Nei (1978) genetic distances between populations determined by GenAlEx. For the cpSSR data set, the Arlequin was used to determine pairwise F_ST values among all populations. A parsimony network illustrating genetic relationships among haplotypes of L. glauca populations was generated using TCS v.1.1 (Clement, Posada & Crandall, 2000).

Considering that isolation by distance (IBD) can be a key factor keeping populations apart by limiting gene flow (Coyne & Orr, 2004), the IBD of wild L. glauca inter-population in mainland China was tested. In view of the potential importance of pollen in the spread of apomixis (Mogie, 1992; Whitton et al., 2008), the pollen/seed migration ratio (r) was calculated. In order to examine IBD, the Mantel test was performed using GenAlEx, correlating the pairwise genetic distances [F_ST∕(1 − F_ST)] with the pairwise geographic distances (in kilometers). To calculate r, we used the followed formula: r = mp∕ms = [(1∕F_ST(n) − 1)(1 + F_IS) − 2(1∕F_ST(cp) − 1)]∕(1∕F_ST(cp) − 1) (Ennos, 1994; Petit et al., 2005), where mp is the pollen migration rate, ms is the seed migration rate, F_ST values (overall F_ST) are population differentiation estimates derived from AMOVA, F_ST(n) is the nuclear F_ST and F_ST (cp) is the chloroplast F_ST.

Population bottlenecks were evaluated using BOTTLENECK v. 1.2.02 (Piry, Luikart & Cornuet, 1990) with the infinite alleles model (IAM) that a single mutation is allocated at a time and the resulting number of alleles is computed, stepwise mutation model (SMM) that is a Bayesian approach and generally more appropriate when testing microsatellite loci, and two-phased model (TPM) that is a modified SMM. According to Piry et al., sign tests, Wilcoxon tests, and mode-shift were applied, excluding standardized differences tests, which are useful when at least 20 polymorphic loci are available.

Genetic variation

A total of 74 alleles and 13 haplotypes were identified at 13 nSSRs and 5 cpSSRs across 300 individuals of L. glauca. For each locus, the number of alleles for 13 nSSR loci ranged from 3 (P-298) to 8 (XBLG-060), with a mean of 5.7 alleles (Table S1). In particular, A_R ranged from 1.807 to 2.774, with a mean of 2.329, and PIC ranged from 0.363 to 0.711, with a mean of 0.556. Ho and He varied between 0.210 and 0.563, with a mean of 0.380, and between 0.380 and 0.754, with a mean of 0.602.

The total number of alleles of 20 wild populations across 13 nSSR loci varied from 22 (population SZY) to 44 (population ZJJ), with a mean of 32.9, and allelic richness ranged from 1.231 to 2.011, with a mean of 1.740 (Table 2). Population Ho ranged from 0.108 to 0.708 and population He from 0.106 to 0.477, with means of 0.328 and 0.376. The PPB ranged from 53.85% to 100%, with a mean of 85% (Table 2). The values of total number of alleles, Ho, He, and PPB were all smaller in the two cultivation populations, with a mean of 18.5, 0.212, 0.145, and 42.31%, respectively. These parameter values of Ho (0.212), He (0.145), and PPB (42.31%) in two cultivated populations were much lower. Significant deviations from the HWE indicating a heterozygote deficiency were detected in 9 of 22 populations. For some populations (KYS, WYS, ZJS, WJS, GJS, NHS, FHS, and ZJJ), negative F_IS values within the populations were observed, indicating more heterozygotes than expected. However, none of the 13 nSSR loci with the heterozygote excesses appear when calculating it on all samples (Table S1), and there is no evidence of private alleles within the data set. Among all wild samples (290), 277 individual plants had a unique multi-locus pattern after PCR amplification with 13 nSSRs primers, indicating that these samples were from different individuals. Of the remaining 13 individual plants, 5 pairs exhibited the same multi-locus pattern in pairs, and 3 individuals exhibited the same multi-locus pattern.

Genetic diversity parameters for cpSSR loci are summarized in Table S2. All 5 cpSSR loci exhibited 2–3 alleles per locus across all samples. Dv ranged from 0.115 to 0.249 per locus, and PIC varied from 0.155 to 0.218. Analyzing combinations of all alleles, there were 22 unique haplotypes (hereafter H; Fig. 1). All populations contained several haplotypes, except for populations ATM, YTH, NTB, FHS, TMS, and HZY (Fig. 1A). The network of plastid haplotypes was complex (Fig. 1B). Haplotype H12 exhibited the highest frequency and was detected in 18 of 20 wild populations (Table S3). Of the 22 haplotypes, 7 were private haplotypes excluding cultivated populations, Hrs per population ranged from 0 (populations ATM, NTB, YTH, TMS, and FHS) to 0.964 (population FJS), and the mean genetic distance between individuals (${\mathrm{D}}_{\mathrm{sh}}^2$) varied from 0 to 39.638 (population JGS) (Table 2).

Genetic clustering and population differentiation

A Bayesian analysis based on 13 nSSRs implemented in STRUCTURE showed the presence of 2 clusters (K = 2), with only slight admixture at the individual level in each population, except for population ATM (Fig. 2A). The ΔK statistic developed by Evanno et al. indicated that the overall differences were not substantial (Fig. 2B). Cluster orange included 13 wild populations (ATM, JGS, LDZ, SJG, NTB, YTH, DBS, HMF, TMS, SQS, LYS, KYS, and WYS), and the remaining 7 wild populations and 2 cultivated populations were assigned to cluster blue. The NJ tree (Fig. 3A) and principal coordinates analysis (Fig. 3B) based on the nSSR dataset supported the results of STRUCTURE analysis, indicating that 22 populations could be grouped into 2 clusters. However, the network diagram of all 22 unique plastid haplotypes revealed by 5 cpSSRs was complex (Fig. 1B), and failed to support the 2 distinct clusters revealed using nuclear data.

The 2 clusters revealed by STRUCTURE analysis were set as groups for AMOVA based on both the nSSR and cpSSR dataset (Table 3). Using the nSSR dataset, the majority of genetic variation was detected within populations (63.82%), indicating a genetic differentiation mostly at the individual level. Nevertheless, a considerable proportion of the total variation (26.42%) was found among populations within groups, and a small amount of variation (9.76%) occurred among the 2 groups. In contrast to the nSSR results, the AMOVA based on the cpSSR dataset showed that a larger proportion of genetic variation could be attributed to variation within populations (68.84%) and among populations within groups (29.37%), and little variation among groups (1.79%). The overall F_ST values calculated by AMOVA were 0.362 (P ≤ 0.0001) for the nSSR dataset and 0.312 (P ≤ 0.0001) for the cpSSR dataset.

Isolation by distance and pollen/seed migration ratios

The estimates of genetic differentiation (F_ST value) based on 13 nSSRs ranged from 0.023 (between WJS and GJS) to 0.427 (between LDZ and FHS) (Table S4), excluding 2 cultivated populations. Only 4 pairwise comparisons (ATM and JGS, WJS and NHS, NTB and YTH, and WJS and GJS) showed significant Fst values (P ≤ 0.05). Adopting a P-value of 0.01, no significant correlation between pairwise genetic distance [F_ST/(1 − F_ST)] and geographic distance (in kilometers) was found using the Mantel test (Fig. 4) for the nSSR dataset (R² = 0.0401, P = 0.068) or the cpSSR dataset (R² = 0.033, P = 0.091), suggesting that L. glauca in mainland China does not exhibit significant IBD. The F_ST values for the nSSR and cpSSR were similar (Table 3), and the r was −1.149, indicating that most gene flow among populations occurs via seed, rather than via pollen.

Population bottleneck effect

Several populations had a significant excess of heterozygosity expected at mutation-drift equilibrium (i.e.,  He > Heq) (Piry, Luikart & Cornuet, 1990) under the 3 models in the bottleneck analysis, which indicated a deviation from mutation drift equilibrium in wild L. glauca populations (Table 4). More specifically, population GJS exhibited a significant bottleneck event according to the sign and Wilcoxon tests in all 3 models, indicating a population size decline (bottleneck effect) in its history. Population WJS experienced a significant bottleneck event according to the sign test and Wilcoxon test for the IAM and TPM methods, and FHS exhibited a significant bottleneck event by Wilcoxon test in all 3 models. SJG, FJS, and ZJS only exhibited a significant bottleneck event based on the sign and Wilcoxon tests for the IAM method (Table 4).

Genetic variation within populations

In this study, sampling covered a large portion of the natural distribution, and overall genetic diversity across wild L. glauca populations exhibited low levels based on both nSSR (mean A_R = 1.74, Ho = 0.33, He = 0.38, F_IS =  − 0.14) and cpSSR (mean Nb = 2.06, Hrs = 0.39) loci. Our estimates of genetic diversity in L. glauca were almost half those of long-lived perennials (Ho = 0.63, He = 0.68), out-crossing species (Ho = 0.63, He = 0.65), and plants with wide distributions (Ho = 0.57, He = 0.62) (Nybom, 2004), which were lower than Laurus nobilis in Lauraceae (A_R = 3.22, He = 0.56) (Marzouki et al., 2009). There are several major factors influencing variation that can result, each by itself or in combination, in the low levels of genetic variation observed in wild L. glauca populations. Asexual reproduction through apomictic seeds can decrease genetic variation in a population especially in apomictic populations (Lo, Stefanovic & Dickinson, 2009). Similarly, effective population size could seem like a limitation, due to the population being established from a limited number of individuals and a small sampling quantity. For two cultivated populations, these values were much lower than in wild populations, perhaps because apomixis reduces genetic diversity, or because of the potential confounding effect of small population sample size. Furthermore, the estimated values of genetic diversity therein are also lower than nSSR-based values found in literature (A_R = 2.61, He = 0.44, F_IS =  − 0.37) (Zhu et al., 2016). Differences in genetic variation between this study and the previous one are likely to be explained by the number of sampling populations and individuals the study of Zhu et al. (2016) included 6 populations, a total of 96 individual plants, while this study included 20 wild populations, making a total of 290 individual plants.

On the other hand, a marked similarity in the molecular variance revealed by the 2 types of markers (overall Fsc = 0.293 and F_ST = 0.362 for nSSRs; 0.299 and 0.312 for cpSSRs) was observed, indicating general consistency between chloroplasts and nuclear DNA. More specifically, predominant apomixis in wild L. glauca could explain how F_ST observed for nuclear markers is a little higher than that for chloroplast markers, and vice versa. In addition, apomictic reproduction of L. glauca could also affect the results. Therefore the use of clone-corrected data (removing the data of clones from the same parent) is necessary if clones are detected, because L. glauca trees could form clones via vegetative reproduction with stolons.

Nine populations had more heterozygotes than expected. STRUCTURE analysis showed that the 9 populations were grouped into one group and the rest into another, suggesting negative F_IS is an important factor affecting group of population difference. Usually, positive F_IS values within populations (JGS, LYS, and FJS) indicated inbreeding, and negative F_IS of wild populations (KYS, WYS, ZJS, WJS, GJS, NHS, FHS, and ZJJ) suggested outbreeding. However, considering that clonality probably generates significant negative F_IS in some wild plant populations with asexual reproduction when considering all individuals (Stoeckel et al., 2006), the observed negative F_IS of wild populations (KYS, WYS, ZJS, WJS, GJS, NHS, FHS, and ZJJ), coupled with the result grouped by NJ tree and the cluster pattern of the NJ tree (blue cluster) (Fig. 3A), suggest that apomixis may have been common in these populations. On the other hand, global multilocus tests indicated that sexual reproduction did exist in these populations that were positive F_IS values (0 < F_IS < 1). There is mixed reproductive system existed in this species, thus explaining the F_IS pattern.

Genetic differentiation among populations

Populations can cluster according to habitat type or geographic distance. However, for species with predominantly asexual populations, like Daktulosphaira vitifoliae (Vorwerk & Forneck, 2006), Crataegus douglasii (Lo, Stefanovic & Dickinson, 2009), and Taraxacum officinale (Majesky et al., 2012), no significant correlation exists between genetic distances and geographic distances. In this study, a correlation between genetic distances (as measured by [F_ST∕(1 − F_ST)] values) and geographic distances (in kilometers) was not detected in L. glauca populations, regardless of marker type, suggesting that weak barriers to dispersal between populations and/or the common founders between neighbors, distant populations, and apomictic populations did not completely limit gene flow. Results indicate that sexual dispersion and apomixis co-occurred in the same natural population. However, further research is needed to investigate the extent to which apomixis limits gene flow, and the exact rate of sexual production and apomixis that occurred among and within populations.

According to Ennos’ formula (1994), the migration ratio of pollen flow versus seed flow (r) in this study was negative (−1.149), suggesting that gene flow between populations is more likely to involve seeds. Usually, the value of r is positive in almost all angiosperm species (Ennos, 1994). The biased r value could be explained by the reproductive mode in different populations being the apomictic seeds of females, coupled with the result that no staminate flower were found at field observation sites for 5 consecutive years. Considering that a situation where a few individuals can self-pollinate is likely, pollen from even just a few staminodes of female flowers may be associated with the negative r value of pollen flow versus seed flow. However, because the formula is derived for a hermaphrodite species, it needs to be modified to account for the disproportionate maternal contribution from the females to the next generation if the exact r value of L. glauca can be expected to be observed.

Source of evolutionary potential of apomicts

Although 22 haplotypes were observed across 20 wild populations and 2 cultivated populations of L. glauca, H12 accounted for 62.03% of overall haplotypes and was detected in all populations (Table S3), except for FHS, LYS, SZY, and HZY. This haplotype existing in many populations separated by considerable geographical distances (e.g., greater than 1,720 km between KYS and GJS), coupled with the complex network (Fig. 1B), suggested three inferences to account for the observed result. First, these individuals from different populations had a relatively recent shared female founder. Second, there was an apomictic lineage for this species within some populations. Third, the genome of this species has genetic traits that lead to apomixis, which might be induced by some factors (e.g., biological stimulation, environmental influence, climate changing, etc.), and coexists with sexual reproduction in identical individual plants.

The first inference could be explained by a hypothesis that the H12 haplotype may be associated with the migratory patterns of some birds responsible for the dispersion of apomictic seeds over a long distance. However, given that H12 is the most frequent and connected haplotype, probably is ancestral, this inference is not reliable. Besides, despite conducting field observations for 5 consecutive years, few birds were found eating grown fruits of L. glauca, as well as few small mammals (e.g., Paguma larvata). The second inference means having an early maternal ancestor through apomictic reproduction for many individuals. However, all the living species are at the tips of the tree of life, and apomixis is a derived condition (APG, 2003; APG, 2016; Horandl, 2006; Thompson & Ritland, 2006; Lo, Stefanovic & Dickinson, 2009). Consequently, having ancestral asexual angiosperms is almost impossible, because asexuals fail to maintain sex and recombination in populations that are limited in size, therefore are unable to bring together high-fitness alleles that reside in different individuals. The third inference means that apomicts of L. glauca may be of very recent origin and have the ability to apomictically reproduce through mutations or losses of some sexual genes. A similar situation exists in some species, such as Taraxacum officinale that exhibit alternations of asexual and sexual histories of apomicts (Majesky et al., 2012), some hawthorns (Crataegus; Rosaceae) that have a population genetic structure of diploid sexual and polyploid apomicts (Lo, Stefanovic & Dickinson, 2009), and a marbled crayfish (Procambarus virginalis) that reproduces through parthenogenesis (Gerhard et al., 2003; Ewen, 2018). Therefore, when coupled with the results of filed surveys that show apomixis occurred in all sampled populations, this haplotype (H12) is probably a predominant genotype existing in populations that have the ability to apomixis.

Above all, the authors reject the former two inferences that most populations across such long geographical distances have a relatively recent founder or an ancient apomictic ancestor, and believe the last inference that apomixis caused by mutations or losses of related genes makes this maternal haplotype (H12) present in many individuals from different populations more likely.

Conjecture about histories of natural populations

According to the relationship Ne = 4(NmNf)/(Nm + Nf), where Nm represents the number of males, and Nf represents the number of females (Beerli & Palczewski, 2010), the accurate values of Ne could not be calculated because there were no males in the samples collected. Even so, based on the rarity of male individuals in mainland China, dioecious reproduction reported in the past several decades (Tsui, Xia & Li, 1982; Wang, 1972), some specimens of branches of male individuals stored in the China National Herbarium (PE), and many males grown on Taiwan (Zhang, 2007), it is inferred that some natural populations of L. glauca recently experienced a severe bottleneck and male individuals experienced a decline, most likely resulting from anthropogenic causes. However, in order to test the above hypotheses and obtain more accurate results, further field investigations are necessary, including samples from male and female individuals in a population (especially males found on Taiwan), correcting for clone and apomictic reproduction, across the full range of L. glauca, including Japan, South Korea, and Taiwan.

Implications for conservation

Genetic diversity is recognized as an important population attribute for both conservation and evolutionary purposes (Cena et al., 2006). The purpose of the conservation of endangered and threatened species is to maintain their contributions to overall genetic diversity. People usually focus only on endangered species and provide protection, whereas some species that are reduced in genetic diversity also require protection. This study detected a lower level of genetic diversity in L. glauca than that of some other species in Lauraceae. The conclusion is that wild L. glauca populations have female-skewed sex ratios, which was consistent with our field survey and sampling for 5 consecutive years. The destruction of habitats in order to plant other commercial or medicinal crops and felling by local farmers may explain the low frequency of male L. glauca in mainland China, although the species is common and widely distributed in other regions. The authors propose finding male individuals and promoting sexual reproduction in order to maintain this species’ overall genetic diversity.