Molecular marker development and genetic diversity exploration in Medicago polymorpha

Sample collection and field experiment

In China, M. polymorpha mainly distributes in five major cultivation areas (Jiangsu province, Zhejiang province, Anhui province, Yunnan province, and Shanghai) with very limited distribution elsewhere. For SLAF-seq analysis, a total of 10 representative M. polymorpha individual accessions with distinct phenotypic characteristics were collected from these five different geographical regions. Among them, six accessions came from Jiangsu province (MP-2, MP-4, MP-5, MP-7, MP-11, and MP-13), and the other four accessions from the other four areas (MP-6 from Zhejiang, MP-8 from Shanghai, MP-3 from Anhui, and MP-9 from Yunnan). For plant growth form, MP-3 and MP-4 grown in lawn were prostrate, MP-13 grown in farmland was semi erect, and the other seven M. polymorpha accessions grown in farmland were erect (Table 1).

Each living accession was taken to the experimental fields of Yangzhou University, Jiangsu Province, China (119°40′E, 32°35′N) for ball planting. After the seeds of each accession were harvested in May, field experiments for phenotypic trait measurement were conducted in a randomized complete block design with three biological replications in October at Yangzhou University. There were two rows for each accession. The seeds of each genotype in each replication were sown in drill with a spacing of 0.5 m in rows and 1 m between two accessions. No fertilizer was used, but timely drainage between the plot ditch, artificial weeding, and timely irrigation were utilized.

Phenotype measurement and data statistical analysis

At the beginning of each accession’s flowering time in April of the following year, six different agronomic traits (leaf length (cm), leaf width (cm), plant height (cm), branch number (n), stem diameter (cm) and stem-leaf ratio) were evaluated on 10 plants per replication. Branch number is the number of branches at the base (including the main stem), and plant height is the absolute height from the base of the main stem to the growing point of the main stem. Stem diameter refers to the diameter of the base of the main stem, measured with a vernier caliper. Analysis of agronomical data was carried out using analysis of variance (ANOVA) and the general linear model (GLM) of SAS (V.9.4). Cluster analysis was conducted using the Squared Euclidean Distance method with SPSS software (V.19.0).

Extraction of genomic DNA and enzyme solution design

Young healthy leaves of each accession were collected, frozen in liquid nitrogen, and then stored at −80 °C before DNA extraction. Total high-quality genomic DNA was extracted using a minor modified CTAB method (Murray & Thompson, 1980). DNA concentration and quality were examined with an ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, USA) and by electrophoresis in 1% agarose gels with a lambda DNA standard. DNA concentrations in all samples were adjusted to the same level for the construction of SLAF libraries. To obtain adequate SLAF tags, the restriction enzyme combination “RsaI + HaeIII” (NEB, Ipswich, MA, USA) was selected based on (a) lots of SLAF tags, (b) fewer restriction fragments with repeat sequences, (c) an even distribution across chromosomes, and (d) simulated fragments aligning uniquely to the reference genome.

SLAF library construction and sequencing

SLAF-seq was performed with slight modifications as described previously (Sun et al., 2013). First, for the in-silico prediction of the number of markers produced by different enzymes, marker-discovery experiments were designed by analyzing the M. truncatula A17 (NCBI: taxid 3880; https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=info&id=3880) using prediction software independently developed by Biomarker Technologies Corporation (Beijing, China). Based on the result of the preliminary experiment, RsaI and HaeIII restriction enzymes were selected as the best combination. To further evaluate the accuracy of the ‘RsaI + HaeIII’ combination enzyme digestion experiment, Arabidopsis thaliana ecotype Columbia was used as the control for sequencing. Genomic DNA was then digested with restriction enzyme RsaI (NEB, Ipswich, MA, USA), ligated to RsaI adapter, and then digested by additional enzyme HaeIII (NEB). Polymerase chain reactions (PCR) were carried out using the diluted restriction-ligation samples, dNTP, Taq DNA polymerase (NEB), and Rsa I-primer containing a barcode. According to the position of Arabidopsis reads in the genome and the length distribution of inserted fragments in Arabidopsis thaliana ecotype Columbia, DNA fragments of 314–414 bp for M. polymorpha were selected using the Blue Pippin automated DNA size selection system. We classified these as SLAF tags and then sequenced them using the Illumina HiSeq 2500 system (Illumina, Inc., San Diego, CA, USA) with 125 bp pair-end. The sequencing was performed by Biomarker Technologies Corporation (Beijing, China).

