Phylogeography of the rare and endangered lycophyte Isoetes yunguiensis

Materials

After obtaining approval from the Guizhou Province Wild Flora and Fauna Management Station, a total of 167 samples were collected from 14 populations in Yunnan and Guizhou Provinces from October 2018 to January 2019 (Table 1). Briefly, 8–15 samples that were distributed 10 m apart were randomly selected from each site (random sampling of minimal population). Healthy young leaves were collected in paper bags, dried immediately with silica, and stored at −20 °C for further use. Voucher specimens were stored in the Guizhou Agricultural College Plant Protection Department Herbarium (GACP). Figure 1 presents a map of sampling locations and geographical distribution of haplotypes.

DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing

DNA was extracted from dried leaves using the TIANGEN Plant Genomic DNA Extraction Kit [DP305, TIANGEN BIOTECH (BEIJING) CO., LTD.]. PCR analyses were performed using the EST-SSR primers SJSSR5, SJSSR12, SJSSR14, SJSSR45, and SJSSR53 and the cpDNA primers trnS-trnG, psbC-trnS, and psbD-trnT (Table 2). PCR reactions were conducted in 25-µL mixtures comprising 12.5 µL 2 × T5 Super PCR Mix (Beijing TsingKe Biotech Co., Ltd.), 40 ng DNA template, and 1 µL each of (10 pmol/µL) forward and reverse primers, topped-up to final volume with double-distilled water. EST-SSR PCR was performed with an initial denaturation step at 94 °C for 5 min, followed by 28 cycles of denaturation at 94 °C for 30s, annealing at respective temperatures for 30 s, extension at 72 °C for 1 min, and final extension at 72 °C for 30 min. Reactions with cpDNA primers were performed with an initial denaturation step at 94 ° C for 5 min, followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at respective temperatures for 30 s, extension at 72 °C for 1 min, and final extension at 72 °C for 10 min. EST-SSR PCR products were visualized as clear bands using 1.5% agarose gel electrophoresis and detected using capillary electrophoresis with ABI 3730 (Thermo Fisher Scientific). GeneMaker 2.2 was used to analyze the original peak patterns and to determine the sizes of EST-SSR marker fragments. Similarly, cpDNA PCR products were visualized as clear bands using 1.5% agarose gel electrophoresis, purified via electrophoretic cutting and using magnetic beads, and sequenced using ABI 3730. All sequences have been deposited in GenBank under the accession numbers MN463102 –MN463268 for psbC-trnS, MN463269 –MN463435 for psbD-trnT, and MN463436 –MN463602 for trnG-trnS.

DNA analysis

Observed allele number (Na), effective allele number (Ne), observed heterozygosity (Ho), and expected heterozygosity (He) were calculated, and principal coordinate analyses (PCoA) and Mantel tests (Mantel, 1967) were performed using GenAlEx 6.5 (Peakall & Smouse, 2012). Polymorphism information content (PIC) was calculated using CERVUS 3.0 (Kalinowski, Taper & Marshall, 2007). Nei’s genetic distances were calculated using POPGene 32 (Yeh et al., 2000), and Bayesian clustering was performed using STRUCTURE 2.3.3 (Pritchard, Wen & Falush, 2009). The original peak images of noncoding cpDNA regions were manually evaluated using Chromas 2.6, whereas the ClustalW program in the MEGA 7.0 software was used to align and truncate the three gene fragments and examine whether poly structure existed or were deleted. Concatenate Sequence program in PhyloSuite v1.1.15 (Zhang et al., 2018) was used to connect the three gene fragments into one fragment. Insertion–deletions are considered mutation sites (Chen et al., 2008). DnaSP v5 (Librado & Rozas, 2009) was used to determine variation sites, haplotype diversity (Hd), and nucleotide diversity (π) of cpDNA sequences, and neutral tests and mismatch distribution analyses of cpDNA fragments were performed at the species level. LAMARC 2.1.8 was used to further verify the group expansion events for I. yunguiensis, and Permut 2.0 (Pons & Petit, 1996) was used to calculate total genetic diversity (Ht) and intrapopulation mean genetic diversity (Hs). To construct haplotype phylogenetic trees, analysis of molecular variance (AMOVA) was performed using Arlequin 3.5 (Excoffier & Lischer, 2010) and IQ-TREE (Nguyen et al., 2014). A median-joining network map was constructed using Network 5.0 (Bandelt, Forster & Röhl, 1999).

