Development of a EST dataset and characterization of EST-SSRs in a traditional Chinese medicinal plant, Epimedium sagittatum (Sieb. Et Zucc.) Maxim

Contig assembly of 454 reads

A total of 228,768 sequences (51.7 Mb) are generated from leaf cDNAsof E. sagittatum. After removing SMART primers, polyA tail, and low quality sequences, a total of 217,380 high quality sequences is retained with a total length of 50.9 Mb and an average length of 224.9 bp ranging from 50 bp to 556 bp (Additional file 1). After clustering and assembly, 76,459 consensus sequences with a total length of 18.8 Mb (6.5 Mb of 17,231 contigs and 12.3 Mb of 59,228 singlets) are obtained (Additional file 1). The length of contigs ranges from 51 bp to 1,911 bp with an average of 375.9 bp, and that of singlets ranges from 50 bp to 556 bp with an average of 207.8 bp (Additional file 1).

Structural and functional annotation of ESTs

The consensus sequences are blasted against non-redundant protein database (NR database) and 29,466 consensuses comprising of 10,245 contigs and 19,221 singlets have significant hits to NR database (Additional file 1). All the consensus sequences are also compared against Pfam domain database and about 50% (38,220) of them are successfully annotated. Gene Ontology (GO) terms are further assigned to 38,220 Epimedium consensus sequences that are successfully annotated based on Pfam matches. Only 22,295 Epimedium consensus sequences are able to be assigned GO numbers. Figure 1 shows the percentage distributions of gene ontology terms according to the GO consortium. Protein metabolic process (23.4%) is the most dominant group out of 12,242 consensus sequences that are annotated to the biological process category. It is followed by the biosynthetic process at 23.2%, photosynthesis at 22.5%, generation of precursor metabolites and energy at 19.3%, cellular process at 15.7%, translation at 13.6%, metabolic process at 10.5%, and carbohydrate metabolic process at 9.3% (Figure 1A). A total of 13,807 sequences could be assigned to the molecular function category. Binding (30.5%) is the most dominant group followed by catalytic activity (20.6%), structural molecular activity (12.7%), hydrolase activity (11.9%), nucleotide binding (10.5%), transferase activity (6.9%), and transcription regulator activity (2.2%) (Figure 1B). With regard to the cellular component, 64.9% of 12,308 consensus sequences are assigned to membrane followed by intracellular (15.1%), ribosome (12.0%) and extracellular region (5.8%) (Figure 1C).

Being famous medicinal plants used in China, Epimedium species have flavonoids as the major phytochemical components. There are also many other unidentified compounds in Epimedium species [29]. Thus, genes related to secondary metabolism, especially for the flavonoid biosynthesis, are under further studies in order to decipher the molecular mechanism of natural variation of phytochemical components and facilitate the molecular breeding of medicinal cultivars. In our EST collection, 29 EST consensus sequences are annotated relating to secondary metabolic process (Figure 1A), including genes encoding key enzymes in the flavonoid biosynthetic pathway such as phenylalanine ammonia-lyase, cinnamate-4-hydroxylase, 4-counaroyl:ligase CoA, chalcone synthase, chalcone isomerase, flavanone 3 hydroxylase, flavanone 3' hydroxylase, flavanone 3'5' hydroxylase, flavonol synthase, flavone synthase, dihydroflavonol 4-reductase, anthocyanidin synthase, rhamnosyl transferase, and UDPG-flavonoid glucosyl transferase. In addition, more than 300 EST consensus sequences were annotated relating to transcription regulator activity (Figure 1B), part of which might play roles in regulating flavonoid biosynthesis and metabolism. For example, regulatory genes are annotated as members of WD40 family, MYB family, and bHLH family that can form a triplex compound to regulate the specific expression patterns of flavonoid structural genes. All these structural and regulatory genes are being cloned and functional studies are carrying out in our lab. Further studies such as molecular modification (or gene transformation) and metabolic engineering of enzymes are also in progress.

