Genome sequencing of herb Tulsi (Ocimum tenuiflorum) unravels key genes behind its strong medicinal properties

Krishna Tulsi, a member of Lamiaceae family, is a herb well known for its spiritual, religious and medicinal importance in India. The common name of this plant is ‘Tulsi’ (or ‘Tulasi’ or ‘Thulasi’) and is considered sacred by Hindus. We present the draft genome of Ocimum tenuiflurum L (subtype Krishna Tulsi) in this report. The paired-end and mate-pair sequence libraries were generated for the whole genome sequenced with the Illumina Hiseq 1000, resulting in an assembled genome of 374 Mb, with a genome coverage of 61 % (612 Mb estimated genome size). We have also studied transcriptomes (RNA-Seq) of two subtypes of O. tenuiflorum, Krishna and Rama Tulsi and report the relative expression of genes in both the varieties. The pathways leading to the production of medicinally-important specialized metabolites have been studied in detail, in relation to similar pathways in Arabidopsis thaliana and other plants. Expression levels of anthocyanin biosynthesis-related genes in leaf samples of Krishna Tulsi were observed to be relatively high, explaining the purple colouration of Krishna Tulsi leaves. The expression of six important genes identified from genome data were validated by performing q-RT-PCR in different tissues of five different species, which shows the high extent of urosolic acid-producing genes in young leaves of the Rama subtype. In addition, the presence of eugenol and ursolic acid, implied as potential drugs in the cure of many diseases including cancer was confirmed using mass spectrometry. The availability of the whole genome of O.tenuiflorum and our sequence analysis suggests that small amino acid changes at the functional sites of genes involved in metabolite synthesis pathways confer special medicinal properties to this herb.


Background
Plants of the genus Ocimum belong to the family Lamiaceae (Order Lamiales) and are widely distributed in the tropical, sub-tropical and warm temperate regions of the world [1]. These plants are known to produce essential oils comprising of a number of aromatic compounds and Tulsi is rightly known as the "Queen of Herbs" for this reason. In India, these plants are mostly grown at homes for worship and as offerings in temples. Among plants with medicinal value, those belonging to the genus Ocimum are very important aromatic herbs or shrubs.
The genus Ocimum is highly variable and possesses wide genetic diversity at intra and inter-species levels.  [2]. It is difficult to distinguish all these species on the basis of leaf morphology alone (Fig. 1). The metabolites (essential oils) of genus Ocimum have been reported to possess antioxidant and antifungal properties and to cure many diseases including bronchitis in Ayurveda, an Indian system of medicine [3]. Plants produce specialized metabolites as part of their defense mechanisms and these metabolites have significant medicinal properties that cure several human diseases. They can be isolated from various parts of the plant, including leaves, flowers, roots, bark, seeds and stem [4]. Pharmacological screening and the systematic study of the chemical constituents of plant metabolites provide a basis for developing new drugs. Some of the important metabolites reported from Ocimum species include linalool, linalyl, geraniol, citral, camphor, eugenol, methyleugenol, methyl chavicol, methyl cinnamate, thymol, safrol, taxol, urosolic acid etc. [4]. These metabolites are of immense value in the pharmaceutical, perfume and cosmetic industries. Metabolites derived from Ocimum species have been found to contain many medicinally relevant properties including anti-cancer, antioxidant, antifungal and anti-inflammatory virtues, and are also recommended for the treatment of malaria, bronchitis, diarrhea, dysentery, etc. [5]. Essential oils produced as specialized metabolites found in leaves, seeds, flowers and roots of Ocimum species are used in pharmaceutics and many systems of traditional Indian medicine [3,4]. Genome and transcriptome sequencing of medicinal plants serve as a robust tool for gene discovery and downstream biochemical pathway discovery of medicinally important metabolites [6]. Recently, an abundance of transcripts for biosynthesis of terpenoids in O. sanctum and of phenylpropanoids in O. basilicum [7] was reported during an attempt to compare transcriptomes of the two species of Ocimum. Despite its important role in traditional Indian medicine and its impressive arsenal of bioactive compounds, our understanding of Krishna Tulsi biology is limited. In this paper, we present the draft genome sequence of the non-model plant O. tenuiflorum (subtype Krishna), along with transcriptomes of two subtypes, Krishna and Rama Tulsi from leaf samples. We have identified a large set of genes involved in the production of specialized metabolites of medicinal interest such as apigenin, luteolin, rosmarinic acid pathway, eugenol, and ursolic acid.

