Abstract: Teak is one of the highly famous woody plant species for its premier quality of wood. Teak has problem on productivity because of long reproductive cycle. The problem is basically related to mechanism of flower development. The aim of this study was preliminary development of expressed gene database to characterize the floral transcriptome in teak. Two subtracted cDNA libraries were constructed from vegetative and generative bud tissues. Libraries were sequenced using Illumina MiSeq technology which generated paired-end read sequences 3,778,316 for vegetative and 3,701,878 for generative. The sequences assembled de novo into 87,365 transcript contigs consisting of 42,435,728 bases with N50 of 498 bp using CLC-Genomics Workbench. 76,169 (87.18%) of the 87,365 assembled contigs exhibited significant similarity BLASTN to Solanum lycopersicum database (www.phytozome.com). The assembled contigs were annotated through high stringency BLASTX analysis to proteome of S. lycopersicum. Distribution of contigs abundance between vegetative and generative stages analyzed using the DEGseq approach. The numbers of contigs distribution are 24,730 in vegetative, 28,912 in generative and 33,723 in both stages. The functionally protein datasets characterized by Gene Ontology (GO) annotation and KEGG metabolic pathways assignments for the result of DEG analysis. These contigs, 18,756 (75.84%) from vegetative, 22,089 (76.40%) from generative and 22,917 (67.96%) from both stages were assigned to GO classes. A total of 1455 (13.77%) were mapped to 30 pathways from vegetative, 1,638 (13.70%) were mapped to 27 pathways from generative and 1,652 (12.20%) were mapped to 30 pathways from both by BLAST comparison against the KEGG database. The biological processes of flowering developments were identified in the biological process dataset and the numbers of contigs were discovered different between stages. This transcriptome dataset information will act as a valuable resource for further molecular genetic studies teak, as well as for isolation and characterization of functional genes involved in flowering development pathways.
INTRODUCTION
Teak (Tectona grandis, L. f) is a tropical tree species distributed naturally in countries including India, Myanmar, Thailand, Myanmar and Indonesia (Orwa et al., 2009; Palupi et al., 2010; Lyngdoh et al., 2010). Teak is one of the world’s premier hardwood tree species, highly famous for its quality, profile and durability of timber. In Indonesia, teak flowering usually appears every year at the beginning of the rainy season (October-November) and only few flowers (about 1%) develop into fruits. Fruits fall gradually during the dry season (Orwa et al., 2009). According to the fact, the main limitations of teak improvement are it has a long reproductive cycle and produces low seeds. Both problems are basically related to mechanism flower development (De Gyves et al., 2007; Rosli et al., 2009; Widiyanto et al., 2009; Palupi et al., 2010). Hence, the determination of the genetic pathways and identifying specific genes involved in teak flowering and flower development could be beneficial for teak productivity improvement. We are interested in studying more about the roles of genes that control development of flowers in teak especially during the transition period between shifts of the vegetative to reproductive phase. This study was preliminary of teak floral transcriptome characterization, before isolation and characterization of functional genes involved in flowering development pathways.
In this study, we sequenced the transcriptome of T. grandis using the next generation of high throughput paired-end RNA sequencing (RNA-seq) technology, Illumina MiSeqTM 2000 (Fan et al., 2013). Then, CLC bio bioinformatics technology tool was used to perform a de novo assembly and annotation without prior genome information (Agusti et al., 2011; Angeloni et al., 2011). This transcriptome database helped to reveal much about the functional genomics of T. grandis and was then used to predict the functional classification of many unigenes using GO and KEGG pathway analysis (Chang et al., 2013; Rosero et al., 2011). These results lay the foundation for understanding the relation between gene expression patterns and plant development, physiology and structure and will be helpful for the molecular approach to improve of T. grandis. Furthermore, we focused on the sequences that are related to flowering developmental biological process in the aim of exploring the relationship between genes in transition development vegetative to generative stage.
MATERIALS AND METHODS
Teak tissues materials and RNA isolation: Vegetative and generative stage shoot tips of teak were collected from a 12 years old teak plant in Institute Technology Bandung, Indonesia for RNA isolation. The following VS tissues were sampled from vegetative apical shoots. LB2 tissues were sampled from lateral (nodal) floral-buds 2nd of generative stage shoots (Fig. 1). Both of teak tissue samples were frozen in liquid nitrogen immediately upon collection and put in Dry Shipper for shipping from ITB-Indonesia to Pennsylvania State University (PSU)-USA. Samples were immediately frozen at -80°C upon arrival at PSU until use.
| Fig. 1(a-b): | Vegetative and generative stage shoot tips of teak, (a) Vegetative stage and (b) Generative stage |
Total RNA was obtained using the method for RNA isolation protocol that developed by Dr. Carlson’s team at Schatz Center Laboratory, PSU-USA. Frozen tissue were ground to a fine powder under liquid nitrogen and dispersed in CTAB buffer. Following 2 chloroform extractions, RNA was precipitated with LiCl2,again extracted with chloroform and precipitated with ethanol. The resulting RNA pellet was resuspended in 20-100 μL of DEPC-treated water (Barakat et al., 2012). RNA concentration analysis on a QubitTM fluorometer (www.invitrogen.com/qubit) to showed a total yield of RNA sample. The concentration of RNA are 555 and 206 ng μL-1 for VS and LB2 sample, respectively. The integrity of RNA was assessed with the Agilent 6000 RNA Nano Chip Kit on 2100 Bioanalyzer (Agilent Technologies).
Paired-end cDNA library preparation and MiSeq Illumina sequencing: Total RNA of teak was extracted from the two tissues using the protocol described previously. The double-stranded cDNA was synthesized using the cDNA Synthesis System using random hexamer primers (illumina) according to manufacturer’s instructions. The paired-end library was developed according to the protocol of the Paired-End sample Preparation kit (lllumina, USA) (Angeloni et al., 2011; Li et al., 2012; Fan et al., 2013). The resulting library was sequenced at Penn State University using Illumina MiSeqTM 2000 (Illumina Inc., USA).
Transcript assembly: Two sequence data in FASTQ files computed with CLCbio for transcript assembly strategy (Angeloni et al., 2011). Paired-End reads were trimmed for quality score and the presence of repeated sequences >50 bp using the modified Mott-trimming algorithm present (default parameters) in CLCbio. We assembled de novo the Illumina-trimmed paired-end reads into transcript contigs using the software CLC Genomics Workbench by setting minimum 95% identity, minimum 40% overlap and 200 bp as minimum contig length.
Contig annotation: The quality of the de novo assembly was assessed with a local BLASTn (e-value<10‾6) alignment of all the contigs against A. thaliana, P. tricocharpa, M. guttatus, S. tuberosum and S. lycopersicum (www.phytozome.com(Last time access on this date 2015-03-17)) using CLCbio Workbench. Top hit species results use for homology based annotations of teak (Rosero et al., 2011; Barakat et al., 2009, 2012).
DEGseq analysis: Comparison of Digitally Gene Expression (DEGseq) between vegetative and generative tissues was done using RNAseq analysis software test developed by CLCbio-Genomic Workbench. DEGseq analysis was used to identify flowering development genes in transcript abundance because it integrates several statistical methods (Barakat et al., 2009, 2012). The number of reads per contig for each gene was compared between vegetative stage as control and generative tissues in teak separately. Similar analyses were performed for gene orthologs from both tissues. Orthologs were identified using a reciprocal best hit approach. RNAseq employs a random sampling model based on the read count in vegetative and generative tissues libraries and performs a hypothesis test based on that model. Genes expression in vegetative, generative and both of them are identified and go to GO enrichment.
GO analysis: Further assessment of the quality of the de novo assembly was carried out as follows. We compared the depth and the length of contig coverage with reference to orthologous genes in S. lycopersicum and A. thaliana, by plotting the ratio of contig length to S. lycopersicum and A. thaliana orthologue coding region length against coverage depth. Orthologous genes were retrieved performing a local BLASTX alignment (e-value<10-6) using CLCbio Workbench with the TAIR9 A. thaliana database and S. lycopersicum predicted proteins (Unipro/Swissprot database). To further assess the coverage and the quality of the assembly, we used BLASTX to align the contigs to the manually curated protein database Uniprot/Swissprot using DAVID Bioinformatics Resources at http://david.abcc.ncifcrf.gov/ (Huang et al., 2007; Jiao et al., 2012). DAVID Bioinformatics is an automated tool for the assignment of Gene Ontology (GO) terms to BLAST hits and it has been designed for use with novel sequence data (Jiao et al., 2012), Assignment of GO terms to contigs with significant BLASTX match with Swissprot (http://www.expasy.ch/sprot/) and the KEGG pathway (http://www.genome.jp/kegg/) were also performed using DAVID Bioinformatics. In addition, we generated GO assignments for A. thaliana and S. lycopersicum annotated proteins to compare the distribution of functional annotation in teak to those plants species with a well-characterized transcriptome, we did the GO analysis for the result of DEG analysis from vegetative stage and generative stage of teak samples.
RESULTS
Illumina sequencing output statistics and reads assembly: T. grandis vegetative and generative cDNA libraries were constructed from a pool of RNA isolated from vegetative and generative bud tissues teak tree using the Illumina MiSeqTM 2000 system at Penn State University. A total of 3,778,316 and 3,701,878 reads were generated from vegetative and generative teak transcriptomes, respectively (Table 1). The average length of the reads was 151 nucleotides (Fig. 2). De novo contig construction of the Illumina reads using the CLCbio assembly software led to the construction of 87,365 contigs from combined vegetative and generative teak (Table 2). Those contigs were having an average length of 486 nt, 225 nt for minimum length and 4,361 nt for maximum length (Fig. 3).
| Table 1: | Summary statistics of cDNA library |
| Fig. 2: | Paired reads distance distribution |
| Fig. 3: | Histogram of the frequency of different contigs sizes in transcriptome assemblies of teak samples |
| Table 2: | Summary statistics of sequencing and de novo assembly results |
| Table 3: | No. of BLAST hit top species for homology based annotations of teak contigs |
Contigs annotations: Collection of 87,365 contigs, enriched in vegetative to generative transition stage related transcripts, was obtained from both vegetative and generative bud subtracted libraries. Top hit species for homology based annotations of teak contigs (Table 3) were: A. thaliana (74.36%), P. trichocarpa (8.56%), M. guttatus (18.91%), S. tuberosum (18.91%) and S. lycopersicum (87.18%). S. lycopersicum was the highest blast hits species for teak. This is a remarkable result when considering the current state of functional annotation of teak to the S. lycopersicum proteome database (www.phytozome.com(Last time access on this date 2015-03-17)).
Transcriptome comparison between vegetative and generative tissues: We compared the transcriptomes from teak vegetative tissues and generative tissues to gain insight into the differences in the gene activity of the transition vegetative to generative stages in teak development. This comparison showed that the distribution of contigs in vegetative stage, generative stage and both using DEG analysis software (Fig. 4).
Detailed comparison of the Gene Ontology (GO) transcriptomes in vegetative stage, generative stage and both showed the different percentage of biological processes, cellular component and molecular function (Fig. 5). Figure 6 showed top 25 of biological process, cellular component and molecular function that occurred in the tissue samples.
| Fig. 4: | Contigs distribution result of DEG analysis |
| Fig. 5(a-c): | Chart of GO categories of teak, (a) Vegetative, (b) Both and (c) Generative |
In the category of biological process, response to abiotic stimulus, phosphorus metabolic process and phosphate metabolic process comprised the largest proportion of sequences, accounting for (6.4% in vegetative stage, 5.53% in generative stage and 5.3% in both), (6% in vegetative stage, 6.3% in generative stage and 6% in both) and (6% in vegetative stage, 6.3% in generative stage and 6% in both) of the total, respectively. According to flowering development biological processes, there are post-embryonic development (4.4% in vegetative stage, 4.1% in generative stage and 4% in both) and reproductive developmental process (4.1% in vegetative stage, 4% in generative stage and 3.8% in both), reproductive structure development (3.8% in vegetative stage, 3.6% in generative stage and 3.4% in both) comprised part of the top ten largest proportion.
In the category of cellular components, plastid comprised (16.20% in vegetative stage, 13.7% in generative stage and 13.8% in both), chloroplast was (15.8% in vegetative stage, 13.3% in generative stage and 13.5% in both) and intrinsic to membrane (11.3% in vegetative stage, 11.6% in generative stage and 11.3% in both) these three subgroups were dominant over the others. In the category molecular function, sequences with the functions of nucleotide binding, metal ion binding and purine nucleotide binding comprised (17.4% in vegetative stage, 18.4% in generative stage and 17.3% in both), (14.7% in vegetative stage, 14.9% in generative stage and 14.5% in both) of the total (Table 4).
On the other hand we also identified the other biological processes of flower development. We compared the biological processes of flower development between vegetative stage, generative stage and both. In the flowering developmental biological processes, positive regulation of developmental process, positive regulation of flower development and regulation of meristem development comprised only in vegetative stage.
| Fig. 6(a-c): | Histogram of GO classification of teak, (a) Biological process, (b) Cellular component and (c) Molecular function |
| Table 4: | GO classification of teak |
The pollen tube development, tube development and negative regulation of flower development accounted only in generative stage (Table 5). Analysis of KEGG metabolic pathway assignments revealed that our contig catalog covers all major plant metabolic pathways with a certain dominance of plant hormone biosynthesis and many alkaloid biosynthesis, indicative that those pathways, seemingly paired in response to reproductive developmental process (Table 6).
| Table 5: | Biological processes of flower development |
| Table 6: | KEGG pathway |
Vegetative and generative stage of teak transcriptome sequencing and annotation: Next Generation Sequencing (NGS) technology during the last decade have dramatically impacted genome sequencing and transcriptome analysis (Fan et al., 2013; Fu et al., 2013; Fox et al., 2013). This technique could be used for model plants with known genome sequences and also has been successfully used to analyze the transcriptome in non-model plants (Collins et al., 2008; Logacheva et al., 2011; Li et al., 2012). However, this technique requires cDNA cloning and individual RNA preparations for each sample stages, is time consuming and very costly. Pyrosequencing like 454 and illumina plat form introduced recently constitutes a better alternative for transcriptomics (Lulin et al., 2012). The high number of reads generated per run together with the low sequencing error rate in the contigs obtained makes it a good tool to deeply sequence the transcriptome of plants. This approach has been used successfully for analyzing the transcriptomes of maize and Arabidopsis and have applied it to the non-model tree species Castanea dentata and C. mollissima (Collins et al., 2008; Barakat et al., 2009, 2012).
Tectona grandis, also known as teak is a tropical deciduous tree native from moist tropical forests of Asia (Palupi et al., 2010; Khanduri, 2012). Tectona grandis is lamiaceae family is known for the quality of its wood (Lyngdoh et al., 2010). Despite its ecological and increasing economic importance, very little is known about the biology of this species at the genetic, molecular and biochemical levels (Borges et al., 2008). Genomic tools have recently increased the numbers and volume of genomic resources for several crop plants and trees and have contributed to enlarge our knowledge on basic aspects of plant biology; furthermore, they represent valuable sources of candidate genes and new molecular markers to assist improvement programs (Collins et al., 2008; Eveland et al., 2010; Huang et al., 2012). Biological sequences reported to date in public databases and belonging to teak do not exceed 20 entries: This very narrow availability of genetic information is the main problem to initiate improvement programs in T. grandis.
Our study generated 3,778,316 and 3,701,878 reads and 87,365 high quality contigs from vegetative and generative teak transcriptomes, respectively. A fraction of teak contigs could be annotated using the S. lycopersicum proteome than Arabidopsis or the others (Table 2). Most of the genes in teak hits to the S. lycopersicum proteome encoded proteins annotated. Those genes could be homology to S. lycopersicum using the Blast algorithm. Over 80% of the teak reads could be annotated using the S. lycopersicum proteome. By taking into consideration the sequences that have homologies in the S. lycopersicum proteome, assuming that the two samples of teak have a similar gene number as S. lycopersicum. cDNA sequences generated from both teak samples cover various biological processes and molecular functions indicating that the technique constitutes a powerful tool for sequencing the transcriptome of non-model species. These results confirm that pyrosequencing constitutes a powerful tool for transcriptome characterization and gene discovery.
Transcriptome comparison between vegetative and generative tissues from Tectona grandis: Gene Ontology (GO) annotation analyses showed that, overall vegetative and generative tissues from teak present a similar transcriptome. Gene function categories associated with response to abiotic stimuli and metabolic process are highly represented in both transcriptomes. The second most highly represented category includes genes involved in reproductive development. The category represented the most is composed of genes associated with various reproductive processes as previously described in other systems such as Gerbera, Fagopyrum and Prunus. Detailed analysis of illumina sequences from both vegetative and generative tissue showed that the tagged genes included a large number associated with response to abiotic stimuli, metabolic process and reproductive development.
These include genes involved in regulation of development, meristem development and reproductive development genes. Comparison of flowering developmental genes highly expressed in the vegetative and generative tissues of teak showed that a fraction were either preferentially expressed in vegetative or in generative stage. Genes of positive regulation of developmental process, regulation of meristem development, meristem development, shoot development, shoot system development, positive regulation of flower development, embryonic meristem development and embryo sac development represented the functional category with the largest number of reads in vegetative stage. Genes of reproductive developmental process, reproductive structure development, regulation of flower development, flower development, gametophyte development, pollen development, pollen tube development, tube development, negative regulation of flower development, fruit development and seed development complied the largest number of reads in generative stage.
Positive regulation of developmental process, regulation of meristem development and positive regulation of flower development genes category expressed only in vegetative stage. Pollen tube development, tube development and negative regulation of flower development genes category founded only in generative stages.
These different suggest that these tissues may modulate the expression of flowering development genes in transition vegetative to generative in teak. The important thing after this step is select the candidate genes involved in regulation of teak vegetative to generative transition. Overall, this study allowed us to conclude that teak tree responds to abiotic stimuli before entering to flowering developmental stage. The different category of flowering developmental processes between vegetative and generative stage showed us the regulation of transitional vegetative to generative.
CONCLUSION
In conclusion, this study allowed us to (1) Obtain 87,365 contigs from vegetative and generative tissue of teak, (2) Transcriptomes of teak could be annotated using the S. lycopersicum proteome according to BLAST result, (3) Compare the transcriptomes of vegetative and generative tissues of teak in flowering developmental stage and (4) Identify potential biological processes involved in teak flowering developmental stage.