Frequency and distribution of EST-SSRs in the Epimedium transcriptome

The Epimedium EST-SSR Transferability

A total of 32 primer pairs (Additional file 2), referred to EsESP01 to EsESP32, have been stochastically selected among the 2,810 microsatellites to amplify 10 dinucleotide SSRs and 22 trinucleotide SSRs for testing the transferability in other Epimedium species (Additional file 3). Of these primer pairs, seven primer pairs (EsESP1, EsESP3, EsESP4, EsESP20, EsESP24, EsESP31, and EsESP32) cannot amplify any fragment in E. sagittatum, suggesting that these primers are not well designed. Four primers, EsESP14, EsESP19, EsESP26, and EsESP27, amplify numerous non-target bands, suggesting these primers are also problematic. Three primer pairs, EsESP05, EsESP25, and EsESP29, succeed in transferring to other Epimedium species, but amplify PCR products larger than the expected sizes. The rest of primer pairs produces polymorphic bands in most Epimedium species as shown in figure 2. Therefore, 85.7% (18 out of 21) EST-SSR primers can be transferred successfully to other Epimedium species. In theoretically, there are 56.2% (18/32) of 2810, namely ~1580, E. sagittatum EST-SSRs markers suitable for genetic studies in genus Epimedium.

Species from other genera, including Vancouveria hexandra, Nandina domestica Thunb., Berberis julianae Schneid., Berberis sargentiana Schneid., Dysosma versipellis (Hance) M. Cheng ex Ying, Mahonia bealei (Fort.) Carr., and Mahonia fortunei (Lindl.) Fedde., have been used to test the transferability of Epimedium EST-SSRs across different genus of Berberidaceae. A total of 8 primers, EsESP6, EsESP7, EsESP9, EsESP10, EsESP11, EsESP13, EsESP23, and EsESP28, can produce PCR fragments in V. hexandra (the most closely related genus reported by Wang et al [30]) with similar sizes and same SSR motifs as those in Epimedium (sequence data not shown). Only four primers, EsESP8, EsESP13, EsESP18, and EsESP21, can amplify PCR fragments in other Berberidaceae species mentioned above. However, sequences of the PCR fragments from genera Berberis, Nandina, Dysosma, and Mahonia, do not contain the SSR motifs found in Epimedium species (data not shown). Therefore, Epimedium EST-SSRs can successfully transfer to Vancouveria, but not other Berberidaceae genera tested.

The genetic diversity and dendrogram in genus Epimedium

A total of 55 individuals representing 52 species (Additional file 3) have been used to evaluate the genetic diversity in genus Epimedium. Two primer pairs, EsESP07 and EsESP08, detect a single band in some species while null in other species. Thus, they are excluded for further study and the remaining 16 EST-SSRs are used to estimate the genetic diversity among 52 species. A high level of genetic diversity in genus Epimedium is discovered (Additional file 2). The number of alleles per locus ranges from 3 to 27 with an average of 11.9 alleles. Ho ranges from 0.04 to 0.6 with an average of 0.35. He ranges from 0.17 to 0.94 with an average of 0.65. Polymorphism information content (PIC) is from 0.17 for EsESP22 locus to 0.93 for EsESP21 locus. The observed length of PCR products vary greatly in a large spectrum compared with the fragment size in E. sagittatum (Additional file 2). The size variation of amplicons reflects the length of repetitive sequences in each SSR locus, suggesting that the polymorphism might be related to the differences of tandem repeat numbers (Figure 2).

In order to study the evolutionary relationship among Epimedium species, those 16 EST-SSR primers successfully transferring in almost all Epimedium species are used for reconstructing Epimedium phylogeny. V. hexandra is used as an outgroup species, because it is more closely related to Epimedium genus than Jeffersonia dubia (Maxim.) Benth. et Hook. f. and Podophyllum peltatum [30]. Furthermore, sequences of some loci in V. hexandra amplified by Epimedium EST-SSR primers confirm that similar single nucleotide repeats exist in V. hexandra as mentioned above. As shown in additional file 4, the phylogenetic relationship among all 52 species are supported by low bootstrap values, especially for Chinese species, partially due to the small number of EST-SSRs and the nature of EST-SSR with low PIC used in this study. Previous studies also found it is hard to differentiate the closely related Chinese species using different types of molecular markers including RAPD [12], RFLP [12], AFLP [13], and rDNA [6, 7, 14]. Therefore, more EST-SSR markers or combining results of different types of markers would provide more opportunities for better resolving the phylogeny among Epimedium species, which is an on-going project in our lab.

