Discovery of Simple Sequence Repeat Markers through Transcriptome Analysis of Baccaurea motleyana

Baccaurea motleyana (rambai) is underutilized fruits that are native to Malaysia, Indonesia and Thailand. In this study, a total of 54,779 unigenes identified from rambai transcriptome were used for simple sequence repeat (SSR) analysis by MIcroSAtellite (MISA). A total of 20,420 SSRs were found to be distributed within 37.27% of the total number of unigenes. Mononucleotide repeats represented the main type, accounting for 64.04%, followed by trinucleotide repeats (20.28%) and dinucleotide repeats (19.94%). Gene annotation to seven databases has a success ratio of 68.53% (National Center for Biotechnology Information (NCBI) protein sequences), 53.68% (NCBI nucleotide sequences), 27.43% (Kyoto Encyclopedia of Genes and Genome Ortholog), 56.0% (SwissProt), 52.44% (Protein family), 53.99% (Gene Ontology) and 26.44% (Kluster of Orthologous Group). Further rambai SSRs were randomly selected and validated to B. motleyana (rambai), B. macrocarpa (tampoi), B. polyneura (jentik-jentik), B. ramiflora (pupor) and B. scortechinii (setambun).


Introduction
Baccaurea motleyana (rambai) is an indigenous fruit species of Malaysia. It belongs to the Phyllanthaceae family, and the five species, namely B. macrocarpa (tampoi), B. polyneura (jentik-jentik), B. ramiflora (pupor), B. lanceolata (asam pahong) and B. angulata (tampoi belimbing) are found in abundance in the forest in Peninsular Malaysia, whereas B. angulate (tampoi) is an endemic to Borneo which can be found wild and also cultivated [1]. Currently, rambai is widely cultivated in the orchard, home gardens, and many of these trees are still found in the forest. Rambai trees are mainly planted from seeds.
They are valued for their shade and fruits, which are eaten raw, stewed or made into jam or wine. The fruit skin dries under the sun for 2-3 days and can be used in dishes for seasoning, especially in Malay community.
An effort in collecting and conserving rambai has been carried out by Industrial Crop Research Centre, Malaysian Agriculture and Research Centre (MARDI) to determine the distribution, diversity and elite accessions. At present, there are 155 accessions of rambai that have been characterized, and 25 accessions have been planted at MARDI Station Jelebu, Negeri Sembilan. There are some variations among the accessions based on morphological characters, especially fruit characters. Normally, those of selected cultivated fruit trees are sweet, sweet sour and pleasant to eat.
Rambai still lacks the required basic genomic resources for efficient breeding. Study of de novo transcriptome assembly of the arillode of B. amotleyana has been carried out, and the assembled contigs of transcriptome are deposited in the DDBJ TSA [2]. To the best of our knowledge, no report on rambai SSR has been published to date.
Simple sequence repeat (SSR) is arrays of short motifs that can be characterized by their hyper variability, abundance, reproducibility, Mendelian inheritance and codominant nature [3,4]. SSRs can be harvested from genomics and transcriptomes that are known as genomic SSRs and EST-SSRs. Since EST-SSRs derived from expressed sequences are more evolutionarily conserved compared to non-coding sequences, this type of SSR has a relatively high transferability than genomic SSRs [5]. Development of next-generation sequencing technologies assists researchers to harvest and evaluate many EST-SSRs from important crops such as Ipomoae batatas [6], Cicer arietinum [7], Epimedium sagittatum [8], and Elymus sibiricus [9].

Plant Material for RNA
Sample used for RNA extraction in this study was grown at Industrial Crop Research Centre, Malaysian Agriculture Research and Development Institute (MARDI), Selangor, Malaysia. The B. motleyana (rambai) apical meristems were collected and frozen in liquid nitrogen until further use.