Genome sequencing and assembly of the non-model plant O. tenuiflorum subtype Krishna
The paired-end (PE; 2x100-bp) and mate-paired (MP; 2x50-bp) DNA libraries were generated for Krishna Tulsi subtype using Illumina protocols. In total we obtained 373 million reads of PE and 166 million reads of MP data for Krishna Tulsi. Low quality (LQ) sequence reads were trimmed (Additional file 1: Figure S1 and Additional file 2: Figure S2) and reads with quality scores of less than Q30 were removed. The good quality reads were used for de-novo genome assembly. Median insert size of PE data was 335 (with a median absolute deviation of 21), whereas median insert size of MP data was 2473 (with a median absolute deviation of 704). K-mer 43 was opted as the best assembly from the statistical analysis of different k-mers. We obtained a maximum scaffold length of 184.7 Kb (Table 1) with an N50 length of 27.1 Kb. This assembly gives rise to a total of 78,224 scaffolds including equal to or more than 100 bp. The current draft assembly of Krishna Tulsi genome is 374.8 Mb in length. The genomic content of Krishna Tulsi is 0.72 pg/2C which is equivalent to 704.6 Mb [8], but the estimated genome size by k-mer method is 612 Mb and 61 % of the estimated genome size was assembled. The genome size reported in the literature [8], may be of a different cultivar. This lower genome coverage may be due to limited sequencing data (only two libraries were used in sequencing) or due to a high percentage of repeats (42.9 %). In terms of depth of sequencing, we sequenced 59× of the genome with paired-end (100 bp) and mate-pair (50 bp) libraries (since one lane can produce approximately 30Gb of data, even assuming that reads cover the entire 612 Mb of the estimated genome size). Ocimum species are characterized by the different basic chromosome numbers x = 8, 10, 12, or 16 [9,10]. In case of O. tenuiflorum individuals with 2n = 32, 2n = 36, and 2n = 76 have been recorded and the chromosome number of O. tenuiflorum is observed to be 2n = 36 [8].
A comparative analysis of the assemblies generated using PE data alone and with both PE and MP data show that the size and quality of the genome assembled using PE data alone improved substantially with the inclusion of MP data (Additional file 3: Figures S3 and Additional file 4: Figure S4, Additional file 5: Table S1 and Additional file 6: Table S2).
Validation of de novo genome assembly, annotation and repeat content of Ocimum tenuiflorum subtype Krishna genome The de novo genome assembly was validated by mapping raw reads to the assembled genome. On an average, 74 % of reads were mapped back to the assembled genome. Almost 83.3 % of the RNA-seq reads were mapped to the assembled genome. The completeness of de novo genome assembly and annotations were also checked with two other approaches i.e., by using CEGMA (Core Eukaryotic Genes Mapping approach) [11] and DEG (Database of Essential Genes) [12] (please see Methods for details). First, we searched for essential eukaryotic genes in the O. tenuiflorum assembly. This resulted in the mapping of 85.1 % of complete core proteins (CEGMA) and more than 95 % including partial genes against our genome assembly (Additional file 7: Table S3). Secondly, we searched for the predicted genes from the final assembly of essential genes recorded in the DEG database. We observed that about 89 % of essential genes were included within the assembly. These genes were also validated using Pfam domain annotation and were of comparable domain lengths as the classical members of that family (Additional file 8: Table S4). Phylogenetic trees for highly conserved essential genes like glyceraldehyde 3phosphate dehydrogenase (Additional file 9: Figure S5), cytochrome P450 (Additional file 10: Figure S6) and actin (Additional file 11: Figure S7) from Krishna Tulsi and their respective homologues were analyzed and compared with other plant species. Krishna Tulsi genes were found to cluster with genes belonging to related species namely, These trends further support the quality of the genome assembly.
Regarding the repeat content of the genome, we identified 78224 repeat regions, with a GC content of 36.1 %, adding to 160889218 bp (160 Mb), which constituted 42.9 % of assembled genome which is 374806882 bp (374 Mb) long (Additional file 12: Table S5). Long terminal repeats (LTRs) are found in large numbers in plant genomes (Schmidt T, 1999) and a similar trend is also found in the type of repeats identified in the Tulsi genome.