E. sagittatum ESTs annotation

Currently, numerous studies are focusing on isolating the Epimedium bioactive components [3–5], but the potential molecular mechanism producing bioactive components is still unclear. Furthermore, little or no information regarding to its biological conservation, population genetics, molecular breeding, and molecular biology is reported for Epimedium species due to the lack of background information. The 454 GS-FLX pyrosequencing technology is a very cost-effective method to obtain large-scale EST sequences and will speed up the researches of less-studied Epimedium species. The technology has been broadly used for bacteria de novo sequencing, cDNA sequencing, small RNA sequencing, metagenomic sequencing and whole genome shotgun survey http://www.454.com/news-events/publications.asp. Margulies et al. [31] demonstrated that the technology had 100 fold sequencing capability at the cost of relative low accuracy than Sanger-based capillary electrophoresis sequencing systems when the GS20 sequencing machine was used. Recently, the new generation machine, GS-FLX, possesses more powerful sequencing capability of up to 200-300 bp per read and higher sequencing accuracy than that of GS20 http://www.454.com. In this study, an Epimedium transcriptome population is sequenced by 454 GS-FLX system, and a total of 226,544 sequence reads are generated. After assembly, 76,459 consensus sequences are obtained.

It is noticeable that only 29,466 Epimedium consensuses have significant hits in the NR protein database (Additional file 1). Thus, there are a large number of sequences that don't have any significant match in other species (46,993). This might be due to the following reasons: (1) the shortness of sequence reads leading to assemble difficultly, (2) inefficiency of blast using short sequences and the incompleteness of the known database, (3) the advantages of deep sequencing, which can discover novel genes with low expression levels, (4) the shortness of sequence reads resulting in low efficiency of annotation, which was confirmed by correlation test between sequence length and percentage of sequences successfully annotated (Additional file 5), and (5) part of consensuses representing UTRs. With regard to the first point, the size of raw reads ranged from 50 to 556 bp with an average length of 224.9 bp. After assembly, most contigs still range from 201 to 300 bp, and the average length only increased up to 245.7 bp. It is also supported by the fact that most contigs were assembled with 2-5 reads (Table 2). Thus, the shortness of sequences is the main disadvantage in this scenario when using 454 pyrosequencing technologies. Concerning on the second point, short query sequences might not be able to have sequence matches with a significant low E-value and high score, thus resulting in a high level of false negative results. In addition, Ranunculales is much less studied when compared with core eudicots and monocots, suggesting that a lot of Ranunculales lineage-specific genes might not be included in the database. As to the third point, EST fragments are directly sequenced instead of being cloned and sequenced, which leads to the truly random sequencing of all the expressed genes and might facilitate the discovery of rare genes. Emrich et al [32] found that coupling of LCM (Laser Capture Microdissection) and 454 sequencing technologies could facilitate the discovery of rare, possibly cell-type-specific transcripts. In summary, the nature of shortness in reads generated by 454 GS-FLX sequencing technology is the origin of low efficiency of annotation. However, it is still possible that new genes and/or Ranunculales- or even Epimedium- specific novel sequences could be discovered due to large-scale sequencing.

Potential applications of Epimedium EST dataset

