Development of Cymbidium ensifoliumgenic-SSR markers and their utility in genetic diversity and population structure analysis in cymbidiums

Cymbidium is a genus of 68 species in the orchid family, with extremely high ornamental value. Marker-assisted selection has proven to be an effective strategy in accelerating plant breeding for many plant species. Analysis of cymbidiums genetic background by molecular markers can be of great value in assisting parental selection and breeding strategy design, however, in plants such as cymbidiums limited genomic resources exist. In order to obtain efficient markers, we deep sequenced the C. ensifolium transcriptome to identify simple sequence repeats derived from gene regions (genic-SSR). The 7,936 genic-SSR markers were identified. A total of 80 genic-SSRs were selected, and primers were designed according to their flanking sequences. Of the 80 genic-SSR primer sets, 62 were amplified in C. ensifolium successfully, and 55 showed polymorphism when cross-tested among 9 Cymbidium species comprising 59 accessions. Unigenes containing the 62 genic-SSRs were searched against Non-redundant (Nr), Gene Ontology database (GO), eukaryotic orthologous groups (KOGs) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The search resulted in 53 matching Nr sequences, of which 39 had GO terms, 18 were assigned to KOGs, and 15 were annotated with KEGG. Genetic diversity and population structure were analyzed based on 55 polymorphic genic-SSR data among 59 accessions. The genetic distance averaged 0.3911, ranging from 0.016 to 0.618. The polymorphic index content (PIC) of 55 polymorphic markers averaged 0.407, ranging from 0.033 to 0.863. A model-based clustering analysis revealed that five genetic groups existed in the collection. Accessions from the same species were typically grouped together; however, C. goeringii accessions did not always form a separate cluster, suggesting that C. goeringii accessions were polyphyletic. The genic-SSR identified in this study constitute a set of markers that can be applied across multiple Cymbidium species and used for the evaluation of genetic relationships as well as qualitative and quantitative trait mapping studies. Genic-SSR’s coupled with the functional annotations provided by the unigenes will aid in mapping candidate genes of specific function.


