Introduction

The genus Diospyros (Ebenaceae) is a widely distributed heterozygous genus in tropic and subtropic areas of Asia, Africa, and America (Central and North) with complex ploidy levels ranging from diploid (2n = 2x = 30) to nanoploid (2n = 9x = 135) (Yonemori et al. 2000). Diospyros kaki Thunb. (classified by a prominent Swedish naturalist Carl Peter Thunberg) or Japanese persimmon is the most economically important climacteric fruit species (having varied levels of ploidy; 2n, 6n and 9n) within this genus. In 2007, gross production of persimmon was estimated to be about 2,340,000 tons, of which 89.8% was produced in China, one of the origins of Japanese persimmon (Wang et al. 1997; FAO 2008). Vitamins A and C constitute the major portion of the vitamins present in fresh persimmon fruit.

On the basis of proanthocyanidins (PAs) (colorless phenolic polymers known as useful agents for human health, which show brown coloration upon oxidation), persimmon are further classified into astringent (A)-type fruits and non-astringent (NA)-type fruits (Dixon 2005; Ikegami et al. 2009). A comparative analysis of catechin composition among the five Japanese persimmon demonstrated that epigallocatechin (EGC) is relatively lower in the non-astringent type persimmon (Suzuki et al. 2005). An AST/ast allele having allelotypes as astringent (A) and non-astringent (NA) controls the expression of the trait. Expression of homozygous recessive ast allelotype at the AST locus results in the non-astringent (NA) genotype (Kanzaki et al. 2001; Yamada and Sato 2002). Recently, Ikegami and her team (2007) isolated seven genes (PAL, C4H, CHI, F3H, F30H, ANS, and ANR) from an astringent-type cultivar using suppression subtractive hybridization. Transcription of PAL, C4H, 4CL, CHS, CHI, F3H, F30H, F3050H, DFR, and ANR genes is high until mid-August, and then declines in October in the astringent-type cultivar (Ikegami et al. 2005a, b). Pang et al. (2007) identified ethylene receptor genes (DkERS1, DkETR1, and DkETR2) homologous to Arabidopsis ethylene receptor genes (ERS1, ETR1, and ETR2) in D. kaki. A Myb transcription factor (DkMyb4) controls proanthocyanidin biosynthesis in persimmon fruit (Akagi et al. 2009). Overall, D. kaki has great potential for becoming a model for understanding important traits in tannins and flavonoid biosynthesis as a fruiting crop species.

Diospyros kaki fruiting genotypes have wide morphological, physiological, and molecular diversity. Dominant and transposon-based markers have been described for Diospyros, including RAPD, RFLP, IRAP, REMAP, SSAP, SRAP, and AFLP (Yamagishi et al. 2005; Guo and Luo 2006; Du et al. 2009a, b).

In recent years, using in silico approach through mining expressed sequence tags (ESTs) have become an effective way for developing molecular markers such as Simple sequence repeats (SSRs), SNPs, SSR-FDMs for developing the saturated genetic linkage maps for various plant species (Hytena et al. 2010). In addition, the markers so developed not only exhibit higher level of intragenic transferability but also transferability to other closely related genera and may serve as potential markers for species discrimination, evolutionary inference and comparative genomics (Varshney et al. 2005). Extensive analysis has been done using ESTs available in the publicly available databases to identify genes temporally or spatially regulated during fruit growth and development in tomato, grape and apple (Fei et al. 2004; Da Silva et al. 2005; Park et al. 2006).