Species distribution prediction

MaxEnt 3.4 (Phillips, Anderson & Schapire, 2006) was used to predict the distribution ranges of I. yunguiensis during the last interglacial (LIG) period and last glacial maximum (LGM) as well as in the present and future. Coordinates of I. yunguiensis populations were determined using field records from the China Digital Herbarium (http://www.cvh.ac.cn/), as described previously (Pang et al., 2003; Yuan et al., 2012; Li et al., 2015). Maps were acquired from the National Catalogue Service for Geographic Information (http://www.webmap.cn). Climatic data were downloaded from the World Climate Data website (http://www.worldclim.org; resolution choice, 30 s). Correlations among 19 bioclimatic variables were calculated with the Raster Correlations and Summary Statistics tool in the SDMtoolbox. The environmental and extreme climatic factors with a correlation coefficient of —r— ≤ 0.8 were selected. Finally, eight environmental variables were selected as the prediction variables of the distribution of I. yunguiensis: average daily range (bio 2), temperature seasonality (bio 4), maximum temperature of the warmest month (bio 5), average temperature of the wettest quarter (bio 8), average temperature of the coldest quarter (bio 11), precipitation of the wettest month (bio 13), precipitation of the driest month (bio 14), and precipitation of the warmest quarter (bio 18). Briefly, environmental data and spatial coordinates were imported into MaxEnt 3.4, parameters were set to default values, and distribution ranges during various periods were predicted. Subsequently, ArcGis 10.2 was used to generate ASCII raster layers and simulated distribution maps.

Genetic diversity of EST-SSR loci

The genetic diversity of EST-SSR loci in I. yunguiensis populations is presented in Table 3. Mean Na was 2.843, ranging from 1.929 (SJSSR53) to 4.000 (SJSSR12), mean Ne was 2.138, ranging from 1.521 (SJSSR53) to 2.981 (SJSSR12), and mean Shannon information index (I) was 0.782, ranging from 0.417 (SJSSR53) to 1.175 (SJSSR12). Ho ranged from 0 (SJSSR14) to 0.960 (SJSSR5), and He ranged from 0.273 (SJSSR53) to 0.642 (SJSSR12). Mean Ho and He were 0.475 and 0.463, respectively. PIC was high at a mean of 0.516, ranging from 0.313 (SJSSR53) to 0.641 (SJSSR12). The genetic diversity was the highest for SJSSR12 loci and the lowest for SJSSR53 loci.

Genetic diversity and structure of populations

The genetic diversity of different populations is summarized in Table 4. Mean Na was 2.843, ranging from 2.400 (LS) to 3.400 (CX and XB), mean Ne was 2.138, ranging from 1.777 (LS) to 2.675 (CX), and mean I was 0.782, ranging from 0.595 (LS) to 1.016 (CX). Ho and He ranged from 0.251 (XB) to 0.673 (DJC) and from 0.335 (XLC) to 0.591 (CX), with mean values of 0.475 and 0.463, respectively. The genetic diversity was high in BSH, CX, DJC, DPJ, GJS, SG, and PB and low in LS, GP, and XB.

AMOVA using EST-SSR markers showed that 22.25% of the total genetic variation was present between I. yunguiensis populations and up to 77.75% was present within populations (Table 5). Mean interpopulation gene flow rate (Nm) was 1.534, which is considered high, indicating frequent gene exchange between populations. These observations are consistent with the lower genetic variation observed between populations than within populations.

Mantel tests showed no significant associations of Nei’s genetic distances with geographical distances (p = 0.242 > 0.05; Fig. 2A) or altitude gradients (p = 0.402 > 0.05; Fig. 2B).

Genetic relationships and cluster analyses

Unweighted pair group method with arithmetic mean (UPGMA) cluster analyses were performed based on Nei’s genetic distances using MEGA 7.0, and 14 UPGMA phylogenetic trees of I. yunguiensis populations were constructed (Fig. 3). With a genetic distance of 0.2, the 14 populations were divided into 2 groups. The first group CX comprised only one unique population from the Yunnan Plateau, whereas the second group comprised the remaining 13 populations from the Guizhou Plateau. With a genetic distance of 0.1, the second group was subdivided into the groups II-a and II-b. The group II-a predominantly comprised populations from the west of Guiyang City, whereas the group II-b comprised populations from to the east of Guiyang City; however, there was no notable difference between both groups. PCoA results (Fig. 4) were consistent with UPGMA results, in which the first and second components explained 51.50% and 18.76% of the variance, respectively.

STRUCTURE 2.3.3 was used to analyze the genetic structure of I. yunguiensis populations. At K = 2, ΔK reached the maximum value (Fig. 5), and I. yunguiensis individuals were divided into two groups (Fig. 6). The group I (blue; Fig. 6) mainly comprised individuals of the CX population, whereas the group II comprised individuals from the remaining populations (red; Fig. 6). Overall, 30 of the 167 individuals showed differences of different degrees, indicating a variable degree of gene exchange between populations.

cpDNA sequence characteristics and haplotype analyses

After aligning and connecting the three cpDNA sequences from I. yunguiensis, the total length was 1,657 bp. Moreover, 36 mutation sites were detected, among which 11 were detected by trnS-trnG, including 8 single-base mutations and 3 insertion–deletion sites. Only two mutation sites were detected by psbC-trnS, and both were single-base mutations. Overall, 23 mutation sites were detected by spsbD-trnT, including 20 single-base mutations and 3 insertion–deletions.

A total of 21 haplotypes were detected using DnaSP v5, and haplotype H1 was the most widely distributed haplotype. With the exception of the CX population, in all the remaining 13 populations, haplotype H1 was the most common, followed by haplotype H2. Haplotype H3 was distributed in GP, BSH, DPJ, LS, and PB populations, haplotype H4 was distributed in GP and SCZ populations, and the unique haplotypes H5–H21 were distributed in only a single population, and only 1 CX population carrying the unique haplotype H7. Genetic diversity, nucleotide diversity, and haplotype compositions and frequencies of the sampled populations are shown in Table 6.

Population genetic structure based on cpDNA sequences

The total genetic diversity (Ht) was 0.631, mean genetic diversity (Hs) was 0.515, and the genetic differentiation coefficients between (GST) and within (NST) populations were 0.183 and 0.344, respectively. The haplotype variation structure of the 14 I. yunguiensis populations was examined using 1000 bootstrap replicates, which revealed a tendency of NST >  GST, although the difference was not significant (p > 0.05); therefore, I. yunguiensis exhibited no distinct phylogenetic structure in the studied areas. Consistently, AMOVA (Table 7) showed that the genetic variation between populations was 34.00% and that within populations was 66.00% (FST = 0.34, p < 0.001, 1000 bootstrap replicates). The genetic variation in I. yunguiensis was predominately higher within populations. Mean Nm was 2.18, which indicated substantial gene flow between I. yunguiensis populations.

Population history dynamics based on cpDNA sequences

DnaSP v5 was used to examine the combined cpDNA sequence of all the individuals in the 14 sampled I. yunguiensis populations. Although the mismatch distribution was bimodal (Fig. 7), Tajima’s D = −2.55856 (p < 0.001), Fu & Li’s D = −8.03665 (p < 0.02), and Fu & Li’s F = −6.94682 (p < 0.02) were significantly negative. The population growth index obtained by LAMARC simulation was 921.590 (G > 200). Taken together, these results indicated that the I. yunguiensis populations have experienced expansion events.

Phylogenetic relationships based on cpDNA haplotypes

IQ-TREE was used to construct a phylogenetic tree of I. yunguiensis haplotypes (Fig. 8). The model with the best HKY+F value was selected. Isoetes flaccida and I. vallda were used as outgroups. The 21 haplotypes were divided into two branches: the first branch included H1–H6 and H8–H21 and the second branch included H7 alone.

A mediation-link (median-joining) network diagram and phylogenetic tree based on the maximum likelihood method were constructed using Network 5.0 (Fig. 9). Because H1 and H2 appeared in the middle of the Network diagram, these were speculated to be relatively old haplotypes; moreover, there were few other haplotypes as well as variation coefficients were small, suggesting that recent haplotypes were derived from H1 and H2.

Niche modeling

MaxEnt 3.4 was used to predict distribution ranges during the LIG period (Fig. 10A) and LGM (Fig. 10B) as well as in the present (Fig. 9C) and future (Fig. 10D). In all models, the areas under the curve were >0.947, indicating that the predicted distribution ranges of I. yunguiensis are located in the central Guizhou Plateau and northwestern and eastern Yunnan Plateau and that the suitable areas of distribution are mainly located in the central Guizhou Plateau. The future distribution pattern was similar to the present distribution pattern. The distribution range of I. yunguiensis during the LIG period was significantly reduced compared with the present distribution range, and the distribution ranges in the Yunnan–Guizhou Plateau were clearly separated. The distribution range during LGM was slightly expanded in the Yunnan Plateau and reduced in the Guizhou Plateau compared with the future distribution range.

Genetic diversity and structure

In the present study, I. yunguiensis showed high genetic diversity compared with other rare, endemic, and endangered plants, such as Isoetes sinensis Palmer (He = 0.118; Kang, Ye & Huang, 2005), Isoetes hypsophila Hand. Mazz. (He = 0.039; Chen et al., 2010), Isoetes malinverniana Ces. & De Not. (H = 0.1491 for ISSR data; H = 0.2289 for AFLP data; Gentili et al., 2010), Ottelia acuminata var. jingxiensis H. Q. Wang & S. C. Sun J. (He = 0.441; Li et al., 2019), Dalbergia odorifera T. Chen (He = 0.37; Liu et al., 2019), and Brasenia schreberi J. F. Gmelin (He = 0.256; Li et al., 2018). Most rare and endangered plants are seed plants that can achieve a higher gene flow via long-distance pollen and seed dispersal than ferns that rely solely on spores. Ma et al. (2018) investigated the genetic structure of I. yunguiensis using inter-SSR markers and showed inter- and intrapopulation genetic variations of 31.99% and 68.01%, respectively. Similarly, Dong et al. (2018) used amplified fragment length polymorphism markers and showed inter- and intrapopulation variations of 40.12% and 59.88%, respectively. However, in these previous studies, samples were collected from Pingba and Hongfeng Lake alone, which limited the representativeness of the genetic variation of I. yunguiensis. The 14 sampled populations covered all areas of I. yunguiensis distribution, and the analyses of genetic structures and variations were conducted using EST-SSR and cpDNA sequences to reveal the genetic diversity of I. yunguiensis. Based on these markers, I. yunguiensis showed high genetic diversity. Moreover, AMOVA using these markers showed that the genetic diversity was higher within populations than between populations, which is consistent with the results of previous studies. Reportedly, in the areas of I. yunguiensis distribution, the observed haplotypes accounted for 80.95% of the total haplotypes, and this species often did not show clear biogeographical patterns (Gao et al., 2012; Gao et al., 2016). Accordingly, in the present study, the population genetic differentiation index was higher than the gene differentiation index, although the difference was not significant (p > 0.05), indicating that I. yunguiensis shows no distinct biogeographical patterns. I. yunguiensis is mainly dispersed through spores released in flowing water (Yang et al., 2011), and it is distributed in the Yangtze and Pearl River water systems. In the present study, a high gene flow was observed between I. yunguiensis populations, which may explain the high genetic variation within populations (Huang et al., 2017).

Using STRUCTURE 2.3.3, the 14 I. yunguiensis populations were divided into two groups, and these results were consistent with those of UPGMA culture analysis, PCoA, and IQ-TREE based on Nei’s genetic distances. The two groups are separated by the Wumeng Mountains, which follow the border between the Yunnan–Guizhou Plateaus. The mountains between these provinces have isolated these I. yunguiensis populations, thereby explaining the genetic distances between the two main haplotype groups. However, based on the STRUCTURE analysis, the genetic structure of I. yunguiensis is relatively complex. Further studies are required to elucidate specific reasons underlying the difference between the CX population from the Yunnan Plateau and the remaining 13 populations from the Guizhou Plateau. However, Mantel tests showed no significant associations of Nei’s genetic distances with geographical distances or altitude gradients.

Population structure and distribution ranges of I. yunguiensis

The genetic structures of extant plants can be traced back to the Quaternary glacial period. Based on this, Yu et al. (2013) have summarized the genetic structures of 36 species of alpine plants from the Qinghai–Tibet Plateau and surrounding areas. Some hardy plants did not migrate to lower altitudes during the glacial episodes; however, plants on the plateau surface experienced small-scale expansions of ranges following LGM. The genetic structures of species are often characterized by high variations or specific haplotypes, with population distributions in the highland mesa of one or more isolated areas. In addition, some species occupied narrow ranges during the glacial period, leading to small areas of expansion following LGM. These populations were characterized by high genetic diversity and unique haplotypes distributed evenly throughout the distribution range. The present study demonstrated that although I. yunguiensis shows high genetic variation, regional genetic variation that is not linked to longitude and latitude or elevation is present, presumably because during the Quaternary glacial period, the Yunnan–Guizhou Plateau remained free of ice. Therefore, the influence of glacial episodes on I. yunguiensis was only characterized by the influences of changes in climate and availability of regional water, which affect the habitat and dispersal of I. yunguiensis. Few ancient haplotypes may have randomly been fixed in the population during migration, and young haplotypes could have subsequently been fixed in the population because of genetic drift and the founder effect, thereby forming the current population structure of I. yunguiensis.

Hewitt (2000) identified a small population that migrated to form populations with greater genetic variation and richer haploid types, potentially reflecting the duration of isolation and accumulation of genetic variants. When diffusion distances increase from outer diffusion migration groups that are prone to genetic drift or the founder effect, haplotype polymorphism is gradually reduced. According to the coalescent theory, populations located in refuge show high genetic diversity and ancient haplotypes are often located in the center of the haplotype network maps. Therefore, the geographical locations of populations with high genetic diversity and ancient haplotypes can be used as bases for identifying potential habitats. In the present study, H1 and H2 were located in the center of the haplotype network map and may thus be ancient haplotypes that are widely distributed in the Guizhou Plateau. The CX population from the Yunnan Plateau comprised the unique and stable haplotype H7 and is located far from the Guizhou Plateau. It is possible that I. yunguiensis occupied specific areas of the Guizhou and Yunnan Plateaus during the glacial period. Because the ancestral population from the Yunnan Plateau is extinct, only one population was noted in the Yunnan Plateau in the present study, and the specific location of the extinct population in this plateau remains unknown.

Conservation strategies for I. yunguiensis

The present study demonstrated that the BSH, CX, DJC, DPJ, GJS, SG, and PB populations showed high genetic diversity. It is important to understand genetic diversity and population structure to establish scientific and effective protection measures. Therefore, to protect the endangered I. yunguiensis plants, a provincial nature reserve has been established; the DPJ population used in the present study was sampled from this reserve. However, it has been difficult to protect this population. The SCZ population in the present study was sampled from a farm where sewage is directly discharged into the habitat. The DJC population was sampled from the middle of Longli Forest Farm; seven populations were sampled from this farm, accounting for half of the surveyed populations. However, although these populations were not under the threat of sewerage exposure, limited attention and protection have led to a state of low survival and even extinction. Therefore, the establishment of a protected area dedicated to I. yunguiensis in Longli Forest Farm and training of staff and residents of the surrounding areas regarding protective measures for I. yunguiensis are warranted. In a previous study, populations with high genetic diversity were selected from similar environments for transplantation (Zhang, 2016). According to the principles of genetic diversity and conservation prioritization in special areas, the HFH populations, which showed the highest haplotype diversity (Hd = 0.8) should be granted the highest priority, followed by the BSH, CX, DJC, DPJ, GJS, SG, and PB populations. As an overall conservation strategy based on the molecular genetic evidence provided in the present study, the Guizhou and Yunnan Plateaus can be divided into two areas to protect populations with high genetic diversity as well as those with high and stable haplotype diversity. The CX populations from the Yunnan Plateau and the BSH, HFH, DJC, SG, DPJ, GJS, and PB populations from the Guizhou Plateau should be primarily protected. Both the Guizhou and Yunnan Plateaus should be protected, with particular focus on the CX, BSH, and HFH populations. Moreover, the water quality in I. yunguiensis habitats should be maintained to avoid contamination, and the growth of companion species should be controlled. Improved knowledge and awareness regarding the areas of I. yunguiensis distribution are paramount. To strengthen the conservation strategies for endangered species, further studies are required to identify threats. Moreover, considering that the number of wild I. yunguiensis individuals included in this survey was <10,000, the populations are declining further.