RNA Extraction
The rambai shoot apical meristems were ground into fine powder with pestle and mortar that had been pre-cooled with liquid nitrogen. One milliliter of CTAB extraction buffer was added into ~100 mg tissue and incubated at 65 ºC for 30 min before being centrifuged at 16,000× g for 20 min. The supernatant collected was used for chloroform:isoamyl (24:1) extraction with equal volume (1 volume supernatant:1 volume chloroform:isoamy (24:1). In total 8 M LiCl (1/3 volume of supernatant), was added into the supernatant and incubated overnight at 4 °C. The solution was centrifuged at 16,000× g for 60 min at 4 °C. Pellets were washed with 70% ethanol and centrifuged for 15 min before being dried at room temperature. The pellets were resuspended in 20 μL nuclease-free water and treated with DNase I Recombinant, RNase-free water (QIAGEN, USA) to prevent the sample from DNA contamination. The quality and quantity of extracted RNA were examined using a 2100 Bioanalyzer (Agilent Technology, Santa Clara, CA, USA). The RNA sample with RNA integrity number > 7 was used for Illumina sequencing.

Plant Material for DNA Extraction
Five Baccaurea species, B. motleyana (rambai), B. macrocarpa (tampoi), B. polyneura (jentik-jentik), B. ramiflora (pupor) and B. Scortechinii (setambun) were collected at MARDI Station Jelebu, Negeri Sembilan and used for DNA extraction. Baccaurea leaves were DNA-extracted using a conventional plant DNA extraction method. Approximately 5 g of leaves were stuffed into a 2 mL tube with stainless steel balls and ground to a powder using Tissue Lyser (Qiagen, GmbH Germany). Next, 2,000 µL of the extraction buffer (final concentration: polyvinylpyrrolidone (PVP) 2%, diethyldithiocarbamate (DIECA) 4 mM, 5 mM ascorbic acid, 1.4 M NaCl, 100 mM tris-HCl (pH 8.0) and 20 mM ethylenedinitrilotetraacetic acid (EDTA)) was added to the ground powder. The mixture was incubated at 65 °C for 1 h with intermittent mixing. The same volume of isopropanol was added to this mixture, and the two were Discovery of Simple Sequence Repeat Markers through Transcriptome Analysis of Baccaurea motleyana 35 homogenized by inverting the plates. The mixture was centrifuged at 5,500 rpm for 15 min. Later the supernatant was transferred into an equal volume of isopropanol and centrifuged at 5,500 rpm for 15 min. The DNA was air-dried and resuspended in 50 μL tris-EDTA buffer; the DNA was quantified for concentration using Thermo Labsystems Fluoroskan Ascent TM (Thermo Scientific, USA); and DNA integrity was tested on 0.8% agarose gel for quality control.

PCR SSR Validation
Primers were successfully designed and oligo synthesized for 28 SSRs after SSR mining; these primers were used for PCR. The M13-tailed primer was used, which includes three primers: a forward primer (with an M13 tail sequence), a reverse primer, and a fluorescently labeled M13 IRDye primer [10]. PCR amplification was performed in a total reaction volume of 10 µL, consisting of 1.0 µL of template DNA (approximately 40 ng·µL -1 ), 1 µL PCR reaction buffer (1.5-3.0 mM MgCl 2 ), 2 mM dNTP mixture, 10 µmol of each primer set, 5 µmol of a fluorescent dye (FAM/VIC/NED and PET), and 1 unit of Taq DNA polymerase (Invitrogen, USA). Amplification was performed using a Peltier Thermal Cycler, DNA Engine Tetrad 2 (BioRad, USA). The PCR profile was an initial denaturation at 94 °C for 5 min, followed by 34 cycles of 1 s denaturation at 94 °C for 30 s, annealing (see Table 7) for 45 s, and 72 °C for 45 s extensions and finally 1 cycle of 72 °C for 7 min. The PCR products were genotyped using the 3730xl DNA Analyzer (Applied Biosystems, USA).

Sequence Data Analysis and De Novo Transcriptome Assembly
The original raw data from Illumina HiSeq TM were filtered first to generate high-quality data. The filtering included removal of adaptor contaminants, removal of reads with more than 10% of N, where N represented the base with uncertainty, and removal of low-quality base reads, with the quality value Q ≤ 20.
The high-quality transcriptome data were then assembled using Trinity [11] version r20140413p1 with parameter minimum k-mer coverage of two and others being default. De novo assembled transcripts from Trinity were clustered using Corset Version 1.05 to remove redundancy.

Molecular Marker Discovery
SSR markers were discovered using MIcroSAtellite (MISA), an identification tool with the following parameters: minimum SSR motif length of 12 base pairs (bp) and repeat lengths of mono-10, di-6, tri-5, tetra-5, penta-5 and hexa-5. The maximum size of interruption allowed between two different SSRs in a compound was 0 bp. Identified SSR markers were then mined using a custom Perl script and R programming to examine the repeat distributions. Mining of SSR markers included choosing SSR markers that represented more than two alleles and represented a single contig. We extracted 200 bp of each side of repeat motif region using BEDTools Version 2.21.0. SSR markers were annotated using the BLAST program (blastx) against a non-redundant (nr) protein database in NCBI with an E-value of 1 × 10 -5 and BLAST2GO software (Supplementary excel file of rambai SSR).

Gene Functional Annotation
The functional annotation was outsourced to Novogen Company. The software and parameters used in each database were, NT using NCBI Basic Local Alignment Search Tool (BLAST version 2.2.28), the E-value threshold of 1.0e-5, non-redundant protein sequences (NR) using SwissProt and Cluster of Orthologous Groups (KOG) of Diamond 0.8.22. NR and SwissProt databases, the E-value thresholds are 1.0e-5 and 1.0e-3 for KOG. Protein family (Pfam), the prediction of protein structure domain was made using HMMER 3.0 package, hmmscan, the E-value threshold is 0.01. Gene ontology (GO) for the protein annotation results of NR and Pfam using Blast2GO v2.5 [12] and a Perl script, the E-value threshold is 1.0e-6 and Kyoto Encyclopedia of Genes and Genome (KEGG) using automatic annotation server, the E-value threshold is 1.0e-10.

Scoring and Data Analyses
The fragment files were analyzed using GeneMapper 5.0 (Applied Biosystems). GS500LIZ was used as the standard control. Allele peaks in the electropherograms were scored and analyzed [13]. The scored peaks were imported into a Microsoft Excel file, and the data were pre-analyzed using MicroChecker v2.2.3 [14] to check for the presence of null alleles, large allele dropouts, and scoring errors due to stutter peaks. The passed data, after MicroChecker analysis, were converted into an input data file using CONVERT 1.2 [13] for further analysis. POPGENE32 version 1.31 was used to calculate the number of alleles, allele frequencies of each variety, and the genetic distance between varieties. NTSYSpc software was used to construct a dendrogram using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) method based on the genetic distance produced from the POPGENE32 software.