Background
Cymbidium is a genus of 68 species in the orchid family [1]. Cymbidium species are mainly distributed in the tropical and subtropical regions of Asia, including northwest India, China, Japan, Korea, the Malay Archipelago, and north and east Australia [2,3]. A total of 49 species can be found in China, including five famous species, i.e., C. goeringii, C. faberi, C. ensifolium, C. kanran, and C. sinense. These cymbidiums comprise some of the rarest plant species, with only a few surviving original populations and some reintroduced plants in the south of China, including Yunnan and Taiwan [4]. The fascinating varieties and shapes of their flowers endow these species with extremely high ornamental value that has attracted the world? s attention.
Knowledge of the genetic diversity and population structure of germplasm collections is an important foundation for plant improvement [5]. Estimation of genetic distance among germplasm is helpful in selecting parental combinations for creating segregating populations so as to maintain genetic diversity in a breeding program. However, genetic diversity may appear spatially structured at different scales, such as population, subpopulation or among neighboring individuals [6]. Population genetic analyses can provide important parameters including standing levels of genetic variation and the partitioning of this variability within/ between populations [7]. The genetic diversity or population structure of C. ensifolium and other cymbidiums have been measured by using different molecular tools, including restriction enzyme polymorphism (RFLP) markers [3], random amplified polymorphic DNA (RAPD) markers [3,4,8], amplified fragment length polymorphism (AFLP) markers [4], polymorphisms of internal transcribed spacers (ITS) of nuclear ribosomal DNA and plastid, inter-simple sequence repeats (ISSR) markers [4,9], and SSRs [10,11]. Compared with RAPD, ISSR and ITS, SSR markers are more reliable, locus-specific, codominant, highly polymorphic, and well distributed throughout the genome [12]. Moreover, SSR? s only require polymerase chain reaction (PCR), which is a big advantage over RFLP and AFLP. These features make SSR? s well suited for marker-assisted selection, genetic diversity analysis, population genetic analysis, genetic mapping, and genetic map comparison in various species [13,14].
The number of SSR is very limited for C. ensifolium, due to limited sequence resources. Until now, the National Center for Biotechnology Information (NCBI) contained very limited Cymbidium sequence information, i.e., 692 nucleotide sequences and 78 expressed sequence tags (ESTs) (http://www.ncbi.nlm.nih.gov/nucest?term=cymbi-dium%5BOrganism%5D, verified 2014). RNA-seq provides a fast, cost-effective, and reliable approach for generating large-scale transcriptome data in non-model species, and also offers an opportunity to identify and develop genic-SSRs by transcriptome data mining [15]. Compared with traditional ? anonymous? SSRs from genomic DNA, these new genic-SSR markers have two advantages, i.e. a wealth of functional annotations and high transferability across taxa [15,16]. Herein, we extracted the total mRNA from C. ensifolium flower buds for RNA-seq, which resulted in 9.52 Gb of transcriptome data. From the C.ensifloium transcriptome, we obtained 55 new polymorphic microsatellite loci after testing their transferability across 59 Cymbidium accessions.

Plant materials
A total of 11 C. ensifolium accessions were employed to test genic-SSRs and additional 47 accessions from C. lancifolium, C. floribundum, C. suavissimum, C. cyperifolium, C. qiubeiense, C. faberi, C. goeringii and C. sinense were used to cross-test these markers among multiple species. The plants were grown and maintained in a greenhouse at the Zhejiang University under natural light (Table 1). Fresh leaf samples were collected from two or three seedling of each accession for genomic DNA extraction.

Genic-SSR search and primer design
Total RNA was isolated from native cultivar of C. ensifolium ? Tiegusu? using TRIzol? reagent (Invitrogen, CA, USA) and treated with RNase-free DNase I (TaKaRa Bio, Dalian, China) for 45 min according to the manufacturer? s protocol. The RNA was used in cDNA library construction and Illumina deep sequencing [17]. The raw sequencing reads were stringently filtered, and high-quality reads were assembled de novo using Trinity with an optimized k-mer length of 25 [18]. MSATCOMMANDER V. 0.8.2 [19] was used to analyze SSR distribution. The minimum number of repeats for SSR detection was as follows: six for di-SSRs, and four for tri-, tetra-, penta-, and hexa-SSRs. The open reading frame (ORF) and untranslated region (UTR) within unigenes were identified using Trinity [18]. Software Primer3.0 [20] was used to design primers for genic-SSR loci with sufficient flanking sequences.
Unigenes containing genic-SSRs were compared with protein databases, including the non-redundant (Nr) database (http://www.ncbi.nlm.nih.gov/), using BLASTX with a significance cut-off E-value of 1e −5 [17]. For the non-redundant annotations, BLAST2GO V. 2.4.4 was used to obtain Gene Ontology (GO) annotations of unique transcripts [21]. Metabolic pathway analysis were performed based on the pathways of Oryza sativa in the Kyoto Encyclopedia of Genes and Genomes (KEGG) [22,23]. The unigene sequences were also aligned to the KOG (Eukaryotic Orthologous Groups) database to predict and classify possible functions [24].

Genotyping
Genomic DNA was extracted from leaf samples as previously described [25]. PCR primers were synthesized by Life Technologies (AB & Invitrogen, Shanghai, China). PCR reactions were conducted based on a previously published protocol [26]. The PCR products were separated through polyacrylamide gel electrophoresis using 8% bis-acrylamide, 0.5% TBE buffer, 0.07% APS, and 0.035% TEMED. The gel was run at constant 120 V for approximately 3 h in 1? TBE buffer. The gel was silverstained according to Li? s procedure [27], and was then documented using a scanner. The genotype was determined by analysis of the bands? pattern, dependent on the number and the position of bands.

Statistical analysis
Genetic distance was calculated using Nei? s distance [28]. Phylogenetic reconstruction was based on the unweighted pair-group method that utilizes the arithmetic average (UPGMA) method implemented in PowerMarker version 2.7 [29]. The tree that was used to visualize the phylogenetic distribution of accessions and ancestry groups was constructed using MEGA version 4 [30]. A model-based program structure [31] was used to infer population structure with 5,000 burn-in and run length. The model allowed for admixture and correlated allele frequencies.
The number of groups (K) was set from 1 to 10, each with 10 independent runs. The most probable structure number (K) was determined through log probability [32]. Principal component analysis (PCA), which summarizes the major patterns of variation in a multi-locus data set, was performed using NTSYSpc version 2.11 V [33]. Two principal components were used to represent the dispersion of the collection accessions graphically [34]. PowerMarker was used to calculate the average number of marker alleles and the polymorphism information content (PIC) values. Fixation index (Fst), which indicates the differentiation among genetic groups, was calculated using an Analysis of Molecular Variance (AMOVA) approach in Arlequin V2.000 [35].

Genic-SSR search and primer design
In C. ensifolium transcriptome, 98,819,349 reads, (9.52 Gb), were obtained after removal of adaptor sequences, ambiguous reads, and low-quality reads (Q-value <25). These reads were used for the subsequent assembly, and then resulted in 101,423 unigenes (139,385,689 residues). The length of unigenes averaged 1,374 bp and ranged from 351 bp to 17,260 bp. The data were uploaded to the NCBI (http:// orchidbase.itps.ncku.edu.tw/est/home2012.aspx) for public use (Accession: SRA098864).
In the present study, 7,936 genic-SSRs were identified, with one SSR locus for every 17.56 kb (kb/SSR). Estimated locations (coding, 5′UTR or 3′UTR) were obtained for 5,524 genic-SSRs. Sequence information could not be determined for the remaining 2,412 genic-SSR regions, because the locations were extended over both estimated coding and non-coding regions. Given such high numbers of SSR, we analyzed the sequence data to isolate high quality SSR loci for further testing. An important factor considered was the locations of SSRs relative to ORFs. SSRs within UTR are exposed to lower selective pressure than those in coding regions and have a higher likelihood of being polymorphic [36]. Another two factors are the length of the motif and the number of the repeat motif, which are often associated with polymorphism [37]. Thus, SSR? s within UTR, with short motifs and high repeat number would be the best marker candidates. Herein, we selected 80 genic-SSRs and designed primers based on their motifs, sizes and locations.

Genetic diversity and population structure
These genic-SSRs revealed genetic variation among accessions. The genetic distance among accessions ranged from 0.016 to 0.618, with an average of 0.391. The model-based clustering method revealed five groups ( Figure 1A and B). Group 2 had the most accessions (26), with the highest mean genetic distance (MGD) of 0.431 among these accessions; Group 4 had 10, with an average distance of 0.236; Group 5 had 7, with MGD of 0.332; Group 1 and Group5 both had 7 accessions, with MGD of 0.155 and 0.332, respectively; Group 3 had 9, with MGD of 0.213. Genetic distance among five groups was from 0.340 (between group 1 and group 5) to 0.176 (between group 2 and group 4, with average of 0.248) ( Table 3).
The five groups revealed by the model-based clustering analysis consisted of different species. Three groups comprised more than one species, whereas the other two only comprised one species. Group 1 included two species i.e. C. cyperifolium and C. goeringii; Group 2 included C.ensifolium, C. lancifolium, C. suavissimum, C. qiubeiense, C. goeringii, C. faberi, and C. sinense; Group 5 included C. floribundum, C. suavissimum and C. goeringii. Goup 3 and Group 4 included only C. faberi and C.ensifolium, respectively ( Figure 2).
The first two components in PCA (47.87% and 21.59% of total variation, respectively) discriminated the five groups at a certain level. Basically, accessions in group 1 and group 3 stayed alone, whereas group 2 overlapped with group 4 and group 5 ( Figure 1C). In the phylogenetic tree, group 2 and group 4 were genetically close, while group 5 was relatively distant from the other groups ( Figure 1A). In addition, a few accessions in group 2 had admixture ancestry from group 3 and group 4, while accessions in group 3 and group 1 had less admixture ancestry ( Figure 1B). AMOVA results showed that 25.34% of the total variation was among groups, while 74.66% of the variation was within groups. The F ST was 0.25, as indicated by the AMOVA approach.

Genic-SSR annotation
Annotations of these unigenes provide biological information for 62 genic-SSRs, such as KOG clusters, GO, and KEGG pathway information. Distinct gene sequences were first searched using BLASTX against the Nr database. The results showed that 53 unigenes had hits that exceeded the E-value threshold. In the present study, 39 unigenes were categorized into 25 GO terms in three GO ontologies ( Figure 3A). Two groups ? membrane? and ? nucleus? , one group ? binding? , and one group ? cellular process? comprised the most representative genes found in cellular components, molecular function, and biological processes, respectively. Out of 53 hits in the Nr databases, 18 sequences were classified into 9 KOG categories ( Figure 3B). Among the 9 KOG categories, ? General function prediction only? and ? Posttranslational modification, protein turnover, chaperones? were the two largest groups. When referenced to rice (Oryza sativa), 15 unigenes were found to be involved in 14 pathways ( Figure 3C). The most highly representative one was ? metabolic pathways? , where unigenes shared similarity with 18 rice sequences.    Note: A total of 62 genic-SSR markers successfully amplified were listed, however 55 polymorphic markers were used in subsequent population analysis or cross species comparison. NULL: no annotation. NA: monomorphic marker.

Diversity
Because genic-SSR markers are derived from transcribed regions of DNA, they are expected to be more conserved and have a higher rate of transferability than anonymous SSR markers [38]. Herein, 55 C. ensifolium polymorphic genic-SSR markers exhibited 100% transferability across the 59 accessions of the 9 Cymbidium species tested. It is common that genic-SSRs possess a high potential for inter-specific transferability [39,40]. Other markers such as RAPD? s, ISSR? s and non-genic SSR? s have also been used with success among C. ensifolium and the Cymbidium species reflecting the genetic similarity among many members of the genus [8,11,15]. The conserved nature of the genic-SSRs may limit their polymorphism relative to randomly selected SSR? s.   In this study, PIC of genic-SSR markers averaged 0.407, lower than 0.782 [5] and 0.639 [11] of anonymous SSR? s tested on Chinese cymbidiums in other studies. The pairwise genetic distance averaged 0.391 among 59 accessions, which is also lower than that from previous studies conducted on Chinese Cymbidiums using other molecular markers [3,8,[41][42][43][44]. Even though genic-SSRs revealed less variability than SSRs, these markers still reveal sufficient levels of variation for population genetic analysis.

Population structure
One of the biggest advantages for genic-SSRs is that they allow one to make direct comparisons among taxa without running the risk that locus-specific differences might mask true species-level differences, such as overall levels of genetic diversity, the extent of population structure, and so on. However, the greatest concern with the utilization of genic-SSRs in genetic studies is that selection on these loci might influence the estimation of population genetic parameters. While a recent study by Woodhead et al. [45] revealed that estimates of population differentiation based on genic-SSRs are comparable to those based on both SSRs and AFLPs in ferns, and large-scale comparative analysis suggest that only a very small percentage of all genes has experienced positive selection [46,47], a small fraction of SSRs will be inevitably subject to selection. The view is consistent with the theory that most mutations are neutral, or nearly neutral, [48] or, at least, do not change the function of gene products appreciably [49].
In the population genetic analysis, almost all accessions from the same species clustered together. C. suavissimum and C. floribundum were clustered into one brand, and clearly distinguished from other cymbidiums. Two of them belong to Section Floribundum, and have a distant relationship with other cymbidiums. However, the genetic relationship between C. goeringii and C. sinensis was close, which was congruent with the previous reports [5,11]. The close relationship was also found between C. ensifolium and C. cyperifolium. In the intersection level, we discovered that two accessions of C. faberi were clustered with C. cyperifoliumm, and accessions of C. lancifolium and C. ensifolium were scattered among ones of C. goeringii. The splitting feature of these clusters might be linked to the non-homologous synapomorphy, even though accessions belonged to different species. The accessions of C. goeringii did not always form a separate cluster in the phylogenetic tree or were not grouped together in structure analysis, suggesting that they were polyphyletic. Previous morphologic, cytogenetic, and molecular studies have shown that the major lineages of Chinese cymbidiums are ambiguous. C. ensifolium and C. sinense are classified in section Jensoa; C. faberi and C. goeringii, are classified in section Maxillarianthe; C. faberi, C. kanran, and C. longibracteatum are classified in one group; C. ensifolium, C. goeringii, and C. sinense are categorized into another group [44].

Genic-SSR annotation
Putative functions were assigned to those unigenes containing SSRs by sequence similarities. These unigenes were involved in a wide range of functions, which indicated that these genic-SSRs were likely important biologically characters. For example, unigene containing SSR47 shares homology with CONSTANS-like protein. In Arabidopsis, the CO (CONSTANS) gene has an important role in the regulation of flowering by photoperiod [50]. Unigene containing SSR43 has homology with a glycinebetaine/proline transporter. The accumulation of glycinebetaine (GB) is one of the adaptive strategies to adverse salt stress conditions [51]. The transporters mediate the uptake of GB and/or proline in many plant species e.g. Arabidopsis thaliana [52], tomato (Solanum lycopersicum) [53], rice (Orazy sativa) [54], barley [55]. Unigene having SSR75, was annotated as mitogen-activated protein kinase (MAPK). MAPK cascades function as key signal transducers that use protein phosphorylation/dephosphorylation cycles to channel information [56]. In the plant, MAPKs have been shown to regulate numerous cellular processes, including biotic stress relief [57,58]. Although some unigenes with SSRs had no match to known genes in current gene database, they will likely gain functional annotations as the knowledge of plant genes increases. Compared with anonymous SSRs, genic-SSR markers have a higher probability of being functionally associated with differences in gene expression, which may be in identifying associations between genotype and phenotype. Mapping of genic-SSRs will also provide a map location, in many cases, for genes with known functions.

Conclusion
In this work, 7,936 genic-SSRs were identified in C. ensifolium transcriptome and their characterizations were further analyzed. A total of 80 genic-SSRs were chosen for validation, and 55 markers successfully yielded polymorphism across 9 Cymbidium species including 59 accessions. The high transferability of genic-SSR will be a powerful resource for molecular taxonomic studies and construction of a reference molecular map of the Cymbidium genome. Since genic-SSR markers belong to generich regions of the genome, some of these can be exploited for use in marker-assisted breeding of Cymbidium. Therefore, the set of genic-SSR markers developed here is a promising genomic resource.