One of the most recently discovered regulatory mechanisms is post-transcriptional and involves 21–24-nt small RNA molecules (sRNAs). Micro RNAs (miRNAs) are non-protein coding, genomic derived small RNAs that participate in regulation of gene expression at a post-transcriptional level. In plants, they are involved in development, responses to biotic and abiotic stress and whilst some appear unique to a species, a large number of miRNA families are highly conserved across a wide range of plant species (Yang et al. 2007; Jian et al. 2010). miRNA transcripts are capped, spliced, polyadenylated and folded into long hairpin stem-loop precursor molecules (pre-miRNA), which are then processed by RNase III enzymes (Dicer-like in plants; DCL1) to form shorter hairpin primary miRNAs (Zhang et al. 2006). A 21–23 base pair double-stranded miRNA: miRNA* duplex is produced by further action of the Dicer enzyme and transported to the cytoplasm where the single-stranded mature miRNA is used as a template for target mRNAs silencing with complementary sequences by cleavage or translational inhibition by an RNA-induced silencing complex (RISC) (Bartel 2004; Zhang et al. 2006). The high sequence conservation of mature plant miRNAs has led to their successful prediction from sequence data using homology-based approaches (Zhang et al. 2005; Sunkar and Jagadeeswaran 2008).

We report here comparative mining of fruit cDNA libraries from different developmental stages of D. kaki. We have identified 506 SSRs primer pairs that can be further utilized for the inference of genetic diversity, species discrimination and studying the phylogeography of Diospyros genus and in particular D. kaki, a potential ortholog of miRNA159 in OYF library, which correlates the potential involvement of miR 159 family in development and relative distribution of SNP and SSR-FDMs markers.

Materials and methods

Sequence source and assembly