Sequencing and De Novo Assembly of Rambai Transcriptome
High-throughput next-generation sequencing generated a total of 81,941,472 paired end raw reads of rambai transcriptome comprising 12.3 Gigabase (Gb) with an average read length of 150 base pair (bp). After quality assessment, approximately 96.5% of high-quality data were trimmed out based on the base quality and adaptor contaminants (Table 1). High-quality cleaned reads (79,121,624 paired end reads) were then assembled into 54,813 transcripts and 54,779 unigenes in the rambai transcriptome (Table 2). 7 100 a Raw data were the original raw reads from sequencer Illumina HiSeqTM 2500. b Clean data were the reads after quality filtering. The filtering included removing the reads with adapter contamination, removing the reads with more than 10% of uncertain nucleotides (N > 10%) and removing low-quality reads (base quality less than 20) with more than 50% low-quality base. c Percentage of high-quality data based on this formula: Total reads of clean data Total reads of raw data × 100 (%) 1,537 1,538 a Transcripts were the transcripts produced by transcriptome assembler Trinity. b Unigenes were the clustering transcriptomes that were not redundant. c The N50 length was defined as the length for which the collection of all transcripts of that length or longer contains at least half of the total of the lengths of the transcripts and for which the collection of all transcripts of that length or shorter contains at least half of the total of the lengths of the transcripts.
The mononucleotide repeats and compound formation were excluded from the analysis. The most frequent repeat motifs for trinucleotide were GAA/TTC, followed by dinucleotide repeats AG/CT and tetranucleotide repeats AAAT/ATTT (Table 4).

