Population structure and diversification of Gymnospermium kiangnanense, a plant species with extremely small populations endemic to eastern China

Population sampling

The entire distribution area of Zhejiang and Anhui Provinces was surveyed from April 2016 to March 2021. Like most of extremely limited distributed species, G. kiangnanense usually have a population distribution range less than 10,000 m². In our field survey, we found that it is reproduced by seedling, without clones. We collected as many individuals as possible for each population. The collected individuals grow apart from each other for at least 10 m to avoid sampling of clones from the same maternal plant. A population generally contains 10–22 individuals. The sampling locations were recorded using a GPS locator (Table 1; Fig. 1). Leaves from each individual were sampled and placed into a tea bag with sampling code. Then, all the bags with leaves were dried with silica gel, and stored at −20 °C.

Microsatellite marker development

Primers of polymorphic microsatellite loci were designed using the CandiSSR program based on transcriptome data of G. kiangnanense (GenBank accession numbers: SRR27230252, SRR27235853). The parameters were set as follows: flanking sequence length 100, comparison value 1 × 10⁻¹⁰, comparison recognition 95, comparison coverage 95. Finally, SSR primers with ideal amplification efficiency and high polymorphism were selected for genetic analysis (Table S1).

DNA extraction, PCR amplification and detection

Total genomic DNA was extracted using the Plant DNAzol (Doyle & Doyle, 1987) from approximately 0.1 g of silica-dried young leaf samples. The quantity of the obtained DNA was examined using a NanoDrop spectrophotometer, and PCR amplification products were detected using agarose gel electrophoresis. PCR amplifications were conducted in 15.0 ml reactions containing 30.0–100.0 ng genomic DNA, 2.0 mM of each primer, and 7.5 ml 2 × Taq PCR Master mix. The PCR program was as follows: Pre-denaturation at 95.0 °C for 3.0 min, denaturation at 94.0 °C for 2.0 min, annealing at 55.6–59.6 °C for 30.0 s, and extension at 70.0–75.0 °C for 10.0 min with a total of 30 cycles (Table S1). After completion of the amplification, the PCR products were run on a 1.5% agarose gel with the addition of 1.5 μL of EB staining. The gel was placed in 1 × TBE buffer and subjected to electrophoresis at 80 V for 20 min. All the electrophoresis images were shown in Fig. S1. DNA fragments were visualized under UV light and photographed using a Gel Documentation System (Bio-Rad, Hercules, CA, USA). The amplified PCR products were sent to Tsingke Biotech Company (Chengdu, Sichuan, China) for sequencing. The raw reads were shown in Table S1.

Data analyses

Population divergence history

The population history of G. kiangnanense was investigated using the Approximate Bayesian Computation (ABC) approach implemented in DIYABC v.2.1 (Cornuet et al., 2014). Based on the results from STRUCTURE and systematic clustering analysis, the studied populations were divided into two major groups: populations from Zhejiang Province (ZJ) and Anhui Province (AH). The ZJ clade was further divided into two sub-clades: populations from Zhuji City (ZJZ) and Hangzhou City (ZJH). Three possible historical divergence scenarios were proposed for these three groups (Fig. 2).

In scenarios 1, populations of the AH group and ZJ group initially diverged from the common ancestor of the two lineages, followed by recent divergence of the ZJZ group from the ZJH group. In scenarios 2, populations of the AH group and ZJ group initially diverged from the two lineages, followed by recent divergence of the ZJH group from ZJZ group. In scenarios 3, populations from the ZJH, ZJZ and AH groups diverged simultaneously (Fig. 2).

To evaluate the optimum parameter values of the DIYABC model, the entire dataset is subjected to initial calculation under the assumption of uniform distribution. The minimum population size is set to 10, and the maximum is set to 100. The minimum population separation time is set to 10 generations, and the maximum is set to 10,000. The average mutation rate per generation for each locus was set to 5 × 10^–4 with a minimum value of 10^–5 and maximum value of 10^–2. A total of one million pseudo observation datasets (PODs) were simulated, then the operating parameters were adjusted and determined based on the 95% confidence interval of the simulated results until the optimal support plan is obtained. The posterior probability and corresponding 95% CI value for each scenario were calculated using a local linear regression program. The optimal calculating method was selected for principal component analysis (PCA) to check the fit degree between simulated values and observed values.

Primer polymorphism analyses