Combining all the available information, we have a set of 22,295 GO-annotated Epimedium consensus sequences, which is very valuable and cost-effective for a non-model organism by sequencing one tissue at certain developmental stage. Additionally, 10.5%, 9.3%, 3.0%, and 0.2% of 12,242 Epimedium GO-annotated consensuses for biological processes are related to metabolic process, carbohydrate metabolic process, amino acid and derivative metabolic process, and secondary metabolic process, respectively (Figure 1A). Many bioactive components of Epimedium might derive from these biological processes (Figure 1A). A total of 13,807 Epimedium consensuses is annotated for molecular function, 2.2% and 20.6% of which have transcription regulator activity and catalytic activity, respectively, suggesting that these genes might take part in regulating and catalyzing the production of bioactive components (Figure 1B). All of these sequences can provide fundamental information for the metabolic-engineering and molecular breeding of Epimedium.

In details, this large EST dataset of Epimedium is a powerful resource for further studies. Firstly, it can facilitate cloning genes of interest, especially for genes involved in controlling metabolic pathways such as the flavonoid pathway. The Epimedium EST dataset generated in our study will provide enough sequence information for designing primers to obtain full-length genes. Therefore, these sequences will promote the molecular study of secondary metabolism in medicinal Epimedium plants. Secondly, molecular markers can be exploited based on EST sequences for researches in the fields of biodiversity, taxonomy, and population genetics. Although genomic SSRs have been developed from Epimedium, they have much lower transferability among different Epimedium species than EST-SSRs [9–11]. Therefore, EST-SSRs can further promote the taxonomic study of Epimedium species, especially those closely related Chinese species. Thirdly, it is noticeable that some EST-based molecular markers might be linked to genes involved in controlling traits that are important for market value [17, 33]. Thus, EST-SSRs can provide a powerful tool for genetic mapping and MAS. In summary, the EST dataset reported here can be used for decoding the molecular mechanism of the flavonoid pathway, as well as lays a foundation for biodiversity, population genetics, genomics, and functional genomics of Epimedium species.

EST-SSR frequency and distribution in Epimedium transcriptome

A total of 2,810 perfect microsatellites over 12 bp have been identified from Epimedium EST dataset, i.e. 3.67% EST sequences possess single sequence repeats. Obviously, the frequency of SSRs detected in Epimedium is similar to that of other dicotyledonous species ranging from 2.65% to 16.82% [34]. The EST-SSR frequency is dependent on the following factors. Firstly, parameters used in exploring microsatellites affect the results dramatically, such as the repeat length threshold and the number of repeat unit. In the present study, mononucleotide repeat motifs are excluded in our strategy for identifying microsatellites. Secondly, the genome structure or composition may also account partially for the EST-SSR frequency, because the genetic-biochemical background of the cells plays important roles in fixation of de novo-generated SSRs [35]. Thirdly, different softwares used in detecting SSRs might also affect the frequency. Some SSR identifying tools can find out imperfect SSR, for example, Sputnik [36], while others can only identify perfect SSRs such as SSRIT [37] and MISA [20, 21, 26].

Theoretically, the frequency of di-, tri-, tetra-, penta-, and hexanucleotide repeats should be in turn decreased. However, trinucleotide repeat unit is the most dominant SSRs, followed by di-, tetra-, hexa-, and pentanucleotide repeat units (Table 3). It is consistent with the results reported for other plant species [20, 23, 38–40]. As shown in table 4, the most dominant di- and trinucleotide repeat motifs are AG/CT (19.3%) and AAG/CTT (23.6%), respectively. These results are in agreement with that of Morgante et al [36], Kumpatla et al [34], and Toth et al [35]. Interestingly, there are only three CCG/CGG trinucleotide repeats in Epimedium, which is the most predominant in monocots [36, 38]. Our results strongly support and extend the notion of rarity of the CCG/CGG repeat units in a large number of dicotyledonous plants (Citrus, Coffea, Medicag, Glycine, etc.) [34]. These might be resulted from the high GC content and consequent codon usage bias in monocots [36, 38]. Selection might also be against the CCG/CGG repeat unit because of the requirements of the splicing machinery [41]. Long CCG/CGG sequences could compete splicing machinery components and give rise to inadequate splicing. Moreover, CCG/CGG repeats, might form potential higher structures such as hairpin and quadruplex, and thus affect the efficiency and accuracy of splicing and influence the formation of mature mRNA [35, 42].