Molecular marker development

SLAF marker identification and genotyping were performed as described by Sun et al. (2013). High-quality SLAF tags were obtained from 10 M. polymorpha genotypes. SLAFs with two to four tags were considered polymorphic SLAFs and potential markers. Polymorphic SLAF tags showed sequence polymorphisms between different accessions. Before mapping, low-quality reads with quality scores less than 30 were removed, and high-quality reads with terminal 5-bp and barcode sequences were discarded. All high-quality SLAF paired-end sample reads of the 10 accessions were then aligned to the newly assembled reference genomes of M. polymorpha ‘Huaiyang Jinhuacai’ (Cui et al., 2021; https://ngdc.cncb.ac.cn/gwh/Assembly/17540/show) using the Burrows-Wheeler alignment tool (BWA) software (Li & Durbin, 2009) with default parameters. Variants at each SLAF locus were called independently for each accession by GATK (v.4.0) (McKenna et al., 2010), and then the genotypes of all accessions (n = 10) were joined into a VCF file with the criteria of integrity >90%, minor allele frequency (MAF) >0.05, and max missing 0.5. The filtered high-quality SNPs were further used to reconstruct the phylogenetic tree with the maximum likelihood estimation method in MEGA-X (Kumar et al., 2018). The principal component analysis (PCA) of total individuals was performed using GCTA v1.93 (Yang et al., 2011). A three-dimensional (3D) point cloud was plotted using the “scatterplot3d” package in R 4.1.0 software (R Core Team, 2021) (https://CRAN.R-project.org/package=scatterplot3d).

Genetic diversity analyses

The assessments of genetic diversity, including minor allele frequency (MAF), observed heterozygosity (Ho), expected heterozygosity (He), Nei diversity index (Nei), Shannon’s information index (I), and polymorphic information content (PIC), were calculated through Perl programming.

Phenotypic characterization

There were significant differences in the phenotypic traits across the different M. polymorpha accessions. As shown in Table 2, the leaf length of M. polymorpha materials in this study varied from 1.38 to 2.40 cm, with a variation amplitude of 1.02 cm, and the leaf width varied from 1.45 to 2.23 cm, with an amplitude of 0.78 cm. Among all accessions, the leaf length and width of MP-6 from Wenling, Zhejiang province were the largest, and this material is also known as “Wenling big leaf” in the local area. We found that the germplasm characteristics of the MP-6 large leaves were not changed due to alloplanting. The leaf length and width of MP-4 from Danyang, Jiangsu province were the smallest, followed by MP-3 from Anhui, which retained more wild characteristics. The plant height of M. polymorpha accessions showed obvious differentiation. For example, the plant heights of MP-3 and MP-4 were only about 4 cm, and the lateral branches crept close to the ground because the main stem growth was significantly inhibited. However, other materials grew upright, ranging in height from 20 to 25 cm, except for MP-13. Moreover, the lateral branches of MP-3 and MP-4 were very developed and totaled up to 4, which was significantly more than that of the other materials. Stem diameter ranged from 1.71 to 2.08 cm, with an average of 1.89 cm. MP-9 from Chuxiong, Yunnan province had a stem-to-leaf ratio, of 1.33, indicating that it had the most leaves, while MP-4 had the fewest leaves. The variation coefficient showed that the plant height of the M. polymorpha accessions changed the most (43.15%), followed by branch number (33.42%), and stem-leaf ratio and stem diameter changed the least (9.60% and 4.84%, respectively). Additionally, as the correlation analysis in Table 3 shows, leaf length and leaf width were significantly positively correlated, and leaf length and leaf width were both significantly positively correlated with plant height. Branch number, stem diameter, and stem-leaf ratio were positively correlated with each other, indicating that the M. polymorpha accessions with thicker stem had greater branch numbers.

Cluster analysis of phenotypic traits

Clustering 10 M. polymorpha accessions based on their phenotypic traits (Fig. 1A) revealed three groups and a high value of genetic diversity among various accessions. The first group contained 70% of the studied accessions from all geographical regions, including MP-5, MP-6, MP-7, MP-8, MP-9, and MP-11. The second group was comprised of 10% of the accessions assessed. The third group was made up of 20% of the genotypes, including MP-3 and MP-4. MP-13 was assigned to a separate cluster due to its approximate average value for all phenotypic traits among the evaluated genotypes.

SLAF library validity and SLAF-seq

Since the genome of M. polymorpha was not published at the time of this study, we determined that the “RsaI + HaeIII” double-enzyme digestion scheme was the best enzyme digestion scheme for constructing the M. polymorpha SLAF tag library based on the prediction results of electronic digestion of the M. truncatula genome by Biomarker Technologies Corporation (Beijing, China). The selected size for this method was 314–414 bp using the Blue Pippin system. We predicted that there were 117,528 SLAF tags distributed in the genome with 6.19% located in repetitive regions. To further evaluate the efficiency of the restriction enzyme digestion scheme, Arabidopsis thaliana was selected as an example. Through the detection of the remaining restriction enzyme restriction sites in Arabidopsis read inserts, we found that the ratio of Arabidopsis paired-ended reads mapped to their reference genome was 83.41%, and the enzyme digestion efficiency was 97.79%. This indicates that the double-enzyme digestion scheme should be effective for the SLAF library construction of M. polymorpha.

SLAF sequencing of the 10 M. polymorpha accessions generated 6.26 Gb clean reads, with a Q30 percentage of 84.97% and GC content of 38.29% (Table S1). Finally, a total of 87,207 SLAF tags were identified, and a single accession contained 79,008 to 85,727 tags (mean 82,336 tags). The read numbers for SLAF tags ranged from 598,481 to 1,869,872 (mean 944,305), and the sequencing coverage ranged from 7.40-fold to 21.88-fold (mean 11.37-fold). Of the identified 87,207 SLAF markers, 25,616 (29.37%) were polymorphic and 70.53% were non-polymorphic, and only 0.1% were repetitive (Table S2). The integrity values of SLAFs from M. polymorpha accessions were all above 90% (Table S3). The high quality of these data promised the reliability of the subsequent analyses.

Development of SNP and SSR markers

Only polymorphic SLAF markers were further selected as high quality SNPs. PCR duplicates were not removed after alignment since the removal would result in many sites in the results not covered by reads, or very low coverage. Moreover, SLAF sequencing does not sequence the entire genome, but sequences the digested fragments according to the restriction sites. Therefore, data filtering strategies used the following parameters: QD >2, MQ >40.0, MQRankSum >−12.5, ReadPosRankSum >−8. After aligning the fragments to the M. polymorpha reference genome, a total of 52,237 SNPs with integrity > 90% and minor allele frequency (MAF) >0.05 were identified in this study (Table S4). The MAF distribution spectrum is listed in Table S5. These SNP variants were mostly distributed on pseudochromosome 3 (22.21%), least on pseudochromosome 7 (7%), and nearly evenly distributed on the five other pseudochromosomes (~14.16% for each pseudochromosome).

In order to provide a theoretical basis for the development of SSR molecular markers in M. polymorpha, the distribution characteristics of the SSR sequence were analyzed using MISA software. Results showed that 195,753 complete SSR sequences were screened, with relative densities of 428/Mb. The total length of SSR sequences was 3,611,698 bp and accounted for 0.79% of the total genome sequence length. Among the one to six different nucleotide repeats units, the single nucleotide repeat units were the most, followed by dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, and hexanucleotide repeat units. SSR repeat unit types in the M. polymorpha genome were abundant and had great potential for developing polymorphic markers and in the genetic diversity and molecular marker-assisted breeding of M. polymorpha.

Genetic diversity of M. polymorpha based on SNP markers and agronomic traits

Phylogenetic analysis based on high-quality SNPs also showed that all M. polymorpha accessions were divided into three clades (Fig. 1B), consistent with morphological clustering (Fig. 1A). The first clade contained MP-2, MP-5, MP-6, MP-7, MP-8, MP-9, and MP-11. The second clade included only MP-13, and the third clade had two accessions, MP-3 and MP-4. When M. truncatula was used as the outgroup, the rooted evolutionary trees (Fig. 2A) suggested that MP-3 and MP-4 were the most related to the putative ancestral species of the M. polymorpha population. Furthermore, PCA was performed and the first, second, and third principal components explained 14.88%, 12.75%, and 11.92% of the genetic diversity, respectively. The first principal component roughly separated the three groups (Fig. 2B). For the three clusters (Fig. 1A), the He, Nei, I, and PIC values of group I were the highest, followed by group III, and then group II (Table 4). These results indicated that group I possessed the highest genetic diversity, while group II had the lowest genetic diversity.