Abstract: Background and Objective: Buasbuas is the one of the medicinal plants in Indonesia that contains bioactive compounds potential as antimicrobial, antioxidant, antidiabetes, antiinflammation and anticancer. Exploring the pathway and gene related of buasbuas bioactive compounds production has led to the renaissance of understanding buasbuas molecular mechanism database. The aim of this study was to developed data-mining framework of buasbuas to study plant specialized metabolism for phytochemical biosynthesis. Material and Methods: This project was started by collecting shoots and leaves of Buasbuas. Focus of the project was exploring the molecular mechanisms on biosynthesis phytochemical of Buasbuas. Illumina Mi-Seq Next Generation Sequencing was utilized to understand the molecular mechanisms of biosynthesis. Transcriptomes then trimmed and assembled with CLCBio genomic software. Assembled contigs then annotated towards Arabidopsis thaliana using CLC Bio genomic software. Digital Gene Expression was performed to analyze the transcriptional changes in control culture and treatment. Results: There were 5.342 unigenes that expressed only in treatment shoot cultures. Annotation with Gene Ontology showed that 57.9% (3.446) unigenes play role in Biological Process, 56.7% (3.375) unigenes play role in Cellular Components and 63.4% (3.772) unigenes play role in Molecular Functions. Annotation with Kyoto Encyclopedia of Genes and Genomes shows 853 unigenes essentially have role in 24 biological pathways. Highest process with highest unigenes involvement is biosynthesis of plant hormones and biosynthesis of alkaloids. Conclusion: This study showed that phytochemical biosynthesis in buasbuas induces level expression of several genes involved in the jasmonic acid, cytokinin, gibberellin, salicylic acid and ethylene biosynthesis pathway.
INTRODUCTION
Buasbuas (Premna pubescens. Blume.) is herbaceous and treedicots in the mint family, Lamiaceae. Buasbuas is native to Southeast Asia and could live in the forests of Sumatra and Malaya peninsula1. P. pubescens have been used in the traditional system of medicine in Malay ethnic and widespread Sumatra to treat various infectious diseases2. The P. pubescens are rich in the various types of phytochemicals such as alkaloids, flavonoids and terpenoids3. A previous research revealed that buasbuas leaves extract contained potential anti-inflammation, antimicrobial, antioxidant, antidiarrheal, antidiabetes, antiinflammation, analgesic, wound healing and anticancer2-4.
The potential possessed by the P. pubescens is very potential but it is not yet accompanied by the presence of a sufficient molecular database. Until now, Genomic analysis of P. pubescens, a non-model medicinal plant, is limited by the small quantity of publicly available sequence data. However, the emergence of next generation sequencing has paved the way for large scale sequencing of several non-model plants which can be valuable in investigating the basis of medicinal properties of such plants5,6. Next Generation Sequencing (NGS) technologies and their potential applications in plant biology including transcriptome investigations have been reviewed5,7. Strategies and tools which can be employed in transcriptome studies of non-model plants using second generation sequencing have beendiscussed6,8,9.
Non-model plants that have been recently sequenced include Olive10, Reaumuria trigyna (desert grass)11, Medicago truncatula12, Salicornia europea13,Gossypium arboreum14, Hevea brasiliensis15, Sesamum indicum16, Ipomoea batatas17, Camelliasinensis18, Acacia auriculiformis19, Acacia mangium19, Cajanus cajan20, Euphorbia fischeriana21, Myrica rubra22, Tectona grandis7-9 and many others are in progress. We have undertaken an NGS based approach to sequence the P. pubescens transcriptome in order to identify and characterize transcripts potentially contributing to the observed medicinal properties. This study will aid in the understanding of the bioactive phytochemical potential of P. pubescens and serve as a valuable resource for numerous researchers working on developing biosynthetic phytochemical. Availability of this transcriptomic data in public domains will also enable genome wide comparative studies of closely related medicinal plants importance.
We are interested in studying more about the roles of genes that control biosynthetic phytochemical pathway mechanism in P. pubescens especially in young and old leaves. This study is the first research on the database of gene sequences in P. pubescens plants. This study was preliminary of P. pubescens transcriptome characterization, before isolation and characterization of functional genes involved in biosynthetic phytochemical mechanism pathways. This study was conducted to understand phytochemical biosynthetic of buasbuas. Further analysis, total RNA of the leaves was extracted and transcriptome analysis was conducted7-9. In this report, the transcriptomes were deeply analyzed to understand the molecular mechanisms of phytochemical biosynthetic occurred in the leaves.
In this study, we sequenced the transcriptome of P. pubescens using the next generation of high throughput paired-end RNA sequencing (RNA-seq) technology, Illumina MiSeqTM 200010-12. Then, CLC bio bioinformatics technology tool was used to perform a de novo assembly and annotation without prior genome information13,14. This transcriptome database helped to reveal much about the functional genomics of P. pubescens and was then used to predict the functional classification of many unigenes using GO and KEGG pathway analysis15. 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 P. pubescens. Furthermore, we focused on the sequences that are related to phytochemical biosynthetic process in the aim of exploring the relationship between genes in phytochemical biosynthetic mechanism pathways. Understanding of P. pubescens phytochemical biosynthetic could be applied to improve P. pubescens production of bioactive compounds.
MATERIALS AND METHODS
Explant preparation: This research project was conducted from November, 2016-November, 2017 in Plant Science and Biotechnology Laboratory School of Life Sciences, Institute of Technology Bandung and Cell and Molecular Biology Laboratory of Medan State University. Young leaves and old leaves of P. Pubescens were obtained from the Laboratory of Biology Medan State University. Sample was in liquid nitrogen then sent to Plant Science and Biotechnology laboratory, School of Life Sciences, ITB.
RNA Isolation: Total RNA was isolated from Buasbuas. Isolation was done with CTAB-LiCl protocol15. Grinded tissues were submerged into CTAB buffer followed by extraction with 2x chloroform. RNAs were precipitated using LiCl and followed by second extraction using chloroform. Least precipitation were done using ethanol. The RNA pellet was suspended in 20-100 μL DEPC water8,15. RNA quality was evaluated using Agilent 6000 RNA Nano Chip Kit on 2100 Bioanalyzer (Agilent Technologies). Qualified RNA sample was sequenced at Macrogen Korea using Illumina MiSeqTM 2000 (Illumina Inc., San Diego, CA, USA) using protocol of the Paired-End sample Preparation kit12,16,17.
Transcriptome data analysis: All sequences were trimmed and assembled with CLCBio Bioinformatics software16. Trimming were done towards paired-end reads for quality score and the presence of repeated sequences >50 bp using CLCBio default parameters. Assembled contigs then annotated with nucleotide and protein sequences of Arabidopsis thaliana (www.phytozome.com). Both process and quality selection of the de novo assembly were done using CLCBio workbench16. Orthologous genes were retrieved and performed a local BLASTX alignment with e-value<10–6. Digital Expressed Gene (DEG) library was evaluated using the same programe. Reads Per Kilobase per million reads (RPKM) was used to analyze the differences of expression of each contig in old leaves and young leaves. Further analysis, contigs were annotated using Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 bioinformatics software at http://david.abcc.ncifcrf.gov/23,24. Contigs were aligned to the Gene Ontology (GO) and Kyoto Encyclopedia of Genes by and Genomes (KEGG) (http://www.genome.jp/kegg/)using the same programme.
RESULTS
Transcriptome statistic of buasbuas: Total transcriptomere adings on buasbuas leaves are 26.272.763 bases (Table 1). After assembling and trimming as much as 48.782 sequences contig obtained. Contigs then annotated to A. Thaliana genome resulting 36.549 contigs were annotated using BLASTn using CLCBio software and BLASTX using NCBI-nt database. Longest contig was assembled by 5.239 bp. The shortest contig was assembled by 225 bp. The mean of contig length was 536 bp. The proportion of contigs that have length more than 3.000 bp were 0,13% from all contigs, 67 contigs (Fig. 1). Figure 1 also informed that more than 24.000 of 48.782contigs were constructed from 200-400 bp. The data also reported that although overall contigs’ average length were 536 bp, the number was abundant in 201-2000 bp. These contigs then further analyzed with DEG to understand the different expression in control and salinity stress condition.
Species distribution of buasbuas transcript was highly varied hit plant species sequences (Fig. 2). It has highest similarity with Vitis vinivera (3.630 hits) (Fig. 2). It also has a high similarity with Oryza sativa species compared with other species distribution. Although the similarity with Vitis vinivera is higher than Musa spp. that show it share some homolog sequences between Musa spp. and Brachypodium distachyon. Increased expression that occurred in leaves further analyzed by comparing the Reads PerKbper Millionreads (RPKM) of both group. The grouping was applied on ortholog contigs that expressed in both group and annotated against A. thaliana. Based on these characteristics, we identified 11.338 contig experiencing up-regulation in old leaves. There are 24.009 unigenes expressed in same rate from both group. Some genes experiencing down-regulation (13.624 unigenes) and only expressed young leaves. Contigs are further analyzed using the DAVID (Database for Annotation, Visualization and Integrated Discovery, v. 6.7) bioinformatics software.
| Table 1: | Statistic of Buasbuas(P. pubescens) leaves Transcriptome |
| Fig. 1: | Distribution of contigs length of sequenced total RNA from buasbuas leaves |
| Fig. 2: | Species distribution of BLASTX annotation using NCBI-nt |
| Fig. 3: | Biological process Gene Ontology (GO) of buasbuas (P. pubescens) in young leaves and old leaves |
Gene ontology profile of the up-regulated contigs: Gene Ontology (GO) analysis indicates involvement of contigs in several physiological functions of the cells (Table 2). A total of 85.4% unigenes play role in the biological process (BP), 67.9% unigenes play role in cellular component (CC) and 84.7% unigenes play role in Molecular Function (MF) in stress culture. Compared to other group, treatment group show higher percentage of genes involved in three gene ontologies.
Different level of gene expression in three condition, old leaves, young leaves and both were showed in Fig. 3. Highest process that happened in old leaves were response to abiotic stimulus, followed by phosphorus metabolic process, phosphate metabolic process, phosphorylation and protein amino acid phosphorylation.
| Table 2: | Buasbuas Gene Ontology (GO) annotated to A. thaliana |
Those process also happened in young leaves with higher number of gene involved. Both leaves also result in same process that have higher genes number involved. There is one process shown only in young leaves that is regulation of transcription. Under biological process of old leaves we could know that the plant (old leaves) were responses some stress such as osmotic stress (2.2%), salt stress (2.1%) and endoplasmic reticulum stress (0.1%). This responses did not came out in young leaves.
| Fig. 4: | Molecular function Gene Ontology (GO) of buasbuas |
| Fig. 5: | Cellular component Gene Ontology (GO) of buasbuas |
One of the response to abiotic stimulus is response to stress (2.1%). Genes that are expressed in response to this stress are Osmotin-like protein (OSM34), Aquaporin (PIP1-2), Sodium/hydrogen exchanger 1 or Na+/H+ Exchanger (NHX), Sodium/hydrogen exchanger (NHX), Cysteine proteinase (RD19A), Heat shock protein 81-2 (HSP81-2) and Delta-1-pyrroline-5-carboxylate synthetase B (P5CSB). We seem to identify seven genes involved in stress response of Buasbuas. Those genes also widely known as response to stress such as cold, drought and salinity in many species. Expression of shared genes from both leaves have higher number of genes than genes expressed only in treatment or young leaves. Highest genes involvement in molecular function is nucleotide binding. Molecular function data also reveal that genes expression in young leaves have higher number specific genes than old leaves across all categories. Genes expression in metal ion binding, ion binding and cation binding not shown specifically in old leaves (Fig. 4).
Different from two ontologies earlier, cellular component of genes expressed in old leaves show three specific component, mitochondrion, non-membrane-bounded organelle and intracellular-non-membrane-bounded organelle. This component highly expressed only in old leaves. Component cell such as envelope and organelle envelope also highly expressed in old leaves than young leaves. Genes expression of plastid and chloroplast activity show similar expression level in both leaves (Fig. 5). Ten of the processes that occur in each ontology function shown in Fig. 5. Many cellular events occur in the explant associated with nucleotide binding to DNA. It is also supported by the low activity of the helicase on leaves. The cell indicates actively make compartment that shown by high activity of plasma membrane, plastid and chloroplast in CC. The result could implicated that the cell is doing signaling cascade highly in the old leaves.
Membrane related activity observed lower expression in old leaves. Similar result also shown in Arabidopsis grown in salinity medium and the activity induce senescence of the explant. It could be an indication that the leaves can tolerate the maturity process. The leaves had a high mitochondrion activity rather than plastid and chloroplast. Mitochondrion activity is recorded increase during salinity and drought stress in Arabidopsis related to proline catabolism. Proline is broken down into glutamate by catalytic enzyme such as P5CS. Another function could be explained by further analyzing from KEGG pathway happened in stress condition.
KEGG pathway of the unigenes: Up-regulated unigenesin old leaves then further analyzed using the DAVID program to understand biosynthetic pathways using Kyoto Encyclopedia of Genes by and Genomes (KEGG) annotation. Mapping results showed that there were24 KEGG pathways that occur in the explants. The highest observed physiological mechanisms are the biosynthesis of the plant hormones. Old leaves resulted in higher production of plant hormones. The subsequent highest mechanism that occurred is biosynthesis of secondary metabolites. Figure 3 illustrated that there was an increase in alkaloid biosynthesis in the leaves that produced from derivative of histidine and purine, sikim at acid pathway, derivative ornithine, lysine and nicotine acid, as well as from terpenoid and polyketide biosynthetic pathways. Secondary metabolites require a precursor from primary metabolites such as starch and fatty acids (Fig. 6). As a signal, many plant hormones are produced to initiate physiological process further defense against salinity. We could presume that those products were important in responding old leaves.
KEGG annotation result indicated that the unigenes activate multiple biosynthetic pathways groups. Despite plant hormones and alkaloids biosynthesis, the unigenes indicate increasing activity of spliceosome, glycolysis/gluconeogenesis, pyruvate metabolisms, endocytosis and arginine and proline metabolism. Non-homologous end-joining (NHEJ) also occurred in the explants during responding to salt stress. This result could inform that old stress also increase another metabolism pathway related to glucose metabolism and osmolytes biosynthesis such as proline.
Gene involved in biosynthesis of phytochemical: There are 199 unigenes that involved in biosynthesis of phytochemical from the all contigs analyzed. Most of genes were involve in glycolysis, citric aced cycle and PEP pathway. There are seven hormone biosynthesis pathways that regulated by the putative contigs such as cytokinin biosynthesis, carotenoid biosynthesis, gibberellin biosynthesis, brassinosteroid biosynthesis, jasmonic acid biosynthesis, salicylic acid biosynthesis and ethylene biosynthesis. Here we focused to observe genes involved in those biosynthesis pathways. A list of genes that involved in phytochemical biosynthesis pathway with comparation of their expression in stress and control condition was displayed in Table 3. There were three genes involved in cytokinin biosynthesis pathway which evince higher expression in control treatment. Three genes involved in carotenoid biosynthesis and all of the genes were chloroplastic genes. Those genes were responsible in producing cytokinin hormones that involved cell division both in root and shoot. The expression of those genes were higher in young leaves than old leaves. This output could informed that cytokinin expression was suppressed during maturity process.
Carotenoid biosynthesis happened with four genes involved, CRTISO (Carotenoid isomerase), CYP97B2 (Cytochrome P450 97B2), PDS3 (Phytoene desaturase) and PSY (Phytoene Synthase). Gibberellin biosynthesis only showed in stress condition and only one genes involved. The gene, with 9,907 RPKM value, was responsible for the 20ox1 (gibberellin 20 oxidase 1) protein. Higher expression in stress treatment also shown in brassinosteroid and jasmonic acid biosynthesis. Four genes involved in brassinosteroid biosynthesis with total RPKM value 141,552. The four genes are responsible for expression 7-Dehydrocholesterol reductase (DHCR7), sterol 14α-demethylase (CYP51G1), cycloartenol synthase (CAS) and probable squalene epoxidase 1 (XF1). Seven genes involved in jasmonic acid biosynthesis with total RPKM value in stress condition 265,789 are 3-ketoacyl-CoA thiolase 2 (KAT2), fatty acid beta-oxidation multifunctional protein (MFP), peroxisomal acyl-coenzyme A oxidase 1 (ACOX1), 12-Oxophytodienoate reductase 2 (OPR2), lipoxygenase 3 (LOX3), lipoxygenase 5 (LOX5) and Lipoxygenase 6 (LOX6). These genes can further inform bigger picture of plant response toward maturity process.
| Fig. 6: | There are 24 KEGG pathways of the genes involved in biosynthetic phytochemical of buasbuas |
The highest hormone production was showed by genes involved in biosynthesis of jasmonic acid (value = 173,52). This result followed by brassinosteroid biosynthesis pathway (value = 74,75). Figure 7 also showed that gibberellin acid only expressed in old leaves. Interestingly, we pictured lower expression of cytokinin and salicylic acid in old leaves. Thus, the old leaves assumed to increase carotenoid/ABA, gibberellin acid, brassinosteroid, jasmonic acid and ethylene.
DISCUSSION
Premna pubescens, also known as buasbuas is a tropical tree native from moist tropical forests of Asia family of lamiaceae1. Premna pubescens have been used in the traditional system of medicine in North Sumatra ethnic to treat various infectious diseases3. The other report indicated P. pubescens are as an antioxidant, anticancer, antiviral, antifungal and antibacterial2-4.
| Table 3: | Expressions of Genes involved in plant hormones biosynthesis pathway from KEGG pathway analysis |
| Fig. 7: | Expressions of genes involved in phytochemical biosynthesis pathway in leaves, CK: Cytokinin, CAR: Carotenoid, GA: Gibberellin acid, BR: Brassinosteroid, JA: Jasmonic acid, SA: Salicylic acid and ET: Ethylene |
Despite its importance, very little is information about the biology of this species at the genetic, molecular and biochemical levels. There is no biological sequences reported to date in public databases and belonging to buasbuas. This research is the first study of sequence gene databases in buasbuas.
Transcriptome statistic of buasbuas: 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 programs7-9. Next generation sequencing (NGS) technology during the last decade have dramatically impacted genome sequencing and transcriptome analysis7,15. Illumina’s MiSeq platform introduced recently constitutes a better alternative for transcriptomics21. 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 like P. pubescens. This approach has been used for analyzing the transcriptomes of P. pubescens. This study generated 48.782 high quality contigs from young and old leaves with average length 536 bp. Similar result also shown in Diningrat et al.7-9 studied Tectona grandis using Illumina’s MiSeq7-9. Transcriptome analysis with de novo method that worked by Diningrat et al.7-9 on Tectona grandis resulted 87,365 contigs with average length 486 bp7-9. A fraction of P. pubescens contigs could be annotated using the Arabidopsis proteome than Vitis vinivera, Oryza sativa, Musa sp. or Brachypodium distachyon not following the results of the BLASTX analysis as it relates to the availability of data in KEGG for biological processes referring to the arabidopsis database23-25. The results of the BLASTX analysis showed that genetically P. pubescensmore resemblance to the plants mentioned above26,27. By taking into consideration the sequences that have homologies in the Arabidopsis proteome, assuming that the two samples of P. pubescens have a similar gene number as Arabidopsis28,29. The cDNA sequences generated from both P. pubescens 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. Thus, this result expected to give deeper understanding of buasbuas phytochemical biosynthesis.
Gene ontology profile of the up-regulated contigs: Gene Ontology (GO) analysis on biological process of the leaves showed a lot of metabolism involved. Response towards abiotic stress was the highest. The NHX protein that expressed play role as K+ homeostasis, pH regulation, regulating process from vesicle trafficking, cell expansion to plant development30. The protein was usually used as stress tolerant determinant in plant. PIP1-2 was water channel and transport water across membrane28. RD19A known as drought response gene and could initiate bud dormancy29,31. P5CS gene had widely studied as proline biosynthesis gene to promote osmoregulation of the cell in response to drought, cold and salinity32-34. Old leaves of buasbuas was expected responding maturity process with altering ion homeostasis, maintaining relative water content, performing bud dormancy and synthesizing osmolytes such as proline.
KEGG pathway of the unigenes: High concentration of jasmonic acid (JA) in plant indicates that the plant is in the biosynthetic of phytochemical mechanism condition34,35. This hormone trigger growth inhibition, senescence, leaf abscission and tuber formation36. Meanwhile, brassinosteroid hormone in plant play role in stem elongation, seed germination, cell differentiation control and deetiolation27,37. Biosynthesis of Jasmonic Acid was the highest expression activity in the leaves of buasbuas. Jasmonic acid is one of stress signal that expressed in abiotic and biotic stress37,38. This signal will conduct senescence of the old leaves. This information could give us assumption that browning leaves is experiencing senescence as a result of jasmonic acid production.
Gene involved in biosynthesis of phytochemical: Level expression of carotenoid in old leaves show higher result than young leaves. This result informed that the old leaves synthesize more carotenoid. Carotenoid expression is related strongly with abscisic acid expression38,39. Old leaves will increase ABA as stress signal promoting cells to express reactive oxygen species. This response could inhibit plant growth and induce senescence40,41. Old leaves tend to produce carotenoid and resulting in decrease chlorophyll concentration. The results of this transcriptome analysis generally showed that the genetic activity associated with phytochemical production occurring in older leaves was more active than that in young leaves of buasbuas.
CONCLUSION
The transcriptome analysis provides a powerful method for analyzing the buasbuas and investigating the transcriptional phytochemical biosynthesis. The de novo assembled buasbuas transcriptome provides an important resource for future investigations about the buasbuas in phytochemical biosynthesis. This study show that phytochemical biosynthesis in buasbuas induces level expression of several genes involved in the jasmonic acid, cytokinin, gibberellin, salicylic acid and ethylene biosynthesis pathway.
SIGNIFICANCE STATEMENT
Premna pubescens is very potential medicinal plant but until now, genomic analysis of P. pubescens is limited by the small quantity of publicly available sequence data. NGS based approach to sequence the P. pubescens transcriptome in order to identify and characterize transcripts potentially contributing to the observed medicinal properties. This study will aid in the understanding of the bioactive phytochemical potential of P. pubescens and serve as a valuable resource for numerous researchers working on developing biosynthetic phytochemical. Availability of this transcriptomic data in public domains will also enable genome wide comparative studies of closely related medicinal plants importance. Understanding of P. pubescens phytochemical biosynthetic could be applied to improve P. pubescens production of bioactive compounds.
ACKNOWLEDGMENT
The authors would like to thank Ministry of research and higher education Republic Indonesia for providing Research Grants (Post-Doctoral Research Grant, grant No. 045A/UN33.8/ LL/2017). The authors would also like to thank Medan State University (UNIMED) and The School of Life Sciences and Technology, Bandung Institute of Technology for providing the necessary facilities.