Primer designing for EST-SSRs

In this study, attempt of diagnosing transferability have been performed by 32 EST-SSR primer pairs designed according to E. sagittatum ESTs. Unfortunately, seven and four primer pairs fail to transfer in Epimedium and basal eudicot, respectively, which might be due to primer(s) across splice sites, the presence of large introns in genomic sequence, or primer(s) derived from a chimeric or questionable consensuses sequence [39]. Some primer pairs succeed in amplifying in some species while fail in other species, suggesting that null alleles exist. Null alleles might be resulted from some mutations, including the deletion of microsatellites, and indels or substitution in primer binding sites [39]. In addition, the observed sizes of PCR products in Epimedium species deviate frequently from the expected product sizes according to the EST sequences. For instance, polymorphic bands amplified by EsESP05, EsESP25, and EsESP29 are larger than the expected sizes. It might be due to the presence of introns in the corresponding genomic regions [27, 39] and/or variation in the repeat numbers.

The EST-SSR transferability in genus Epimedium

In the present study, transferability of EST-SSRs derived from E. sagittatum is high (85.7%) across 52 species in genus Epimedium. The high transferability of EST-SSRs is also reported in other plants. For instance, 90% EST-SSRs in Citrus clementia can be transferred successfully to ten Citrus species and three related genera [20], and 89% EST-SSRs in Medicago truncatula succeed in transferring to eight genotypes [23]. About 92% Festuca arundinacea EST-SSR primer pairs can produce characteristic SSR bands in at least one of the six related species [27]. Numerous studies have also demonstrated that EST-SSRs possess higher transferability due to the conservation characteristics of EST-SSRs when compared with SSRs derived from genomic libraries [43, 44]. In genus Epimedium, the transferability of genomic SSRs derived from E. sagittatum [9] and E. koreanum [10] ranges from 21.4% to 42.9% and from 42.1% to 47.4%, respectively. Genomic SSRs isolated from E. brevicornu also possess low transferability ranging from 33.3% to 66.7% across other four medicinal Epimedium species [11]. Comparing with the high transferability of EST-SSRs (85.7%), the transferability of Epimedium genomic SSRs is remarkably lower.

Plant materials

A total of 52 Epimedium species including 55 individuals (Additional file 3) were selected and investigated for the transferability of EST-SSRs derived from E. sagittatum. Leaf samples were collected from Wuhan Botanical Garden (Wuhan, Hubei Province, China) and Garden Vision Nursery (Boston, MS, USA). According to the taxa system reported by Stearn [45], samples cover two subgenera and three sections, and represent most of Epimedium germplasm.

cDNA preparation

The red-magenta fully expanded leaves of E. sagittatum Maxim. grown in Wuhan Botanical Garden, Chinese Academy of Sciences, were harvested and prepared for RNA isolation. Total RNA was isolated according to the instruction of TRIzol kit (Invitrogen) and then purified to exclude the tRNA and rRNA and to enrich mRNA, according to the instruction of mRNA purification kit (Promega). mRNAs were reverse-transcribed by Powerscript™ II (Takara) with PCR primers SMART IV™ Oligonucleotide (5'-AAGCAGTGGTATCAACGCAGAGTGGCCATTACGGCCGGG-3') and CDS III/3' PCR Primer (5'-ATTCTAGAGGCCGAGGCGGCCGACATG-d(T)30_N-1N-3'). Long Distance PCR for double strand cDNA amplification was performed with LA Taq enzyme (Takara) for 25 cycles (95°C for 30 s, 68°C for 8 min) according to the SMART™ cDNA Library Construction Kit User Manual. Finally, double strand cDNAs were purified using DNA purification kit (Qigen) to generate high quality cDNAs.

EST processing, assembly, and annotation