After screening for loci polymorphism, a total of 21 pairs of SSR primers were used in the study. Agarose gel electrophoresis images showed that the amplification products of the 21 pairs of SSR primers had clear bands (Fig. S1). The number of alleles per locus ranged from two to six with a mean value of 2.900. The expected heterozygosity (H_E) ranged from 0.007 to 0.762 with a mean value of 0.267. The observed heterozygosity (H_O) ranged from 0.007 to 0.979 with a mean value of 0.212 (Table S3). The Ho values were generally smaller than He values, indicating the presence of different sets of heterozygous deletions. The polymorphic information content (PIC) ranged from 0.007 to 0.721, with a mean value of 0.236. Among these, 12 loci showed low polymorphism (≤0.25), eight loci showed moderate polymorphism (0.250 < PIC ≤ 0.500), and one loci showed highly polymorphism (PIC > 0.500). The total populations inbreeding coefficient (F_it) ranged from −0.736 to 0.866, with a mean value of 0.209. The genetic differentiation coefficient (F_st) among populations of G. kiangnanense ranged from 0.021 to 0.780, with a mean value of 0.261. The intra-population inbreeding coefficient (F_is) ranged from −0.949 to 0.839, with a mean value of −0.050. Among these, nine loci had heterozygous deletion, and the remaining twelve loci had some heterozygous excess. The interpopulation gene flow (N_m) ranged from 0.071 to 11.491, with a mean value of 2.559.

Genetic diversity of G. kiangnanense

The genetic diversity of eight G. kiangnanense populations covering its natural distribution areas was analyzed, revealing different degrees of genetic diversity (Table 2). The number of alleles (N_A) and effective alleles (N_E) ranged from 1.619 (ZJH1, ZJH2, ZJH4) to 1.952 (ZJZ1) and from 1.247 (ZJH3) to 1.445 (AH1), respectively. The mean value of N_A was 1.714, and the mean value of N_E was 1.330. The H_O ranged from 0.172 (ZJZ1) to 0.300 (AH1), with a mean value of 0.222. The H_E ranged from 0.143 (ZJH3) to 0.240 (AH1), with a mean value of 0.179. The H_O of ZJZ1 was less than H_E, indicating a low genetic diversity in the population. The Shannon information index (I) ranged from 0.228 (ZJH3) to 0.366 (AH1), with a mean value of 0.286. The mean value of the fixed coefficient (F) for each population was −0.104 (F < 0). The mean value of the inbreeding coefficient (F_IS) for each population was −0.206 (F_IS < 0), indicating an unbalanced distribution of genotypes within the population, the presence of heterozygotes in each population, the presence of inbreeding in ZJZ1, and the existence of inbreeding and distant crosses within the remaining seven populations. The mean percentage of polymorphic loci (PPL) was 57.14%.

Genetic structure of G. kiangnanense

In order to further investigate the genetic clustering relationship and genetic structure among populations of G. kiangnanense, principal coordinate analysis (PCoA) was conducted using GenAlEx. The results showed that the eight natural populations were divided into two groups (Fig. 3). UPGMA cluster analysis based on genetic distance using Poptree also supported the above results (Fig. 4).

STRUCTURE analysis yielded the highest likelihood when samples were clustered into two groups (K = 2; Figs. 5A and 5B). The majority of individuals (96.20%) from Zhejiang Province populations (ZJZ1, ZJZ2, ZJH1, ZJH2, ZJH3, and ZJH4) were assigned to cluster I (‘light blue pool’, Figs. 5 and 6), while individuals (97.90%) from Anhui Province populations (AH1 and AH2) were assigned to cluster II (‘orange pool’, Figs. 5 and 6). In addition, STRUCTURE analysis of the six G. kiangnanense populations in Zhejiang Province showed that K = 2 is the best grouping, with ZJZ1 and ZJZ2 from Zhuji City forming one lineage and ZJH1, ZJH2, ZJH3 and ZJH4 from Hangzhou City forming the other lineage (Figs. 5C and 5D).

Genetic differentiation and genetic variation of G. kiangnanense

Population subdivisions into the Zhejiang (ZJ) and Anhui (AH) groups explained 36.59% of the total variation in AMOVA analyses, while 47.84% of the variation was within individuals. The overall F_ST value was 0.521, indicating high genetic differentiation among groups (Table 3). We also examined the genetic differentiation among populations located in Zhejiang Province, which were further divided into two lineages, ZJZ and ZJH. AMOVA results showed an overall F_ST value of 0.301, with 13.68% of the variation existed among groups and 69.87% occurred within groups (Table 3). Accordingly, the pairwise F_ST values among populations were significant, ranging from 0.062 to 0.122 (Table 4). Furthermore, the gene flow between the two lineages of G. kiangnanense was 2.792, lower than the gene flow between the two sub-populations, ZJZ and ZJH in Zhejiang Province, which is 4.323.

Population divergence history

