Development of polymorphic EST-SSR markers and characterization of the autotetraploid genome of sainfoin (Onobrychis viciifolia)

Tissue sampling and total RNA isolation

Sainfoin callus cells, emerging tidbits (<2 cm), young inflorescences (2–4 cm), inflorescences (4–6 cm), mature inflorescences (6–8 cm), developing seed pods (16 days after pollination (dap)), mature seed pods (24 dap), roots, germinated seeds (36 h after seed germination), young stems (less lignified), stems (moderately lignified), mature stems (highly lignified), young compound leaves, and mature compound leaves were harvested (Fig. 1). A total of 14 tissues in the Fig. 1, The callus cells were induced from inflorescences at 25 °C on solid MS medium containing 2,4-dichlorophenoxyacetic acid (3.0 mg/L) for 30 days under a 16/8 h (light/dark) light cycle (Ma et al., 2012). (A–G and J–N) were harvested from the same individual plant, these tissues used in this study were from the same 2-year-old plants grown inside a greenhouse set at 22 °C and a 16/8 h (light/dark) light cycle at Lanzhou University, radicle and germ (H and I) were harvested from three seedlings germinated from seeds. The sampled tissues were immediately frozen in liquid nitrogen and stored at −80 °C until use. Total RNA was isolated from 14 individually collected samples using the RNeasy Plant Mini Kit (Qiagen, Cat. #74904) as instructed. Concentrations of the RNA samples were determined using a NanoDrop ND1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA).

cDNA library construction and sequencing

To better elucidate tissue-specific RNA transcription, each RNA sample was adjusted to 400 ng/µL. Twenty micrograms of total RNA was taken from each RNA sample and pooled prior to cDNA library preparation using the mRNA-Seq Sample Preparation Kit (Illumina Inc., Beijing, China). Briefly, poly (A) mRNA was isolated from the pooled total RNA sample using magnetic oligo dT beads and used to synthesize first-strand cDNA with random hexamer primers and reverse transcriptase (Invitrogen, Beijing, China). Short cDNA fragments were purified using a MinElute PCR Purification Kit (Qiagen, Beijing China), resuspended in an EB buffer (Qiagen), and poly A was added. Sequencing adapters were ligated to the short cDNA fragments, and the libraries were sequenced using the Illumina HiSeq2500 sequencing platform at the BioMarker Company (Beijing, China). Processing of fluorescent images for sequence base-calling and calculation of quality values were performed using the Illumina data processing pipeline, which yielded 100-bp paired-end reads.

Sequence assembly and annotation

Before assembly, the raw reads were filtered to remove poly A/T, low-quality sequences, and empty reads or reads with more than 10% of bases having Q < 30. The assembly of a de novo transcriptome using clean reads was performed using the short-read assembling program Trinity. Contigs were generated after combining the reads with a certain degree of sequence overlap. Paired-end reads were used to detect contigs from the same transcript and the distances between contigs. Scaffolds were produced using N, representing different sequences between two contigs but connecting these two contigs together. Gaps between scaffolds were filled with paired-end reads and the reads with the lowest numbers of Ns. The resulting sequences were referred to as unigenes. The EST-SSR annotation positions of these unigenes were determined using BLASTX alignment (e-value < 10-5) against the sequences in the databases: NR, Swiss-Prot, KEGG, COG, and unigene sequences. After NR annotation, unigene GO annotations were conducted using the Blast2GO algorithm. GO functional classifications of the unigenes were performed using the WEGO software.

Detection of EST-SSR markers and designing of primers

EST-SSR markers were detected in the assembled unigenes using the Simple Sequence Repeat Identification Tool. The EST-SSRs were considered to contain mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides with a minimum of ten, six, five, five, five, and five repeats, respectively. EST-SSR primers were designed using the BatchPrimer3 software and synthesized by the Shanghai Sangon Biological Engineering Technology (Shanghai, China).

Sources of wild sainfoin populations

To produce results representing a wide range of sainfoin populations, we collected 40 individual wild sainfoin plants from five different locations (eight plants per location): Minqin, Jingyuan, Yuzhong, Huining and Maqu (Table S3). After air drying, leaves were taken from each plant and stored at room temperature until use.

DNA extraction

DNA was extracted from the harvested and dried leaves using the Nucleon Phytopure Genomic DNA extraction kit (Ezyp column plant genomic DNA extraction kit, Sangon Biotech Shanghai, China) on the samples described above using a protocol reported previously (Mora-Ortiz et al., 2016). The quality of each isolated genomic DNA was examined using 0.8% agarose gel electrophoresis, and the concentration of each genomic DNA sample was determined using a NanoDrop 8000 spectrophotometer (Nanodrop Technologies Wilmington, DE). The DNA samples were diluted individually in TE buffer to 25 ng DNA/µL prior to PCR amplification.

Amplification of EST-SSRs using polymerase chain reaction (PCR)

A total of 200 primer pairs were produced for this study and used to amplify EST-SSRs from the genomic DNA samples from the 40 wild sainfoin plants. PCR amplifications were performed in 5-µL reactions (0.5 µL of DNA, 2.5 uL of 2 × mix (e.g., 0.5 µL of 2 × PCR buffer, 1 µL of 1 mM dNTPs, 0.4 µL of 25 mM MgCl₂, 0.1 µL of Taq DNA polymerase), 0.5 µL (5 pmoL/µL) of forward and reverse primers, and 1uL of sterile distilled water). Three microliters was added to the PCR reaction in each tube, and PCR was performed using a PCR-100TM Thermal cycler set at 94 °C for 4 min, followed by 35 cycles of 94 °C for 1 min at a specific annealing temperature for 30 s and 72 °C for 20 s. The final extension was at 72 °C for 7 min. The resulting PCR products were resolved using 8.0% non-denaturing polyacrylamide gels (Lot# I20826, GelStain; Beijing TransGen Biotech. Co., Ltd., Beijing, China) after electrophoresis.

Construction of T vectors

Fragments of 200–500 bp were PCR amplified from sainfoin genomic DNA using the gene-specific primers. The PCR products were purified and these fragments were then cloned into pGEM-T Easy Vector (Promega, Madison, WI, USA). The constructs were transferred into E.coli DH5 α using the freezing/heat-shock method, and sequenced by Shanghai Sangon Biological Engineering Technology using Sanger dideoxy (Shanghai, China).

Diversity analysis

The observed heterozygosity (Ho) was calculated as previously shown (Liu, Liu & Yang, 2007), and the corrected heterozygosity (He), corrected for sample size, and the average polymorphism information content (PIC) were analyzed using the ATETRA 1.2.a software program. Only specific bands that could be unambiguously scored across all individual plants were used in this study. A clustering analysis was used to generate a dendrogram using the unweighted pair-group method with arithmetic mean (UPGMA) and Nei’s unbiased genetic distance with the NTSYSPC 2.0 software package. We used STRUCTURE 2.3.4 to generate a structure map.

RNA sequencing and de novo assembly

A cDNA library representing 14 different sainfoin tissues (Fig. 1) was sequenced, and a total of 26,912,927 raw reads were obtained (Table 1). After rigorous quality checks and data filtering, a total of 24,630,711 high-quality clean reads were obtained, and about 92% of them had quality (Q)-scores above Q30. These clean reads together contained a total of 6,264,706,761 nucleotides (nt), without N, and about 45% GC content. The high-quality reads were deposited in the U.S. National Center for Biotechnology Information (NCBI) sequence read archive (SRA) database (SRX3763386). From these high-quality reads, a total of 2,678,687 contigs with a mean length of 72.26 bp and an N50 length of 69 bp were generated using de novo assembly. The total number of unigenes from the paired-end reads was 77,764, and 27,437 of them had distinct clusters. A total of 50,327 unigenes had distinct singletons, and the total length of these unigenes was 53,035,704 bp. The average length of a unigene was 682.01 bp with an N50 value of 1,209 bp. Of the 77,764 unigenes, the length of 50,327 unigenes ranged from 200 to 500 bp, the length of 22,096 unigenes ranged from 500 to 2,000 bp, and the length of 5,341 unigenes was above 2,000 bp (Table 1). Also, 36,353 of the 77,764 unigenes were successfully annotated according to the NR, Pfam, Swiss-Prot, KEGG, COG, and GO databases (Table 2), and 10,387 unigenes were assigned to COG classifications.

Frequency and distribution of EST-SSRs

A total of 6,752 potential EST-SSRs were identified in the 77,764 unigenes (Table 1) and used to design 2,469 primer pairs. Since 1,271 unigenes contained more than one EST-SSR, the types and distributions of the total 6,752 potential EST-SSRs were analyzed. As shown in Fig. S4, the density of single nucleotide repeats was the highest (88 of SSRs per Mb), followed by trinucleotide repeats (64 of SSRs per Mb). The most abundant repeat type was the mono-nucleotide repeat (2,906 repeats or 43.04% of the total repeats), followed by tri-nucleotide repeats (2,262, 33.50%), di-nucleotide repeats (1,287, 19.06%), quad-nucleotide repeats (263, 3.90%), hexa-nucleotide repeats (19, 0.28%), and penta-nucleotide repeats (15, 0.22%) (Table 3). Moreover, the EST-SSRs with five tandem repeats were the most common EST-SSRs (24.88%), followed by ten tandem repeats (22.07%), six tandem repeats (16.88%), eleven tandem repeats (10.19%), seven tandem repeats (7.09%), and the remaining tandem repeats (<5%) (Table 3). In addition, the GO enrichment for the 6,752 SSR-containing unigenes was done using the agriGO algorithm (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) and the 971,445 unigenes available in the database as the reference. The results of the GO enrichment analysis indicated that the proportion of the “transcription”-related (GO:0003674) unigenes was significantly enriched (Fig. S5).

Development of EST-SSR markers

Using the SSR-containing unigene sequences, 200 primer pairs were randomly chosen from the 2,469 identified primer pairs, synthesized, and used to determine if the EST-SSR loci identified in this study were true-to-type EST-SSR loci in the sainfoin populations. Of the 200 primer pairs, 178 of them successfully amplified fragments from sainfoin genomic DNA during PCR, while the other 22 primer pairs failed. Also, 132 of the 178 PCR primer pairs amplified products of the expected size. Using genomic DNAs from the 40 different wild sainfoin plants (Table S3) as templates, 61 of the 132 primer pairs were found to be polymorphic (Fig. S6 and Table S4) and the other 71 primer pairs were monomorphic (Table S2).

Assessment of sainfoin genetic diversity

The 61 primer pairs mentioned above were used to analyze the genetic diversity among the population comprising 40 wild sainfoin plants from five different geographic locations. The result showed that a total of 459 alleles were present in the 61 polymorphic loci in the 40 different individuals, and the number of alleles per loci ranged from three to twelve with an average number of 7.52. The Ho, He, and PIC were estimated from 0.05 to 1.0 (mean value = 0.67), 0.55 to 0.91 (mean value = 0.77), and 0.51 to 0.88 (mean value = 0.74), respectively (Table S4). These 61 polymorphic loci or EST-SSR markers are unlinked and have high degrees of universality among the assayed germplasms. Therefore, they are useful for studying biogeographic processes that shaped the current disjunctive distributions of sainfoin. Furthermore, PCR amplicons representing EST-SSRs from different sainfoin individuals were sequenced and the results showed that all of the sequenced alleles were homologous to the original locus from which the marker was designed.

Using NTSYSPC 2.0 software and the UPGMA method, a dendrogram was obtained (Fig. 2). The 40 sainfoin individuals can be divided into five distinct wild populations. The individuals in population 1 originated from Minqin, the individuals in population 2 were from Jingyuan, the individuals in population 3 were from Yuzhong, the individuals in population 4 were from Huining, and the individuals in population 5 were from Maqu (Table S3). The individual plants from the Jingyuan sampling site were assessed on related taxa within the sainfoin genus. Most of them were easily amplifiable and detectable across all genotypes, and just a few showed problems of amplification or scoring, most likely due to polymorphisms (insertions/deletions or base mutations) in primer regions or in regard to subspecies. However, most EST-SSR markers found in this study are highly useful for discriminating related sainfoin populations or related taxa, including wild species, subspecies, and subgenera. Consequently, we propose that these EST-SSR markers can be used in future studies on comparative genomics, genetic differentiation, and evolutionary dynamics within the sainfoin genus.

STRUCTURE 2.3.4 was used to generate a structure map. STRUCTURE analysis based on 459 loci representing EST-SSRs was performed to evaluate the genetic structure of the 40 wild sainfoin individuals. The highest ΔK was observed for K = 6 [ΔK(6) = 144]. ΔK values for K = 3–5 and K = 7 were not significant (ΔK = 0.078–5.632). The mean value of the log probability of the data was higher for K = 6 than for K = 4, and K = 5 [LnP (D) K = 6 =  − 13190.22, LnP (D) K = 4 =  − 14038.11, LnP (D) K = 5 =  − 13599] (Table S6). Therefore, six clusters were chosen as the most probable genetic structure of the wild sainfoin individuals. With K = 6, seven individuals from site 1 were assigned to cluster 1 with coefficient Q values ranging from 0.782 to 1.000; one individual from site 2 to cluster 2 with a Q value between 0.900 and 1.000; one from site 2 to cluster 3 with a Q value between 0.900 and 1.000; three from site 2 and three from site 3 to cluster 4 with Q values from 0.782 and 0.982; seven from site 4 to cluster 5 with Q values from 0.683 to 0.973; and eight from site 5 to cluster 6 with Q values from 0.714 and 0.965 (Fig. 3). Ten sainfoin individuals could not be assigned to any of the clusters due to high levels of admixture (Q < 0.6).

Furthermore, PCR amplicons of Vo61 and Vo157 EST-SSRs from single individuals were sequenced to check the authenticity of the SSR locus (Fig. S6 and Fig. 4). For individual 2 in accession II at the Vo61 locus, four amplicons with different AGAA repeats, (AGAA)₃, (AGAA)₄, (AGAA)₅, (AGAA)₆, were found by pGEM-T easy vector sequencing, while three amplicons with different ACC repeats, (ACC)₃, (ACC)₅, (ACC)₆, were found for individual 1 in accession II for the Vo157 locus.