Gene Functional Annotation
To get comprehensive gene function of rambai genes, gene annotation to seven databases, namely National Center for Biotechnology Information (NCBI) non-redundant protein sequences (NR), NCBI nucleotide sequences (NT), Kyoto Encyclopedia of Genes and Genome Ortholog (KO), SwissProt, Protein family (Pfam), Gene Ontology (GO) and Cluster of Orthologous Groups (KOG), was used as reference.
The NCBI non-redundant protein sequences (NR), include protein sequence information from GenBank, Protein Data Bank (PDB), SwissProt, Protein Information Resource (PIR) and Protein Research Foundation (PRF). The NCBI nucleotide sequences (NT) are the nucleotide sequence database that includes nucleotide sequence from GenBank of the European Bioinformatics Institute (EMBL) and DNA Data Bank of Japan (DDBJ). KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as cell, organism and ecosystem, from molecular-level information, especially for large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. KEGG is an established Cluster of Orthologous (KO) annotation system that can accomplish the function annotation of the genome/transcriptome of a newly sequenced species. SwissProt is a manual annotated and reviewed protein sequence database that has a high-quality protein sequence database from experimental results, computed features and scientific conclusions. Pfam is comprehensive collection of protein domains and families, represented as multiple sequence alignments and as profile of hidden Markov models. Many proteins are composed of structural domains, and the protein sequence of a specific structural domain possesses a certain degree of conservative property. GO is the established standard for the functional annotation of gene products and controlled vocabulary used to classify the functional attributes of gene products of a biological process, a molecular function and a cellular component.
Gene annotation to GO has grouped rambai genes into three main GO domains: biological process, cellular component and molecular function. In the biological process, more than 75,000 genes were identified, and three highest gene functions were involved in the cellular process, the metabolic process and the single-organism process; in the cellular component, more than 46,000 genes were identified, and three highest gene functions were involved in cell, cell part and organelle; in the molecular function, more than 10,000 genes were identified, and three highest gene functions were involved in binding, catalytic activity and transport activity (Fig. 2).
The Cluster of Orthologous Groups (KOG) to rambai genes has been classified into 26 KOG groups (Fig. 3). The three highest groups with more than 9% of genes annotated were attributed to general function prediction, posttranslational modification/protein turnover/chaperones and translation/ribosomal structure/biogenesis. KEGG has annotated rambai genes into five KEEG pathways: the cellular process, environmental processing, genetic information processing, metabolism, and the organismal system (Fig. 4). The highest genes were in metabolism (6,741 genes), followed by genetic information processing (3,247 genes), the cellular process (706 genes), the organismal system (550 genes) and environmental information processing (458 genes).
Genes containing SSRs were annotated into the Cluster of Othologous (KO) (Fig. 5). In the cellular process, four major annotate genes were involved in endocytosis, peroxisome, phagosome and regulation of autophagy; in environmental information processing, three major annotated genes were involved in plant hormone signal transduction, the phosphatidylinositol signal system and ABC transports; in genetic information processing, three highest annotated genes were involved in protein processing in the endoplasmic reticulum, ribosome and spliceosome; in  The number of genes successfully annotated in all the seven databases 7,545 13.77 9

Discovery of Simple Sequence Repeat Markers through Transcriptome Analysis of Baccaurea motleyana 42
The number of genes successfully annotated in at least one database 41,547 75.84 10 Total unigenes 54,779 100 Table 6 Prediction pathway levels of gene-containing SSRs from rambai. The highest pathway level is from the cellular process, environment information processing, genetic information processing, metabolism and the organelle system with ten selected genes.
1st-level (pathway) 3rd-level (pathway) Description of gene-containing SSRs metabolism, three highest annotated genes were involved in carbon metabolism, biosynthesis of amino acids and starch/sucrose metabolism; in the organismal system, two major genes were involved in circadian rhythm and plant-pathogen interaction. Further, we made analysis to predict what genes were involved from the first level pathway of the cellular process, environment information processing, genetic information processing, metabolism, and the organismal system and the third level pathway of endocytosis, plant hormone signal transduction, protein processing in the endoplasmic reticulum, carbon metabolism and plant-pathogen interaction. A list of ten selected genes is provided in Table 6.

Rambai SSR Validation
Twenty-eight rambai SSRs were randomly selected and used for validation of five Baccaurea species, namely B. motleyana (rambai), B. macrocarpa (tampoi), B. polyneura (jentik-jentik), B. ramiflora (pupor) and B. cortechinii (setambun). Out of 28 SSRs, only four were polymorphic, and the call rate was more than 90% across Baccaurea species. The remaining microsatellites were monomorphic, exhibited low call rates, and produced multiple bands across the rambai species. Analysis of four microsatellites revealed a total of 52 alleles, which ranged from 8 (B9) to 21 (B89) alleles per locus and had a mean of 13.75 alleles per locus. The range of the polymorphic information content (PIC) value was from 0.5656 (B9) to 0.9402 (B89), with an average of 0.8087, whereas the value of gene diversity and heterozygosity ranged from 0.5851 (B89) to 0.9432 (B89) and from 0.1667 (B89) to 0.4545 (B89), respectively (Table 8). Dendogram analysis of 15 accessions from five Baccaurea species was classified into two main groups (Fig. 6)