Genome annotation
We identified 36768 putative gene models in the initial genome draft (version 1.2) of O. tenuiflorum genome. At least one gene was observed in each of the 10012 scaffolds, with an average of three to four genes per scaffold. During the process of refined gene prediction, 16384 gene models were observed to have expression evidence (RNA-Seq data from leaves of Tulsi (Krishna and Rama)). A total of 19384 gene models have been identified by ab initio means (without any RNA or protein evidence) ( Table 2).
All the gene predictions, with or without RNA/protein evidences, were screened based on length (>100 bp). In case of sequential overlaps between different gene models, the gene models which are of longer length and with RNA or protein evidence for a given region of scaffold were preferred over the ones without any evidence.
There are 31,020 genes with at least one homologue in NRDB and 24,607 genes which contain at least one Pfam domain. In total, 3929 unique Pfam domains were identified for all the predicted genes in Tulsi (please see URL: http://caps.ncbs.res.in/Ote for the full list of predicted genes). A majority of domains identified were protein kinases or LRR-containing domains (Additional file 13: Figure S8). Further comparison of Pfam results, with assembled plant genomes of similar size, reveals that the number of predicted gene models is in overall agreement in numbers as well as gene boundaries.

Orthology of Tulsi genes
The orthology relationships were deduced between Krishna Tulsi (O. tenuiflorum; Ote) and four other species viz. Arabidopsis thaliana (Ath), Mimulus guttatus (Mgu), Solanum lycopersicum (Sly) and Oryza sativa (Osa) (please see Methods for details). We observe 8370 clusters which contain a total of 89922 gene products from the five plant species (Fig. 2a). M. guttatus and O. tenuiflorum share the same order (Lamiales), but belong to different families (Phrymaceae and Lamiaceae, respectively), which was evident from the presence of the highest number of common gene families (11707) between them. This was followed by Solanum lycopersicum (11022), Arabidopsis thaliana (10206) and Oryza sativa (9154) as expected from taxonomic hierarchy (Fig. 2a). We found 17584 genes to be orthologous to any of the above four species. Considering all the 36768 Ote genes, 1282 groups contained only Ote Krishna Tulsi genes (3302). We obtained 16 Ote genes which lack traceable orthology to 22 other plant species and homology relationships (list of these genes is available on the database). Few of these unique Ote genes are transposons.
In order to inspect in detail the distribution of the orthologous relationship of Ocimum genes across different species and taxonomic levels, 22 fully-sequenced plant genomes (Additional file 14: Table S6) were considered. The orthologous groups from all 23 species were organized according to the clustering. Three hundred and thirty four clusters of genes are present across all the 23 species chosen for the study. Common genes across all species, comprising of their respective orthologous group, are plotted as a horizontal stacked bar plot (Fig. 2b). The pattern of sharing orthologous groups is quite unique to primitive plant genomes (like lycophyte and bryophyte) and monocots. However, the pattern observed in the Tulsi genome is quite similar to that of M. guttatus (Mgu). Interestingly, this pattern is somewhat different for two members of Solanacea, which have more genes shared only in two out of 23 genomes, perhaps due to other features such as polyploidy.
Genes involved in synthesis of specialized metabolites of medicinal value: comparative analysis between O. tenuiflorum (Ote, Krishna Tulsi) and other plant genomes Next, we performed a restricted analysis of the genes involved in the metabolite production in Ote and the genomes of a few plant species that are either closely-related (S. lycopersicum, V. vinifera) or well-characterised (M. truncatula, and A. thaliana). We observed 121 (72.45 %), 130 (77.84 %), 106 (63.47 %) and 94 (56.28 %) scaffolds and contigs from the selected four representative genomes associated with 167 metabolite-related scaffolds and contigs in Ote Krishna Tulsi (Fig. 3) respectively. In terms of the number of orthologous genes from this selected plant genome associated with metabolite genes of Ote, we observed a similar trend of association as 601, 620, 570 and 556 genes in S. lycopersicum, V. vinifera, M. truncatula, and A. thaliana respectively. These numbers agree with the taxonomic phylogeny and hierarchy, suggesting that the evolution of genes involved in metabolic pathways is not a cause of recent expansions or sudden drifts. When compared against 11,389 scaffolds (greater than 10Kb in size) from Ote, 10032, 9997, 8648 and 8277 scaffolds were found to be associated with the four reference plant genomes (Additional file 15: Figure S9, Additional file 16: Figure S10 and Additional file 17: Figure S11 for three genomes and Additional file 18: Table S7 for four genomes). Further, most of the metabolite-related scaffolds in Ote Krishna Tulsi were associated with chromosomes 1, 6, 8, and 10 of tomato (Fig. 4). In particular, gene products that are likely associated in luteolin synthesis pathway are observed to cluster in scaffolds, which are similar to nucleotide stretches in Chromosomes 3, 5, 6, 8 and 10 of the tomato genome (Fig. 4).
Transcriptome de novo assembly of Krishna and Rama Tulsi mature leaf samples De novo transcriptome assembly was performed for the mature leaf samples of subtype Krishna Tulsi. The best assembly resulted in 109291 contigs with N50 of 893 bp and longest sequence of 12.1 Kb. All these contigs added up to 49.5 Mb with a GC content of 42.9 %. Scaffolding of these contigs resulted in 89878 scaffolds with N50 of 1597 bp and longest sequence of 12.7 Kb. All these scaffolds added up to 56.3 Mb with a GC content of 42.9 % (Table 3). Similarly, assembly was performed for the subtype Rama Tulsi and combined reads (Krishna and Rama Tulsi) as well (Table 3).

Differential expression of transcripts
The differentially expressed genes found in the transcriptomes of both the Tulsi subtypes were analysed. We observe a substantial number of genes up-regulated and down-regulated in Krishna Tulsi, compared to Rama Tulsi. Some of the highly expressed genes were also confirmed by q-RT-PCR technique in different tissue samples i.e., stems, leaves and flowers and also in five species viz. O. tenuiflorum subtype Krishna and Rama, O. gratissimum, O. basilicum, and O. kilmand. The numbers indicate percentage of association of these genomes with the metabolite genes of Ocimum genome. These percentages agree with the taxonomic phylogeny and hierarchy, suggesting that the evolution of genes involved in metabolic pathways is not a cause of recent expansions or sudden genome drifts. The inner circle represents chromosomes from respective homolog genome. Each scaffold is organized in the middle circle and is represented in chronological order as per position on chromosomes. The line represents location of each scaffold on the respective chromosome. Colors indicate = < 2 genes, =2 genes, = > 2 genes, = Metabolite related genes. Height of orange columns in outermost circle represents amount of repeats in corresponding scaffolds For a comparison, we generated a heat map of the top 50 differentially more abundant genes in Krishna Tulsi samples (Fig. 5a). Similarly, top 50 differentially more abundant genes in Rama with respect to Krishna sample were also plotted ( Fig. 5b). Gamma-cadinene synthase is one of the top 50 differentially expressed transcripts with RPKM values of 577.0 and 31.7 in Krishna and Rama Tulsi samples, respectively (please see below for details). Other highly expressed transcripts in Krishna Tulsi sample are Heat shock cognate protein 80, Cellulose synthase A catalyic subunit 6 (UDP-forming), Fructose-biphosphate aldolase (chloroplatic), Phototropin-2, and Rubisco activase 1 (chloroplatic). The chalcone synthase or naringenin-chalcone synthase (CHS) is one of the enzymes important for coloration of plant parts, which is observed to be highly expressed. Abundance values of all the transcripts, along with their functional annotations by NCBI BLAST results and their corresponding Krishna Tulsi genomic scaffold, show several genes involved in the synthesis of specialized metabolites implicated to be of medicinal value (Additional file 19: Table S8).  Dark purple coloration of the leaves and stem of subtype Krishna Tulsi is one of its characteristic phenotypes, which distinguishes it from other subtypes and species of genus Ocimum. Chalcone synthase (CHS) is an enzyme belonging to a family of polyketide synthases which catalyzes the initial step for flavonoid biosynthesis. Flavonoids are important plant specific metabolites that perform various functions such as pigmentation, antifungal defense etc. Reviewed protein sequence for CHS from UniProt (Universal Protein resource) database [13] was employed to search against annotated protein sequences of Krishna Tulsi genome and six transcripts were obtained as possible hits. The best hit could be identified with 95 % query coverage and 99 % sequence identity. The extent of abundance of this hit (protein sequence) was checked in the leaf transcriptome of both the Tulsi subtypes viz. Krishna and Rama. Abundance (in terms of RPKM) of the six transcripts was, on an average, two times more in case of Krishna as compared to Rama (please see Fig. 5), and may be involved in the coloration phenotype of Krishna subtype plants [14]. For further confirmation of expression of these transcripts, q-RT-PCR was performed. As expected, anthocyanin producing gene was observed to be more abundant in Krishna young leaf samples and mature leaf samples (used as control) ( Fig. 6a and b). In contrast, the chlorophyll binding protein was more abundant in Krishna mature leaf samples. In addition, we also examined the presence of gamma-cadeninene synthase gene which is responsible for aroma [15]. This gene was found to be more abundant in Rama root sample and young leaf samples of O. Saccharum, but not observed in higher quantities in O. kilmund.

Specialized metabolites detection and validation
Nearly 30 specialized metabolites (Fig. 7a) are reported form the genus Ocimum which are found to have medicinal values or properties [4]. Amongst these, 14 metabolites belonging to five basic groups were found to have complete pathway information in PlantCyc database (http://www.plantcyc.org/) [16] (Additional file 20: Figure  S12). Hence, genes involved in these pathways were chosen for further analysis and searched against the assembled genome of O. tenuiflorum. Figure 7b highlights the distribution of the genes identified in various classes of metabolites of disease relevance (i.e., these metabolites are well-known as drugs in the cure of human diseases).
A total of 458 genes were identified in Ote genome, which are either homologous or directly code for enzymes involved in the synthesis of specialized metabolites (Fig. 8) (details of gene IDs of these proteins are provided in Table 4 and Additional file 21: Table S9). Twenty eight O. tenuiflorum gene products were annotated as putative terpene synthases using BLAST sequence searches with E-value of 10 −4 and query coverage filter of >75 % (Additional file 22: Table S10).
Among these specialized metabolites, we focused on ursolic acid, belonging to sesquiterpenes, since it is known to have anti-inflammatory, anti-microbial, anti-tumour and anti-cancer properties. The synthesis of ursolic acid from squalene is a three-step process starting from squalene (Fig. 9). α-Amyrin is formed by concerted cyclization of squalene epoxide, while ursolic acid is eventually synthesized by the catalytic activity of multifunctional cytochrome P450. The enzymes involved are, therefore, squalene epoxidase, alpha-amyrin synthase and alpha-amyrin 2, 8 monoxygenase. Sequence search algorithms were employed to search for the three enzymes of this pathway in the Tulsi genome, starting from protein sequences for each of these enzymes from PlantCyc database as queries. The search for squalene epoxidase in Tulsi, using the sequence of this enzyme in Oryza sativa japonica (LOC_Os02g04710.2) as a query, gave rise to a hit (C3776143), with 50 % sequence identity covering 80 % of the query length (Additional file 23: Figure S13). Using Amyrin synthase LUP2 from A. thaliana (Q8RWT0) and 13 other well-accepted alpha/beta amyrin synthases as a query, four hits were identified in the Tulsi genome (scaffold16333, scaffold20801, scaffold12312 and maker-C3776143). In classical amyrin synthases, a QW structural motif repeats six times in the entire sequence [17,18], while there are two functional motifs, viz., a well conserved SDTAE [19] motif which is believed to form the catalytic pocket and the MWCYCR [20] motif that is shown to play Fig. 6 Expression quantification of selected genes by q-RT-PCR method. a. Fold changes of genes involved in color production, obtained through q-RT PCR. Blue colour horizontal bar is for chlorophyll a-b binding protein, red to denote Gamma-cadenine synthase and green to denote Anthocyanin. Mature leaf of Krishna subtype was used as control. It can be seen that, genes responsible for color production such as Chlorophyll a-b binding protein and gene in anthacyanin pathway are down-regulated as compared to mature Krishna leaf, which corresponds to phenotypic characteristics. b. Fold changes of genes involved in ursolic acid biosynthetic pathway, as obtained through qRT-PCR for 5 different Tulsi subtypes. Blue colour horizontal bar is for squalene epoxidase, red to denote alpha-amyrin synthase and green to denote Cytochrome P450 monooxygenase.  Further, a phylogenetic tree was constructed using 16 query sequences and these four hits (Fig. 10). One of the Tulsi hits, (scaffold 16333_mrnal) clusters with a well-characterized alpha amyrin synthase from C. roseus (H2ER439) suggesting that this particular scaffold might indeed retain an alpha amyrin synthase.
Interestingly, many genes involved in the synthesis of specialized metabolites of relevance in the treatment of diseases are also more abundant, as observed in the assembled transcriptome (Additional file 21: Table S9). Similarly, genes involved in the synthesis of 16 other specialized metabolites (Additional file 25: Table S11), are also equally interesting. However, this requires detailed understanding of the mechanism of synthesis and enzymes involved in the pathways. We analysed RNA-Seq data of two leaf samples in order to compare the genes related to important metabolite pathways and the peculiar phenotype of O. tenuiflorum subtype  In the case of the multifunctional Cytochome P450 (which catalyses the last three steps in the synthesis of urosolic acid, Fig. 9), a predicted gene from scaf-fold2032 was obtained as a hit, when a reviewed Uni-Prot entry F1T282 from V. vinifera was considered as query and searched in the Tulsi genome assembly using BLAST. This hit retains 61 % sequence identity and the alignment covers 90 % of the length of the query (alignments are shown in Additional file 23: Figure S13). This scaffold contains a total of three predicted genes viz., Ote100020320011, Ote100020320001 (similar to UHRF1-binding protein) and Ote100020 320031 (gene of interest).
From the available transcriptome assembly, these genes, identified as involved in the synthesis of urosolic acid, were analysed for their levels of expression. The RPKM values were also high for these three genes (please see Additional file 21: Table S9). To further validate the levels of expression of these genes, q-RT-PCR was performed using sequence-specific primers. The presence of these three enzymes is generally high in all the mature leaf samples and highest in Rama subtype (using Krishna subtype Table 4 The enzymes involved in metabolite biosynthesis were identified in the assembled genome and these genes were analyzed for their expression level in the transcriptome. The RKPM value signifies the level for expression (Continued)   Table 5).

Discussion
O. tenuiflorum subtype Krishna Tulsi is one of the nonmodel plants of great medicinal value, for which there has been no genomic information available till date. We have performed genome sequencing of O. tenuiflorum subtype Krishna of the paired-end (PE; 2x100-bp) and mate-paired (MP; 2x50-bp) DNA libraries by Illumina Hiseq 1000. The best de novo assembly was obtained at k-mer 43 by SOAP-denovo2, an eukaryotic de novo genome assembler. Repeats were identified and masked, and gene prediction and annotation was carried out using the MAKER annotation pipeline by using genomic, transcriptomics and EST data. The nearest species whose genome has been sequenced is the monkey flower (M. guttatus), which shares its order Lamiales with O. tenuiflorum   to as paralogs by OrthoMCL. Of the remaining Ote genes, 17,584 genes were found to be orthologous to any of the other four species studied in this case. We performed an analysis of the genes involved in the metabolite production in Ote and the genomes of a few other related plant species. Based on the direct evidence or homology a total of 458 genes were identified in Ote genome, which are involved in coding of enzymes implied in the synthesis of specialized metabolites. Comparative analysis of transciptomes of O. tenuiflorum subtype Krishna and Rama was performed to detect potential differentially-regulated genes and their involvement in metabolite synthesis. On comparing both the transcriptomes, differentially expressed genes were observed with a substantial number of genes more abundant and others less abundant in either subtypes. Gamma-cadinene synthase is more abundant in Krishna sample (RPKM value 577.047) as compared to Rama sample (RPKM value 31.73). To confirm some of the more abundant genes along with Gamma-cadinene synthase, we performed q-RT-PCR in different tissue samples i.e., stem and leaves and also in five species viz. O. tenuiflorum subtype Krishna and Rama, O. gratissimum, O. basilicum, and O. kilmand. Expression of Gamma-cadinene synthase is found more in Krishna samples as compared to Rama by q-RT-PCR also. Likewise, Chalcone synthase (CHS) is an anthocyanin-producing gene, which is observed to be more abundant in young leaf samples of Krishna and mature leaf samples in transcriptome data. Subsequently, this has been confirmed by q-RT-PCR and from mass spectrometry readings of ursolic acid and eugenol from different tissue samples and from different species.

Conclusion
We present a draft genome of O. tenuiflorum Krishna Tulsi subtype Krishna Tulsi. The habitat of genus Ocimum is tropical climate and it is wide spread over Asia, Africa, Central and South America. High RNA-seq expression values of the genes responsible for the purple coloration of the plant parts in Krishna subtype, as compared to Rama subtype, were observed. We also identified a fFew unique genes (16) of Ote, which lack any traceable orthology and homology relationships from all the 22 species used in this study.
Krishna Tulsi is described in the Vedas and Puranas (ancient scriptures of Hindus) and has a long history of cultivation, of roughly 3000 years, and is therefore assumed to be of Indian origin [21]. In literature, it is also referred to as the "Queen of Herbs". Major genes involved in the synthesis of medicinally important specialized metabolites in the plant could be unraveled despite limited data on sequencing and coverage [22]. Expressions of these genes were confirmed by complementing with RNA-seq data and q-RT-PCR method. We also investigated one of the important metabolic pathways involving the production of ursolic acid in detail, by massspectrometry and q-RT-PCR methods. Synthesis of specialized metabolites or their precursors appear to begin in the young leaves of Tulsi. Subsequently, the mature leaves retain the medicinally relevant metabolites. O. tenuiflorum Rama subtype retains the high abundance of key medicinally relevant metabolites like eugenol and ursolic acid, as observed in the transcriptome, metabolite quantifications and q-RT-PCR expression values consistent with its high medicinal values. Our main emphasis was to unravel the important metabolite genes by using genomic and transcriptomic data despite limited sequencing information.

Isolation of genomic DNA from O. tenuiflorum subtype Krishna Tulsi
Young leaves of Tulsi subtype Krishna and Rama were used for genomic DNA isolation. About one gram of leaves were crushed using liquid Nitrogen and DNA extraction buffer (200 mM TrisHCL [pH-8.0], 200 mM NaCl, 25 mM EDTA and 1 % PVP) was added [23]. The ground material along with 1/10th volume of 20 % SDS solution was incubated at 65°C for 30 min. The tubes were centrifuged at 14,000 RPM for 10 min at room temperature to remove the debris. The supernatant was transferred into a fresh tube and treated with equal volume of phenol: chloroform: isoamyl alcohol (25:24:1) and mixed gently for 5 min. The mixture was centrifuged at 12,000 RPM for 10 min to separate the phases. The aqueous phase from the centrifuged tube was transferred to a fresh tube and DNA was precipitated with 1/5th volume of 2 M NaCl and 2 volumes of ice-cold ethanol. The DNA was pelleted by centrifugation at 12,000 RPM for 10 min. Precipitated DNA pellet was taken as a starting material for purification using the Sigma Genelute plant DNA isolation kit (G2N70, Sigma). The DNA was run on a 1 % agarose gel to assess the quality. The A260/280 ratio and quantity were determined using the nanodrop.

Genome sequencing, assembly and annotation
Genome sequencing was performed by using Illumina HiSeq 1000 technology in the Next Generation Genomics Facility at Centre for Cellular and Molecular Platforms (C-CAMP). Genomic DNA paired-end and gel free mate-pair library preparation was performed for Krishna Tulsi using TruSeq DNA sample preparation kit (FC-121-2001) and Nextera mate-pair sample preparation kit (FC-132-1001) from Illumina (www.illumina.com). FASTX-Toolkit [24] and FastQC tools [25] were used for pre-processing of raw reads and for quality check of the reads. Genome assembly from reads of PE and MP together was done by using SOAPdenovo2, a de novo draft genome assembler [26]. Preliminary assemblies were performed based on k-mers from 21 to 63 with an interval of two. Gene prediction and annotation was carried out using the MAKER annotation pipeline [27] with predicted gene models using AUGUSTUS [28] and A. thaliana genes as reference for initial prediction. The gene models were refined using homology searches against all protein sequences from Viridaeplantae kingdom.

Validation of genome assembly and annotations
To validate genome assembly, we have mapped raw reads on to the de novo assembled genome by using REAPR (SMALT) [29], SAMtools [30] and Picard tools (http://broadinstitute.github.io/picard/). Maximum and minimum insert size of 500 bp and 0 bp respectively were selected for mapping. We report an alignment pairing with best score, using standard Smith-Waterman scores. The threshold minimum score used was calculated by the formula to be: < Minimum score > = < world length > + step size -1. Here the word length of 13 is used with a step size of 6. Estimation of the genome size of the Tulsi genome was done using the k-mer distribution analysis by Jellyfish [31]. Essential genes implicated in the regulation, assembly and functioning of plant cells, have been identified in the Krishna Tulsi assembled genome using a two-way approach. Firstly, using CEGMA which was derived from the KOG database [32] (for eukaryotic genomes) and core proteins in any eukaryotic genome (including ones in draft stages), essential genes were annotated. Secondly, a subset of A. thaliana genes were extracted from a well-characterized Database of Essential Genes (DEG) and compared against Krishna Tulsi assemblies. Validation of the extracted genes was performed by Pfam domain annotation approaches. Putative essential genes from the Krishna Tulsi dataset were further searched using BLASTP [33] against the NCBI (NR) database and closely-related homologues were aligned and phylogenetic tree constructed.

Repeat identification
Repeat elements in the assembled genome were identified using RepeatScout (version 1.0.5) [34] and Repeat-Masker (version 4.0.3) [35]. The library of ab initio repeats generated by RepeatScout was classified into known repeat classes using the RepeatClassifier module of RepeatScout (Additional file 12: Table S5). The RepBase library of RepeatMasker and the non-redundant library of ab-initio classified repeats were then used to mask the repeat elements in the assembled genome. The repeatmasked genome assembly was then used for genome annotation.

Genome annotation
The repeat-masked assembled genome of Krishna Tulsi was processed through the MAKER annotation pipeline [27]. AUGUSTUS [28] was used for gene prediction, trained on A. thaliana gene models. RNA-seq data obtained from leaf samples was used as EST evidence to refine the gene models. Initial gene models of protein sequences belonging to Viridaeplantae kingdom, obtained from the NCBI database, were used as protein evidence for refining gene prediction. Both EST and protein evidence were prepared using EXONERATE [36] and used for gene prediction refinement through AUGUSTUS. All the protein sequences of these gene models were subjected to validation based on identification of homologues through BLASTP search against NRDB at E-value cutoff of 10 −3 . Pfam release 27 was consulted for all domain predictions with an E-value cutoff of 10 −5 using HMMER3 package [37].

Orthology detection
All the predicted gene models from Krishna Tulsi were used with OrthoMCL tool [38] [40]. Phylogenetic tree reconstruction was carried out using 'RbcS' (Rubisco small subunit) coding sequences from all 23 species. CLUSTALW [41] and Phylip package [42] were employed for multiple sequence alignment (MSA) and subsequent clustering using Neighbor Joining (NJ) method, respectively. Distant homology relationships were verified through PSI-BLAST [33] at different set of E-value cutoffs. Gene products for which we were unable to establish any homology or orthology relationships, but consisted of a Pfam domain, were referred to as unique genes specific to Ote.

Comparative analysis between Krishna Tulsi and other plant genomes
The most recent version of whole genome sequences of S.

Transcriptome sequencing and assembly
We assembled all the mRNA reads having HQ scores of all the bases more than 20, of Krishna and Rama subtype separately and also by combining the reads from both of these subtypes by using SOAPdenovo-trans [26] at different K-mers starting from 19 to 63 at an interval of two. An insert size of 350 was used for the assembly of transcriptomes. RNA-seq reads were mapped to the assembled genome by Tophat2 [44], which uses Bowtie2 [45] as a mapping tool. We used a minimum and maximum intron length of 50 and 500000 bp respectively. Maximum multi hits (parameter that dictates the number of alignments to the reference for a given read) was assigned as 20 and transcriptome max hits (maximum number of mappings allowed for a read, when aligned to the transcriptome) of 60 was used.

Transcript differential expression comparison
To quantify expression in terms of reads per kilo base per million (RPKM), non-redundant combined assembled transcript sequences (at 90 % sequence similarity by CDhit EST [46]) were taken as reference. This non-redundant transcriptome was used as the reference transcriptome to calculate differential expression of transcripts in both the samples [6,47]. The reads of RNA-seq experiments from Krishna and Rama subtypes were mapped back on to the reference transcriptome by using SeqMap (version -1.0.12) [48] and RPKM values were determined by using rSeq: RNA-seq analyzer (version 0.1.1) [49].

Specialized metabolites detection and validation
The dataset obtained after gene prediction on the assembled genome was employed to search for enzymes involved in secondary metabolite production. There are 14 metabolites (flavonoids (2), phenylpropanoids (4), terpenes (2), sesquiterpenes (5) and sterols (1)), which are reported to be present in Ocimum and have known pathway information in PlantCyc (http://www.plantcyc.org/) [16]. Reviewed entries from the UniProt database and all the known sequences of the enzymes from other species possessing these enzymes were used as queries to search in the full dataset of scaffolds and contigs, using PSI-BLAST at E-value of 10 −5 and three iterations. The protein hits obtained in our dataset were further subjected to validation using a query coverage filter of 75 %. In order to study the expression of genes involved in the synthesis of specialized metabolite (s), the assembled transcriptome of both Ocimum species were searched, employing the reviewed entry corresponding to each enzyme in the UniProt database. These searches were performed using TBLASTN at an E-value of 10 −3 , and the best hit in our dataset was selected based on the least E-value. If the reviewed entry for any of the enzyme was not present, unreviewed entries from PlantCyc database were employed.
Quantification of eugenol and ursolic acid using UHPLC-MS/SRM method A Vantage TSQ triple stage quadrupole mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA) equipped with a heated electro spray ionization (HESI) source was used for the analysis of eugenol and an APCI probe was used for the ursolic acid analysis. The mass spectrometer was interfaced with an Agilent 1290 infinity UHPLC system (Agilent Technologies India Pvt. Ltd., India) equipped with a column oven (set at 40°C), auto sampler and a thermo-controller (set at 4°C). The needle was washed from outside with acetonitrile (0.1 % formic acid) before every injection to avoid any potential carry-over problems. Separations were performed using a shim-pack XR-ODSIII column (2 × 150 mm, 2 μm). For Eugenol: Mobile phase A was water (10 mM Ammonium acetate) containing 0.1 % formic acid, and mobile phase B was acetonitrile containing 0.1 % formic acid. For Ursolic acid: Mobile phase A was water (10 mM Ammonium acetate), and mobile phase B was acetonitrile: methanol (3:1). Injections of 10 μL were performed using flow through a needle (A)Eugenol: Eugenol was quantified after derivatizing with pyridine sulfonyl chloride using estrone-d4 as an internal standard. Methanol was used to extract eugenol from fresh leaves (2 mg/mL) and dried stem powder (20 mg/ml). Briefly 10 μL of extract and 10 μL of internal standard (from 2.5 μg/mL) were added into 200 μL of buffer [acetone: NaHCO3 (1:1)]. To this 10 μL of pyridine sulfonyl chloride (10 mg/mL) was added and incubated at 60°C for 15 min. After incubation the derivative was extracted with 800 μL of MTBE and the organic layer was dried and reconstituted in 50 μL of methanol followed by 10 μL injection for the analysis. A gradient (0-2 mins:30 %B, 2-5 mins:30-90 %B, 5-7 mins:90-100 %B, 7-10 mins:100 %B, 10-10.1 mins:100-30 %B, 10.1-15 mins:30) was then initiated at a flow rate of 200 μL/min. Operating conditions were as follows: spray voltage, 3000 V; ion transfer capillary temperature, 270°C; source temperature 100°C; sheath gas 20, auxiliary gas 5 (arbitrary units); collision gas, argon; S-lens voltage was optimized for individual metabolites; scan time of 50 millisec/transition; and ion polarity positive. A standard curve was constructed from 0.078 to 5ngon column to quantify eugenol. The SRM transition used for the analysis of eugenol is (306.1 → 79) and for estrone-d4 (416.3 → 274.1).

(B)Ursolic Acid:
Ursolic acid was quantified using estrone-d4 as an internal standard. A brief extraction was done from 2 mg/ mL of dry powder using 1 mL of methanol (sonication-3 min, centrifugation −5 min). The extract was further diluted to 0.2 mg/mL in methanol. From this extract 10 μL was added along with 10 μL of internal standard (0.1 ug/mL) to 30 μL of methanol and 10 μL was injected for the analysis. A gradient (0-2 mins:20 %B, 2-8 mins:20-100 %B, 8-14.5 mins:100 %B, 14.5-14.6 mins:100-20 %B, 14.6-20 mins:20 %B) was then initiated at a flow rate of 200 μL/min. Operating conditions were as follows: Discharge current 4 μA; ion transfer capillary temperature, 270°C; source temperature 300°C; sheath gas 20, auxiliary gas 5 (arbitrary units); collision gas, argon; S-lens voltage was optimized for individual metabolites; scan time of 50 millisec/transition; and ion polarity positive. A standard curve was constructed from 0.034 to 2.5 ng on column to quantify ursolic acid. The same standard curve was used for the analysis of oleanolic acid. The SRM transition used for the analysis of both ursolic and oleanolic acid is (439.4 → 119) and for estrone-d4 (275.3 → 257.1).