Based on results from genetic structure analysis and phylogenetic reconstruction, three hypothetical evolutionary scenarios were designed and further evaluated using DIYABC. The optimal model was determined by continuously correcting the parameter ranges, and the relative posteriori probabilities for each scenario were estimated using multivariate logistic regression methods and direct estimation methods in comparison with the observed dataset. The comparison between the posterior probabilities of the three scenarios using local linear regression showed that model 1 was the most likely scenario, with a posterior probability of 0.570 (95% CL [0.530–0.610]), which is much higher than the model 2 (0.374; 95% CL [0.335–0.414]) and the model 3 (0.054; 95% CL [0.031–0.078]).

According to the identified best scenario, the G. kiangnanense populations from Anhui Province (AH, effective population size of ~617 individuals) most likely diverged firstly from the common ancestor of the group from Zhejiang and Anhui province around 932 generations ago (95% HPD [highest posterior density] 235.000 to 1,900.000; Table S4; Fig. S1). Within Zhejiang Province, divergence between populations from Zhuji (ZJZ, effective population size of ~1,170 individuals) and from Hangzhou (ZJH, effective population size of ~825 individuals) occurred around 352 generations ago (95% HPD [highest posterior density] 58.600 to 965.000; Table S4; Fig. S1). According to previous observation on the life cycle of G. kiangnanense, the species starts to bloom in about 6 years and bear fruits in the same year. Therefore, a 6-year generation time was assumed, which places the divergence between the AH and the common ancestor of ZJZ and ZJH populations at 5.592 ka (95% HPD [1.410–11.400 ka]), and the divergence between the ZJZ and ZJH populations at 2.112 ka (95% HPD [0.351–5.790 ka]). The changes in the number of effective populations in AH, ZJZ and ZJH may indicate an expansion trend of G. kiangnanense.

Genetic diversity of G. kiangnanense

Genetic diversity is an important indicator of a species’ ability to adapt to environmental changes and is considered as important reference for the development of conservation strategies for endangered plants. This study sampled eight G. kiangnanense populations with a total of 143 individuals across the entire distribution range and obtained molecular data from 21 microsatellite polymorphic loci to examine the genetic diversity of G. kiangnanense. Cervus analysis showed that the polymorphic information content (PIC) of G. kiangnanense ranged from 0.007 to 0.721 with a mean value of 0.236, indicating the polymorphic nature of these loci. The average number of alleles (N_A) was 1.714, which was close to the average number of effective alleles (N_E = 1.330), indicating that the alleles of the 21 loci used in this study were evenly distributed.

The expected heterozygosity (H_E) values of the eight populations of G. kiangnanense ranged from 0.143 to 0.240 with a mean value of 0.179, indicating a relatively low genetic diversity. However, genetic diversity at the species level was relatively high (H_E = 0.266). This result is consistent with Lee, Yeon & Shim (2016), who reported that seven populations of G. microrrhynchum distributed in Manchuria of Korean Peninsula had low genetic diversity and relatively high genetic diversity at the species level. The study of five G. scipetarum populations distributed in Italy using RAPD genetic markers also showed relatively low levels of genetic diversity (Marzario et al., 2022). The average H_E values of G. kiangnanense, G. microrrhynchum, and G. scipetarum subsp. eddae populations were very close (G. microrrhynchum = 0.188; G. eddae = 0.187; G. kiangnanense = 0.179).

The genetic diversity of plant populations is influenced by factors such as the extent of regional distribution, age and size, breeding system, seed dispersal, and life history (Hamrick, Godt & Sherman-Broyles, 1992). Similar to G. kiangnanense, some other endangered plants exclusively or partially distributed in Zhejiang Province, such as Sinojackia microcarpa C.T. Chen & G.Y. Li and Heptacodium miconioides Rehder, also exhibit low levels of population genetic diversity (Zhong et al., 2019; Li & Jin, 2006). A low value of Shannon information index (I) indicated lower inter-population biodiversity and the geographical proximity of the eight G. kiangnanense populations. Additionally, Yu et al. (2021) observed a declining age structure in all populations of G. kiangnanense, indicating a decline in the breeding system. G. kiangnanense is heterogamous, depending on pollinating insects for reproduction, yet it exhibits a low seed setting rate and weak seed dispersal ability. Consequently, the species faces challenges in reproduction due to its narrow regional distribution, limited seed dispersal distance, and potentially increasing anthropogenic disturbances. These factors collectively pose serious threats to the survival of G. kiangnanense, thereby impacting species diversity.

Genetic differentiation and structure of G. kiangnanense

Results from AMOVA analysis revealed a high level of population differentiation and a significant population structure among all the populations (F_ST = 0.521). The mean value of gene flow (N_m) among populations was 2.559. Gymnospermium kiangnanense is an entomophilous plant, primarily relying on insects such as Apis cerana Fabricius, Episyrphus balteatus De Geer, Scaeva pyrastri Linnaeus, and Halictus sp. for pollination (Yu et al., 2021). Considering the abundance of pollen grains and the ease of pollen dispersal, gene flow via pollen may be more prevalent than that via seed, resulting in a relatively lower ratio of genetic variation among populations within groups and a low degree of habitat fragmentation.