Around 10 μg cDNAs were used for a half-plate sequencing run using 454 GS-FLX pyrosequencing technology. The sequencing was done at the Oklahoma University genome center. Considering the shortness of reads generated by 454 sequencing technology, full-length cDNA was sheared stochastically prior to sequence in order to randomly represent all transcripts. Sequence data for Epimedium Short Read Archive (SRA) described in this paper can be found at the public database [NCBI: SRA008151]. A total of 228,768 raw sequences was obtained. Lucy software [46] with default settings, except that minimum good sequence length cut-off was set to 50 bp, was used to remove low quality regions. SeqClean software http://compbio.dfci.harvard.edu/tgi/software/ was then used to eliminate poly A/T sequences, vector sequences, adaptor sequences and bacterial sequences. The resulted high quality sequences were assembled to consensus sequences using CAP3 software [47] with parameters as -p 95 -y 15 -s 251 -o 25. All consensus sequences were compared against Pfam database for protein domain analysis [48]. Based on the domain annotation, GO accessions of the consensus sequences were successfully assigned using Pfam2 go conversion file from the GO consortium http://www.geneontology.org/external2go. GO annotations were formatted for input into the GOSlim program and the output was parsed to count the occurrence of each GO category.

EST-SSR detection and primer designing

The EST-SSR detection was performed using Perl program MISA http://pgrc.ipk-gatersleben.de/misa. Since it was hard to distinguish real mononucleotide repeats and single nucleotide stretch error generated by 454 sequencing, mononucleotide repeats were excluded from this study. The parameters were designed for identifying perfect di-, tri-, tetra-, penta-, and hexanucleotide motifs with a minimum of 6, 5, 4, 4, and 4 repeats, respectively. Primer 3 software [49] was used to design primers. The major parameters for designing the primers were set as follows: primer length ranging from 18 bases to 28 bases with 20 as the optimum, PCR product size ranging from 100 to 300 bp, optimum annealing temperature 60°C, and GC content from 40% to 70% with 50% as optimum. According to these parameters, 32 primer pairs (Additional file 2) were synthesized and were used to detect the transferability in genus Epimedium.

Survey of EST-SSR polymorphism

Total DNA was isolated from leaf samples according to the CTAB method [50]. PCR amplifications were carried out in a final volume of 10 μL containing 50 ng of genomic DNA, 1× PCR buffer [10 mM Tris-HCl (pH 8.4), 1.5 mM MgCl₂], 0.2 mM dNTPs, 0.2 μM of each primer, and 0.5 U Taq polymerase (Biostar). The PCR reactions were performed under standard conditions for all primers in a thermocycler (Eppendorf). The annealing temperature was fixed for all primer pairs at 60°C. After 5 min at 95°C, 35 cycles were carried out with 30 s at 95°C, 30 s at 60°C, 40 s at 72°C for extension, and a final extension step of 10 min at 72°C. The separation of alleles was performed on a 6% polyacrylamide gel. PCR products were mixed with an equal volume of loading buffer. The mixture was denatured at 95°C for 5 min before loading onto the gel. Gels were stained with silver nitrate following the protocol as described in Bassam et al [51], and 25 bp DNA marker (Beijing Yuanpinghao Biotech Co., Ltd) was used for calculating the length of EST-SSR amplicons. All primers were detected firstly in E. sagittatum to check if the primers can successful amplify PCR products. Then primers succeeding in amplifying bands were tested for transferability. For detecting whether E. sagittatum EST-SSRs transferred successfully or not, the PCR products were firstly tested using 2% agarose gel and then polyacrylamide gel.

Genetic diversity and data analysis

Arlequin [52] was used to calculate the genetic diversity and average allele number. Genetic distance of the SSRs genotype was calculated using the log-transformed proportion of shared alleles ('Dps') by MICROSAT [53]. Polymorphic information content (PIC) was calculated by PIC_CALC http://dl.getdropbox.com/u/695591/PIC_CALC.rar combination with GenAlex6 [54]. A dendrogram was constructed based on Neighbor-Joining method [55] in the PHYLIP package [56].