Genetic structure and designing a preliminary core collection of Zizania latifolia in China based on 12 microsatellites markers

Sample collection strategy

From 2019 to 2021, we conducted surveys and collected samples based on existing digital herbarium records (Chinese Virtual Herbarium, https://www.cvh.ac.cn) at 40 potential distribution areas of Z. latifolia (including habitats such as lakes, rivers, and swamps) across 28 provinces. Among these areas, wild Z. latifolia was found in only 25 survey areas, predominantly concentrated east of the Heihe-Tengchong line (Fig. 1). In each survey area, we selected a number of evenly distributed sampling sites according to the actual size of the area, collecting three live plants from each site, with a spacing of not less than 1,000 m between individual plants. These plants were then transplanted to the aquatic plant resource nursery at the Lushan Botanical Garden in Jiujiang, Jiangxi Province, for ex situ conservation. Once the collected plants were established and had survived, 3 to 5 healthy, young leaves were selected, placed in a sealed bag with an appropriate amount of color-changing silica gel, and stored in a refrigerator at −80 °C after complete drying for DNA extraction. Ultimately, we obtained a total of 357 wild Z. latifolia samples from 132 sampling sites across 25 populations (survey areas).

DNA extraction and PCR assays

Total genomic DNA was extracted using a plant genomic DNA kit (Tiangen Biotech, Beijing, China). DNA quality and concentration were then examined using NanoDrop ND2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Twelve microsatellites with clear bands and relatively high polymorphism were selected from 24 Z. latifolia-specific and cross-specific microsatellites after the pre-experiment, and the specifically developed markers were generally better than the cross-specific microsatellites (Quan et al., 2009; Richards et al., 2004, 2007). PCR amplification was performed following the protocol by Quan et al. (2009), and the PCR products were separated on a 6% denaturing polyacrylamide gel (Fig. S1). The fragments were then visualized by silver staining, and alleles were scored in reference to a DNA marker (100 bp ladder; Tiangen Biotech, Beijing, China). To ensure the high quality of the codominant SSR dataset, several precautions were taken, as suggested by Bonin et al. (2004). Then we manually converted the raw data into the input format of the software GENALEX 6.5 (Peakall & Smouse, 2012).

Genetic data analyses

The parameters of population genetic variation were estimated using mean allele number (N_a), mean effective allele number (A_e), Shannon information index (I), polymorphism information content (PIC), expected heterozygosity (H_e), observed heterozygosity (H_o), and Wright’s fixation index (FI) with GENALEX 6.5 and PowerMarker v3.25 (Liu & Muse, 2005; Peakall & Smouse, 2012). F-statistics were calculated to estimate the extent of genetic differentiation among populations. Analysis of Molecular Variation (AMOVA) was conducted to examine genetic variation both among and within populations. A Principal Coordinates Analysis (PCoA) was then performed to illustrate the pattern of genetic divergence. F-statistics, AMOVA, and PCoA analyses were all carried out using GENALEX 6.5 (Peakall & Smouse, 2012). The Unweighted Pair Group Method with Arithmetic Mean (UPGMA) dendrograms of populations and individuals based on Nei’s distance were respectively constructed and visualized in MEGA X (Kumar et al., 2018). The genetic structure was further explored using the Bayesian clustering algorithm implemented in STRUCTURE 2.3.4 (Pritchard, Stephens & Donnelly, 2000). The program was provided with no prior information on ancestral populations and was run 10 times for each value of K (K ranges from 1 to 25), under the admixture model with independent allele frequencies, using 1,000,000 Markov chain Monte Carlo iterations and a burn-in of 100,000 iterations. We inferred K using the ad hoc statistic ΔK (Evanno, Regnaut & Goudet, 2005). The resulting matrices of estimated cluster membership coefficients (Q) were permuted with CLUMPP (Jakobsson & Rosenberg, 2007). The final matrix for each K value was visualized with DISTRUCT (Rosenberg, 2004). Isolation by Distance (IBD) among populations was tested using the Mantel test implemented in GENALEX 6.5. The significance of IBD values was assessed using 9,999 permutations.

Design the core collection

The core collection of Z. latifolia germplasm was constructed using the following methods: First, the program CoreHunter3 (De Beukelaer, Davenport & Fack, 2018) was employed to develop independent core collections based on genotypic data. Various cutoff values, including 10%, 15%, 20%, 25%, 30%, 40%, and 50% of the initial collection, were used to design the core collections in CoreHunter3, based on Modified Rogers distance (MR) with default parameters. Next, PowerMarker v3.25 software was utilized, and core collections were constructed using the simulated annealing algorithm (SA), specifically with maximizing allelic richness (SANA) and maximizing genetic diversity (SAGD). The cutoff values were set identically to those used in CoreHunter3. Finally, the core collections extracted by each method with the 10% cutoff value were combined to form the final core germplasm of Z. latifolia.

Genetic diversity

The genetic variation among 357 wild Z. latifolia accessions was estimated using 12 SSR markers. Generally, all 12 SSR loci were polymorphic, each generating two to 12 alleles, indicating relatively low SSR variation. The mean number of alleles (N_a), mean number of effective alleles (A_e), polymorphism information content (PIC), Shannon information index (I), observed heterozygosity (H_o), and expected heterozygosity (H_e) for each SSR locus are listed in Table S1. The values for heterozygosity varied among these loci (H_o = 0–0.101, H_e = 0.028–0.373). The Z. latifolia populations exhibited a relatively high level of mean genetic diversity, with N_a = 1.823, A_e = 1.443, I = 0.345, H_o = 0.012, and H_e = 0.207.

The sampling information, population size, and genetic diversity parameters for each population are listed in Table 1. The overall genetic diversity of wild Z. latifolia accessions used in this study was high, with N_a = 5.000, A_e = 1.963, H_o = 0.013, and H_e = 0.439. However, three populations (SH, HZ, and CD) showed no SSR variability, and all samples within these populations were genetically uniform. The unexpectedly high mean fixation index (FI = 0.976, Table 1) indicated significant heterozygosity deficiency, although these values varied among populations.

Genetic structure

The F-statistics indicated high genetic differentiation among populations, with an overall genetic differentiation coefficient F_st = 0.452 (Table S1). This was further confirmed by AMOVA, which revealed that 60.35% of the diversity occurred within individuals, while 37.11% was attributed to genetic differentiation among groups (Table S2).

A model-based Bayesian approach was employed to infer the genetic structure of individuals under K = 2, which resulted in the best fit to the model as estimated by Evanno’s ΔK statistics (Fig. S2). According to the genetic assignment analysis of STRUCTURE, when K = 2, the populations of wild Z. latifolia could be divided into two genetic groups. Populations from high latitudes (QQHE, JMS, TH, SY, BJ, ZZ, and HZH) were clustered into Group I (North group), while the remaining populations were clustered into Group II (South group) (Fig. 2). This pattern was also supported by the phylogenetic trees (Fig. 3; Fig. S3). We also demonstrated the genetic divergence pattern illustrated by STRUCTURE when K = 3, where intermediate populations between the north and south groups were further clustered, showing genetic subdivisions within genetic groups (Fig. 2). The PCoA generated a similar pattern of population genetic structure. The first two components of PCoA explained a total variation of 39.52% (Fig. 4). The first principal coordinate accounted for 31.73% of the total variation and allowed the discrimination of the north and south genetic groups. The second coordinate, which accounted for 7.79% of the total variation, showed no clear differentiation pattern. The Mantel test provided significant evidence for isolation by distance (IBD) in Z. latifolia (R² = 0.14, p < 0.01) (Fig. 5).

Development of a preliminary core collection

The CoreHunter3 software was utilized by setting the weights of MR to 0.5 each and applying sampling proportions of 10%, 15%, 20%, 25%, 30%, 40%, and 50%. This resulted in core germplasm samples of 36, 54, 71, 89, 107, 143, and 179, respectively. The sampling proportions in the SA method, which includes the SANA and SAGD methods, were kept consistent with those in CoreHunter3.

The diversity parameters for each sampled core collection are presented in Table 2. For SANA, the H_e peaked at 0.444 when the core germplasm number reached 89 (25%), while the genetic diversity under SAGD increased linearly, peaking at 0.426 when the sample size was at its maximum. In contrast, the highest genetic diversity (0.582) was observed when the sampling proportion was set to 10% in CoreHunter3.

The core germplasms generated by different methods showed minimal overlap (Fig. 6). After removing redundant individuals, we obtained a merged core collection containing 92 accessions. Table 3 illustrates that there were no significant differences in genetic diversity parameters between the 92 core accessions and the original 357 accessions. Figure 3 also provides insight into the distribution of individuals as determined by principal coordinates analysis, showing that the 92 core accessions were evenly distributed across the entire set of 357 samples, indicating effective preservation of genetic diversity.

The distribution and habitat changes of Z. latifolia in China

According to the results of our nationwide survey of wild Z. latifolia populations during the 3-year period from 2019 to 2021, no wild Z. latifolia populations were found in 15 of the 40 survey areas that were pre-set based on the herbarium’s historical distribution data. These areas are mainly located in northern and northwestern China. In contrast, wild Z. latifolia populations were found in the remaining 25 survey areas. Currently, wild Z. latifolia is concentrated in the Yangtze River Basin, and some populations also exist in the Haihe River Basin and Northeast China. However, wild Z. latifolia populations in areas south of the Yangtze River Basin are difficult to find, which is consistent with the findings of Zhao et al. (2018). Moreover, we found large habitat loss in less than 5 years when we revisited some populations, mainly due to the destruction of wetlands during urbanization. In fact, in recent decades, a large number of freshwater wetlands have been reclaimed and their areas drastically reduced due to factors such as population growth and a lack of understanding of the ecological functions of wetlands. As a result, suitable habitats for wild Z. latifolia have been lost. According to incomplete statistics, since the end of the 1940s, more than one-third of the lakes in the middle and lower reaches of the Yangtze River have been reclaimed, with a total area of more than 13,000 km², and over 1,000 lakes have disappeared due to reclamation. About two-thirds of the lakes in the middle reaches of the Yangtze River have vanished due to reclamation, with the area of Dongting Lake shrinking from 4,350 to 2,625 km², and the area of Poyang Lake decreasing from 5,200 to 2,933 km² (Yang et al., 2010). Additionally, the blockage of the water system by water conservancy projects, such as the construction of artificial gates and embankments, may further limit the natural expansion of wild Z. latifolia populations. For example, the middle and lower reaches of the Yangtze River, which harbor the largest and most concentrated wild populations of Z. latifolia, have only three lakes—Dongting Lake, Poyang Lake, and Shijiu Lake—that remain naturally connected to the river.

Genetic diversity of wild Z. latifolia populations

Genetic diversity is an important component of biological diversity, as richer genetic diversity within a species enhances its potential ability to adapt to environmental changes (Hughes et al., 2008). As a perennial aquatic plant, Z. latifolia exhibits strong clonal growth (Xu et al., 2008; Chen et al., 2012). We observed this characteristic during the sampling process and took necessary measures to avoid collecting too many clones. The phylogeny of the genus Zizania suggests that the ancestors of Z. latifolia originated in the Americas and migrated to the Far East and East Asia through the Bering Strait during the Tertiary period (Xu et al., 2010). Xu et al. (2008) investigated the phylogeny and genetic structure of Z. latifolia in China based on the nuclear gene sequence Adh1a, finding that populations at high latitudes (Northeast China) possessed a richer range of haplotypes. Xu et al. (2015) also used three pairs of universal microsatellite markers to compare the genetic diversity of four species within the genus Zizania. They found that the genetic diversity of Z. latifolia was relatively low (H_e = 0.374) compared to Z. palustris (H_e = 0.630). However, this work may have underestimated genetic diversity due to the use of too few molecular markers. Subsequently, Zhao et al. (2019) and Wagutu et al. (2022) examined the genetic diversity and genetic structure of wild populations of Z. latifolia in China using SSR molecular markers. Both studies found high genetic diversity (mean diversity: H_e = 0.426 in Zhao et al., 2019, H_e = 0.173 in Wagutu et al., 2022; overall diversity: H_e = 0.432 in Zhao et al., 2019, H_e = 0.241 in Wagutu et al., 2022). In comparison, the present study found a similar level of overall genetic diversity (H_e = 0.439) as Zhao et al. (2019), but higher than Wagutu et al. (2022); and the mean diversity (H_e = 0.207) was at the same level as Wagutu. These findings may be due to sampling differences between studies. Although Wagutu et al. (2022) surveyed more populations, their sampling was limited to watersheds on five latitudinal gradients. Populations within the watersheds may trend to be genetically identical due to the presence of gene flow, thus making relatively low overall genetic diversity. In contrast, both Zhao et al. (2019) and our work sampled populations extensively across the country. Unfortunately, previous studies (Xu et al., 2015; Zhao et al., 2019; Wagutu et al., 2022) did not calculate the genetic parameters of I and PIC, preventing us from directly assessing differences in genotypic diversity across studies. However, the very low observed heterozygosity (H_o) coupled with an unusually high fixation index (FI) suggests a heterozygote deficiency in most populations. To collect more representative genotypes for constructing a germplasm repository, we sampled an area (population) with multiple sample sites spread far apart to cover as much of the area as possible, spacing samples approximately 1,000 m apart within each sample point. This sampling strategy may have two potential consequences: First, multiple populations might have been mistakenly sampled as one, leading to genetic structure subdivision and high overall genetic diversity (Wahlund effect). Second, if some wild populations of Z. latifolia are already experiencing heterozygote deficiency (Zhao et al., 2019), the sign of heterozygote deficiency could be further amplified because of the higher likelihood of collecting homozygous individuals. The increased probability of biparental inbreeding due to extensive clonal growth in Z. latifolia may also raise the proportion of homozygous genotypes in the population. Additionally, heterozygote deficiency may result from the insufficient polymorphism and number of molecular markers used in this study. For example, Wagutu et al. (2022) used 46 pairs of SSR markers and found heterozygote excess (FI = −0.245), whereas only 12 pairs of SSR markers were used in our study. Although conventional gel electrophoresis may underestimate the actual number of alleles, the relatively low PIC (mean PIC = 0.383) also suggests that this earlier batch of developed SSR markers used in our study may have insufficient polymorphism and that new high-polymorphic SSR markers need to be developed using the genome or transcriptome data (Quan et al., 2009; Wagutu et al., 2022; Wang et al., 2023). Now that a high-quality chromosome-level genome of Z. latifolia has been released (Yan et al., 2022), subsequent studies can employ whole-genome resequencing to more accurately evaluate the genetic diversity of Z. latifolia (Zou et al., 2024).

Genetic population structure of wild Z. latifolia

The results of this study also showed high genetic differentiation (F_st = 0.452) among Z. latifolia populations in China, suggesting restricted gene flow among these populations. Specifically, we found significant isolation by distance (IBD) through the Mantel test, and STRUCTURE, PCoA, and UPGMA clustering analyses indicated that the wild Z. latifolia populations can be divided into north and south genetic groups. The natural habitat of Z. latifolia is generally fragmented, forming isolated habitats that limit gene exchange between populations. This aligns with the observation that aquatic plants typically exhibit high genetic differentiation (Barrett, Eckert & Husband, 1993). Even though these habitats may sometimes achieve intermittent connectivity through water regimes, this connectivity does not necessarily coincide with the seed maturity period required for effective gene flow via waterborne seed dispersal. Additionally, whether birds are significant mediators of seed dispersal in Z. latifolia remains unclear, and pollen flow is only effective in localized areas. Overall, Z. latifolia lacks the ability to disperse gene flow over long distances, resulting in isolation by distance on a large scale (Willson, 1983). Furthermore, the genetic structure of aquatic plants is influenced not only by geographic isolation but also by the connectivity of water systems (Wang et al., 2008; Honnay et al., 2010; Davis et al., 2018). Both Zhao et al. (2019) and Wagutu et al. (2022) have previously reported a north-south genetic differentiation in wild populations of Z. latifolia, and our results also support this conclusion. However, Wagutu et al. (2022) concluded that populations across different latitudinal gradients were more genetically similar within watersheds, whereas genetic subdivision occurred between watersheds in a latitudinal direction. In contrast, Zhao et al. (2019) did not observe latitudinal genetic subdivision, instead pointing out genetic subdivision in an east-west direction within the Yangtze River Basin. To further explore these patterns, we analyzed the genetic structure of Z. latifolia populations using STRUCTURE software by applying K = 3 (Fig. 2). However, we did not find any clear genetic subdivisions associated with the geographic distribution of the populations. Interestingly, in both STRUCTURE and UPGMA analyses, the DL populations located in high northern latitudes were clustered with populations in the south, raising questions about the origin of this population that should be addressed in subsequent work. Zhao et al. (2019) has reported the existence of Z. latifolia populations formed by the escape of cultivars into natural habitats. If feral Z. latifolia samples were incorrectly collected, this could lead to misclassification of the genetic structure of the populations.

The designing of a preliminary core collection for wild Z. latifolia

Core germplasm is a key component in plant breeding that can increase the effectiveness of utilizing plant genetic resources for genetic improvement and new variety development. The purpose of core germplasm is to maximize the representation of the genetic diversity of the entire species with a minimum number of genetic resources. To accomplishing this goal, the sampling strategy is important, and it determined the results of core germplasm construction. A suitable sampling plan can optimize the preservation of genetic diversity while minimizing genetic redundancy. Of these, stratified sampling is often considered superior to random sampling in many aspects, as demonstrated by numerous studies (Brown, 1989; Casler, 1995; Escribano, Viruel & Hormaza, 2008; Kumar et al., 2020). Stratified sampling can be implemented using various methods, including proportional, logit, square root, and genetic diversity-based sampling, each with its own advantages and disadvantages. When selecting a sampling strategy, it is important to consider factors such as the amount of germplasm material, grouping, and genetic diversity.

Core germplasm construction involves extracting a portion of the total samples from existing germplasm resources that meet specific criteria. Achieving this goal hinges on selecting appropriate sampling methods and determining suitable sampling ratios. Currently, the commonly used software for constructing core germplasm includes PowerMarker, PowerCore, and CoreHunter. The proportion of core germplasm relative to the original materials in various plants, both domestically and internationally, ranges from 5% to 40% (Wang et al., 2007; Escribano, Viruel & Hormaza, 2008; Liang et al., 2015; Dervishi et al., 2021). However, none of these studies provided a uniform sampling rate, primarily due to differences in sampling completeness, levels of genetic diversity, and the biological characteristics specific to each plant. In our research, we used PowerMarker v3.25 (SANA and SAGD) and CoreHunter3 (MR) to construct a preliminary core collection of wild Z. latifolia. Our results revealed very limited intersection between the core germplasm constructed by different methods using a 10% sampling ratio (Fig. 6), likely due to differences in the algorithms underlying each method. As a result, we combined the core germplasm constructed by each of the three methods. Evaluating the representativeness of core samples in relation to the diversity of the entire germplasm resource is also a key focus in core germplasm studies. Referring to previous studies, we first created a principal component map based on the final core germplasm and compared it with the principal component distribution maps of all materials, finding a high degree of agreement (Fig. 4). Additionally, no significant differences were observed between the genetic diversity indices of the final core germplasm and the original germplasm materials (Table 3). The results demonstrated that the final core germplasm was well-matched with the original germplasm, indicating the rationality of the final core germplasm screening results. This core collection will facilitate the utilization of Z. latifolia in developing new crop varieties with enhanced yield and resistance traits.