The genetic structure of a species is crucial for evaluating its genetic resources. Analysis of genetic differentiation patterns can provide valuable insights for the effective conservation and utilization of these resources. Our analyses, including PCoA, UPGMA, and STRUCTURE, consistently identified two major genetic groups within G. kiangnanense. The eastern populations in Zhejiang Province formed one distinct lineage (ZJ), while the western populations in Anhui Province comprised the other lineage (AH) (Figs. 3–5). The hierarchical AMOVA analysis indicated a high level of genetic differentiation (F_ST = 0.521), with 36.59% of the variation occurs between populations in the Zhejiang (ZJ) and Anhui (AH) groups (Table 3). Within Zhejiang Province, the six G. kiangnanense populations were further divided into two lineages. The southeast populations in Zhuji City formed one lineage (ZJZ), while the northwest populations in Hangzhou City comprised the other lineage (ZJH). The F_ST value between the ZJZ and ZJH lineages was 0.301, supporting a significant genetic structure within Zhejiang Province (Table 4).

Due to the complex effects of mountains, rivers, and climate in Southeast China, the geographic region in Southeast China displays high levels of ecosystem and species diversity. It harbors a large number of endemics, endangered, and tertiary relict species. The fruit of G. kiangnanense lacks a seed dispersal vector and relies solely on gravity for dispersal. Consequently, seed dispersal can be easily hindered by mountain slopes, streams, and rivers. Therefore, as one of the main mountain ranges in East China, the Huangshan-Tianmu Mountain Range located on the boundary of Zhejiang and Anhui Provinces may serve as a natural boundary separating populations of G. kiangnanense. Within Zhejiang Province, populations distributed in the northwestern part and the remnants of Tianmu Mountain were clustered together, while populations distributed between the Kuaiji and Longmen Mountain Range formed a different lineage (Fig. 1). The Longmen Mountain and Fuchun River may have contributed to the genetic differentiation between the two lineages (Fig. 1). Overall, the phylogeographic pattern of G. kiangnanense suggests a vicariant event that divided an ancestral stock into different sub-groups. The barrier generally conforms to the geographical boundaries caused by mountains and water systems in Zhejiang and Anhui Provinces.

Demographic history

The DIYABC analysis supported Scenario 1, indicating that the AH lineage of G. kiangnanense first separated from the common ancestor of current populations. The divergence time (t2) was predicted as 5,592 years. Such demographic events may be linked to post-Holocene climate change and human activities. Similar to G. microrrhynchum, G. kiangnanense is found in temperature-sensitive habitats. It is presumed to start migration to subtropical regions during the middle Miocene (15.680 Mya) when the global climate became cold and arid, and was thus diverged from its close relatives at around 14.670 Mya (Song et al., 2022). During the middle Holocene (8.450–2.450 ka), the weather became warm and humid, with a significant increase in summer precipitation, exceeding 0.500 mm/d (Zheng et al., 2004). These warmer and wetter conditions, coupled with gradually increased human activities, triggered the growth of terrestrial herbs and the gradual opening of the original evergreen and deciduous mixed broad-leaved forest between ca. 7.550 and 3.750 ka, as supported by archaeological evidence (Chen, Wang & Dai, 2009). Hence, the restricted distribution of G. kiangnanense in mountainous areas and its high inter-group genetic differentiation may be due to its inefficient seed dissemination, as well as its inability to cope with climatic amelioration and compete with other herbaceous plants under warm conditions and human interference. Examples of human activities affecting population structure can also be observed in Parrotia subaequalis (H. T. Chang, R. M. Hao & H. T. Wei) (Hamamelidaceae). A study of the species proposed that declined forests due to farmland expansion over the last 6,000 years have disrupted recent population connectivity (Zhang et al., 2018).

In the following 4,000 years, the climate fluctuated between warm and cold until 221 B.C. However, the impact of human activities during this period on natural vegetation has exceeded the influence of climate in some areas, even became the main driving force of vegetation change (Boivin et al., 2016). The broad-leaved forest disappeared after 3.750 ka, and anthropogenic deforestation and farmland management provided more space for herbaceous plants (Chen, Wang & Dai, 2009; Wang et al., 2023; Ke et al., 2023). It is hypothesized that the ZJZ lineage diverged from ZJH at 2.112 ka. Overall, our results indicate that climate change and human activities have significant impacts on the current genetic variation pattern of G. kiangnanense populations.