Diospyros kaki ESTs sequences were downloaded from GenBank (dbEST http://www.ncbi.nlm.nih.gov/dbEST) to give a total of 5,053 D. kaki ESTs from OYF library and 4,404 D. kaki ESTs from MF library. Mature miRNA sequences for all plant species were retrieved from the miRBase Registry (Release 14, September 2009, http://microrna.sanger.ac.uk/) and were used to generate a non-redundant reference set of 1,064 mature miRNA sequences. EST sequences were clustered using CAP3 program to prepare a tentatively consensus (TC) set (Huang and Madan 1999). To compare the relative richness of gene diversity sampled from each library, library-specific contigs, and singletons were compared.

SSR identification

The identification of SSR containing ESTs was carried out using in-house written program in C, which gives perfect as well as compound SSRs. Repeat patterns ranging from mono- to hexa-nucleotide were identified and systematically analyzed. The parameters defined for the identification of simple sequence repeats were seven minimal repeats for di-, five minimal repeats for tri-, four minimal repeats for tetra- and penta-, and three minimal repeats for hexa-nucleotide. The minimal length of mononucleotide simple sequence repeat was fixed at 14 bp. The poly A and poly T repeats were not considered as SSRs as they exemplify the 3′ end of mRNA/cDNA sequences, thus they were removed. Compound microsatellites were defined as repeats interrupted by a non-repetitive stretch of a maximum of 100 nucleotides.

SNP identification

Expressed sequence tags sequences were trimmed and a redundancy-based method for SNP confidence measurement, combined with SNP co-segregation (an independent confidence measure) was used to mine SNPs (Barker et al. 2003). The co-segregation score is a measure of whether a predicted SNP contributes to the definition of haplotype. The transition (Ts) versus transversion (Tv) ratio was also calculated for both the libraries to find the DNA substitution dynamics in the D. kaki genome.

Locus-specific primer designing and prediction of SSRs in open reading frames to identify relative biasing

Primer 3 software was being used to design a pair of primers flanking each SSR. The following parameter were used while designing the SSRs primers—optimum primer size was set to 20 where the range was between 18 and 27, optimum annealing temperature was set to 60.0 (the range was between 57.0 and 63.0), and the range of GC content was 20–80% (Shanker et al. 2007). Custom scripts and the standard genetic codes were applied to predict ORFs for all SSR-ESTs. SSR-ESTs were translated in all six ORFs and the longest fragments uninterrupted by stop codons were taken as the putative encoding segment (ORF) of the query SSR-ESTs sequences.

Annotation of SSR containing sequences, GC3 biology and gene ontology

Functional annotations of the SSR-ESTs sequences were determined on the basis of similarity using BLASTX program, available at NCBI (http://www.ncbi.nlm.nih.gov/blast) against non-redundant (nr) protein database entries and the best matches (E value <10−10) were compared to terms of the Gene Ontology (GO) Consortium (The Gene Ontology Consortium 2000). The resulting proteins obtained through similarity search by BLASTX were allotted to their respective classes. Using GO/UniProt comparison tables, candidate GO assignments were predicted on the basis of EST matches to the UniProt reference sequences.

Coding sequences of four additional Ericales species, such as: Actinida deliciosa, A. chinesis, Vaccinium corymbosum and Camellia sinensis were obtained from NCBI; in-house C++ code was used to compute position-specific nucleotide composition. In case of D. kaki, open reading frames and corresponding proteins were predicted using the assembled contigs of ESTs and nucleotide composition and sequence length was computed for each of the two EST libraries separately.

Using Gene Ontology (The Gene Ontology Consortium 2000) annotation of Arabidopsis thaliana (available at http://www.arabidopsis.org), all D. kaki protein sequences were aligned to A. thaliana using NCBI blastp with E value cut-off of 10−30, and the GO annotation of the best hit was used to annotate D. kaki genes. Chi-squared test (α = 0.05) was used to identify significant enrichment of different GO categories in high- and low-GC3 genes (Tatarinova et al. 2010). Categories were assigned on the basis of biological, functional, and molecular annotations available from the GO website (http://www.geneontology.org/).

Identification of functional domains markers (SSR-FDMs)

Using a python script sequences were translated into all six reading frames. In addition, Inter pro scan tool was used to analyze protein domain maintaining default parameter value (Quevillon et al. 2005; Yu et al. 2010). The sequences that contained both SSRs and functional protein domains were selected as SSR-FDMs; however, absence of predicted protein (as non-functional domain) caused exclusion for the sequences from further analysis.

Homology search and secondary structure prediction for miRNA identification

Candidate miRNA precursor sequences within the EST data were identified using BLAST and MFOLD RNA folding algorithms with parameters described elsewhere (Nasaruddin et al. 2007). Briefly, standalone BLAST (ver. BLAST-2.2.16) was used for local alignment of the EST against the non-redundant query set of 1,064 plant mature miRNA sequences. Default settings were as described elsewhere (Zhang et al. 2005). ESTs sharing homology with miRNAs in the reference set were defined as those containing a predicted mature miRNA with less than four (<4) mismatches compared to a known mature miRNA sequence in the reference set. Putative miRNA orthologs were analyzed using MFOLD RNA folding program (http://mfold.bioinfo.rpi.edu/cgi-bin/rna-form1.cgi) and candidate precursor miRNA (pre-miRNA) were filtered using the characteristics described elsewhere (Zuker 2003; Nasaruddin et al. 2007; Qiu et al. 2007; Xie et al.2007). Briefly, (1) the composition of the RNA sequences needs to be folded into a hairpin structure as per the stem-loop precursors. According to this process, each arm of the hairpin will contain ~22 nt mature miRNA sequences; (2) the lower minimal-free energy (MFE) and minimal-free energy index (MFEI) should be compulsorily present in the predicted secondary structure of the miRNA precursors than the tRNA or rRNA; (3) 30–70% of A + U content should be present in the predicted mature miRNA; (4) the mature miRNA sequence is the integral part of the hairpin loop segment. This mature miRNA should have less than six mismatches to the opposite miRNA* sequence of the other arm; (5) any part of mature miRNA:RNA*dimer loop or bulge should contain three nucleotides (maximum). This nucleotides should not be involved in canonical base pairing.

Prediction of miRNA targets

Potential targets of strong candidate miRNA from D. kaki EST were anticipated using RNA hybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/) (Rehmsmeier et al. 2004). The mature miRNA sequence was used to query the complete EST dataset using the following parameters: helix constraint (−f) of 8–12; maximum internal loop size (−u) of one and maximum bulge loop size (−v) of one (Rehmsmeier 2006). Good candidates were considered those with a negative folding free energy (MFE; ∆Kcal/mol) value below 70% of the MFE value for perfect complementarily and with end overhangs of no more than two nucleotides (Alves et al. 2009). Function of ESTs were predicted using BLASTX program by comparing the sequences against the non-redundant NCBI protein database with a cut offE value of 10−4 and 40% minimum identity score.

Results and discussion

Sequence assembly

Expressed sequence tags, represent partial and redundant cDNA sequences, making it difficult to analyze them effectively for putative mining of markers. To construct longer and less redundant sequence sets, we assembled ESTs by library, using CAP3 (Huang and Madan 1999). In OYF library, clustering of 5,053 sequences yielded 658 tentative consensus (TC) sequences with 3,117 sequences remained unclustered. In MF library, clustering of 4,404 sequences yielded 604 tentative consensus (TC) sequences with 1,925 sequences as singletons. The average length of the tentative consensus (TC) was 521 and 368 bp in OYF and MF library respectively. The diversity in ESTs libraries was confirmed by diversity index depicting higher degree of transcript diversity (Table 1).

Table 1 Summary of in silico mining of Diospyroskaki cDNA libraries for assembly and repeat analysis
Full size table

Screening, frequencies, primer designing and annotation of D. kaki SSRs-ESTs

In the present study, library-specific tentative consensus (TC) set of D. kaki were mined for SSRs with a minimum length of 14 bp. A total of 407 and 325 SSRs were detected in the OYF and MF libraries respectively; excluding poly A and poly T. Poly (A/T) were excluded (Karaoglu et al. 2004). In the OYF library, 5,053 sequences represent 407 SSRs with an average density of one SSR per 2.92 kb whereas in MF library from a number of 4,404 sequences screened only 325 SSRs were detected demonstrating average density of one SSR per 3.16 kb. The frequencies of SSRs with mono-, di-, tri-, tetra-, penta- and hexanucleotide repeat units are shown in Table 1. The most frequent repeat type found among different developmental libraries analyzed were di-nucleotide repeats (53.8%; 53.6%) followed by tri-nucleotide (30.2%; 24.8%), hexa-nucleotide (6.8%; 12.8%), tetra-nucleotide (6.2%; 4.4%), penta-nucleotide (1.8%; 1.0%), and mono-nucleotide repeats (1.2%; 0.7%), respectively (Fig. 1).

Fig. 1
figure 1

Frequency distribution of library specific repeat types identified in Diospyros kaki

Full size image

In both cDNA libraries surveyed, the mono-nucleotide repeats were relatively low when compared with other repeats. We further analyzed the observed abundant dinucleotide and trinucleotide repeat patterns (Figs. 2, 3) and reduction in the frequency of SSRs before and after assembly is both the libraries (Table 2). Similar patterns have been observed in the mining of the EST-SSRs markers in cereal species (Varshney et al. 2002). In case of the OYF library, out of 407 SSRs detected, primers could be designed only for 286 (70.2%) SSRs, whereas for the MF library, out of 325 SSRs detected, primers could be designed only for 220 (67.6%) SSRs (Supplementary Table 1). SSRs with primer pairs with respect to ORF were predicted in both MF and OYF libraries. In OYF library, out of 407 SSRs identified, 220 SSRs with primer pairs (54.0%) were found with respect to ORF. In the MF library, out of 325 SSRs identified, 153 (47.0%) SSRs with primer pairs were found with respect to ORF (Fig. 4).

Fig. 2
figure 2

Frequency distribution of library specific dinucleotide repeat types identified in Diospyros kaki

Full size image
Fig. 3
figure 3

Frequency distribution of library specific trinucleotide repeat types identified in Diospyros kaki

Full size image
Table 2 Frequency of SSRs in the EST sets before and after assembly
Full size table
Fig. 4
figure 4

Relative distribution patterns of SSRs with primer pairs in ORF and Non-ORF among di-, tri-, tetra-, penta- and hexa-nucleotide repeat in OYF and MF library

Full size image

Genic as well as the intergenic regions displayed the presence and absence of SSRs (Katti et al. 2001). A higher percentage of dinucleotide repeats in non-ORF regions may reflect natural evolution to maintain the conservation of functionality of all genes and their products (Fig. 4). Nevertheless, recent studies indicate that SSR expansions and/or contractions in protein-coding regions may cause a gain or loss of gene function through frame shift mutations (Fondon et al. 2008).

Simple sequence repeat-containing sequences (SSRs-ESTs) were annotated against the non-redundant (nr) protein database available at http://www.ncbi.nlm.nih.gov. Out of 296 MF library-derived SSRs-ESTs, 131 were found to homologous (44.3%) whilst for the OYF library, out of 352 SSRs-ESTs, homologs were available for only 146 (41.5%) sequences.

SNP identification

Redundancy-based SNPs mining resulted in identification of 68,067 SNPs and 4,273 indels in the D. kaki transcriptome. SNPs occurred at a frequency of one out of every 10 bp and indels at one in every 152 bp. A total of 28,232 transitions and 39,835 transversions were reported in this study. For explaining the nucleotide substitution dynamics, transition (Ts) to transversion (Tv) ratio was calculated because it provides insights into the process of molecular evolution. The transition/transversion ratio is relatively low in the OYF library (0.69) compared to the MF library (1.39) the overall transition (Ts) to transversion (Tv) ratio is 0.70, which indicates an relative increase of transversion (Tv) over transitions (Ts) (Table 3).

Table 3 SNP analysis
Full size table

Earlier studies have demonstrated higher rate of transitions over transversions due to abundant hypermutable methylated dinucleotides (5′-CpG-3′) (Ching et al. 2002; Strandberg and Salter 2004; Newcomb et al. 2006). However, in the present study, transversions prevail over the transitions. Neighbouring nucleotide effect demonstrates that the probability of transversion increases when the number of purines increases at the immediate adjacent sites (Zhongming and Boerwinkle 2002). Similar patterns of transversions over transitions were observed for genes on rice chromosome 8 (Wu et al. 2004).

In plant chloroplasts, an increase in transversions with increase in the A + T content of adjacent nucleotides has been observed (Morton 1995). These studies illustrate that the transition bias is not universal and supports the findings of the present study. However, the frequency of SNPs detected in this study is higher than the frequency of EST derived SNPs generally reported in earlier studies; 1 SNP/61 bp in Zea mays, 1 SNP/540 bp in Triticum aestivum, 1 SNP/123 bp in Sorghum bicolor and1 SNP/58 bp in Secale cereale transcriptome (Ching et al. 2002; Somers et al. 2003; Hamblin et al. 2004; Varshney et al. 2007). Possible reasons for variation in SNP density may perhaps be due to dissimilarity in the quantity of data analyzed.

Functional domains markers (SSR-FDMs)

Tentative consensus (TC) from the respective libraries was analyzed for functional domain markers excluding the mononucleotide repeats from this analysis. The translation of the sequences was performed in all six reading frames. InterProScan tool was used to analyse the resulting amino acid sequences from the longest reading frame (http://www.ebi.ac.uk/Tools/pfa/iprscan/). In the case of the OYF library, four potential SSR-FDMs were observed and Vps4 oligomerisation domain and the C2 calcium/lipid-binding domain were identified as major functional domains. The MF library displayed 10 potential SSR-FDMs but Basic-leucine zipper (bZIP) transcription factor, Glycoside hydrolase were observed as major functional domains. Therefore, this strategy not only implicates the evaluation of SSR polymorphisms, but also predicts function viability of these marker sequences. Association between candidate functional markers and trait of interest can be investigated by mapping SSR-FDMs.

miRNA and miRNA target identification

After removal of redundant EST sequences, a total of six ESTs from the MF library and seven from the OYF library were found to align with a known mature miRNA from the plant reference set with fewer than 4 mismatches within the mature miRNA sequence (Table 4). Of these candidates, one EST fulfilled the criteria for miRNA precursors based on the MFE (−137.40 kcal/mol) and secondary structure as predicted by MFOLD RNA folding program (Zuker 2003). A potential ortholog of miR159 was identified from OYF library (Fig. 5). A discovery rate of one miRNA precursor from a set of 9,457 ESTs of D. kaki lies within the expected range as per previous reports, which ranges from 0.83 per 10000 EST reported for Malus domestica to 1.69 per 100,00 EST from Gossypium hirsutum (Qiu et al. 2007; Gleave et al. 2008).

Table 4 ESTs from Libraries for ovary and young fruit and for mature fruit less than four mismatches to mature miRNA from rice or from Arabidopsis
Full size table
Fig. 5
figure 5

The predicted secondary structure of candidate miR159 precursor from Diospyros kaki. The secondary structure for EST DC584412.1 (candidate miR159) predicted using MFOLD (Zuker 2003). The sequence encoding the predicted mature miR159 is indicated by the line below the strand in which the mature miRNA is located

Full size image

Possible targets of the potential D. kaki miR159 (EST DC584412) were identified using RNA hybrid (Table 4). The candidates were screened and those having an MFE of −29.7 kcal/mol or lower (i.e. a minimum of 70% of the MFE for a perfect match) and two or fewer nucleotide overhangs at either end of the duplex were selected. A perfect match for the mature miRNA159 sequence has a predicted MFE value of −42.4 kcal/mol using RNA hybrid. Whilst six ESTs were identified as potential targets of predicted D. kaki miR159, only one of these (EST DC590670) matched an identified protein glutathione S-transferase (Table 5).

Table 5 Predicted miRNA159 targets and functions
Full size table

Previous reports show that miR159 targets include MYB transcription factors (Jones-Rhoades et al. 2006; Mallory and Vaucheret 2006), however, more recent studies suggest that, along with other miRNA families that are highly conserved across plant species, targets of miR159 are involved in diverse biological processes including gametogenesis, anther development, gibberellins signaling, and ethylene biosynthesis (Alves et al. 2009). Thus, the six EST identified as possible targets in this study may also represent a similar range of diversity.

GC3 biology

Unimodal GC3 profile of D. kaki CDS is typical for the Ericales order and other dicot plants (Fig. 6) (Tatarinova et al. 2010). Both ESTs libraries (OYF and MF) have similar GC content (Table 6). However, the contigs assembled for the mature fruit ESTs library are, on average, approximately 100 nucleotides longer. In order to analyze dependence between GC3 and GO, we took 10% of highest and lowest genes by GC3 in D. kaki and four other Ericales genomes (A. chinesis, Actinidia delicosa, V. corymbosum and C. sinesis). According to GO classification, high GC3 genes are over-represented in stress response genes, kinases, transcription factors and located in apoplast, membranes and cell wall (Table 7). Low GC3 genes are over-represented in genes involved in protein and nucleotide binding and located in nucleus, cytosol, and cytoplasm.

Fig. 6
figure 6

GC3 distribution for selected representatives of the Ericales order

Full size image
Table 6 GC variations across the different libraries of Diospyroskaki
Full size table
Table 7 GO-term enrichment for high- and low-GC3 Diospyroskaki genes
Full size table

The present study was aimed to generate resources that can be utilized for the identification and characterization of D. kaki germplam. The markers identified here can be used for subsequent prediction of germplasm diversity, phylogeography and species discrimination among the Diospyros genus.