Small RNA profiling for identification of miRNAs involved in regulation of saponins biosynthesis in Chlorophytum borivilianum

Background MicroRNAs act as molecular regulator of cell signaling, plant growth and development, and regulate various primary and secondary plant metabolic processes. In the present study, deep sequencing of small RNAs was carried out to identify known and novel miRNAs from pharmaceutically important plant, Chlorophytum borivilianum. Results Total 442 known miRNAs and 5 novel miRNAs were identified from young leaf small RNA library. Experimental validation with stem loop RT-PCR confirmed the in silico identification. Based on transcriptome data of root and leaf of C. borivilianum, Oryza sativa, and Arabidopsis thaliana target gene prediction was done using psRNAtarget and mirRanda. BLAST2GO helped in localization of predicted targets and KEGG (Kyoto Encyclopedia for Genes and Genomes) pathway analysis concluded that miR9662, miR894, miR172, and miR166 might be involved in regulating saponin biosynthetic pathway. The correlation between miRNA and its target gene was further validated by RT-qPCR analysis. Conclusion This study provides first elaborated glimpse of miRNA pool of C. borivilianum, which can help to understand the miRNA dependent regulation of saponin biosynthesis and to design further metabolic engineering experiment to enhance their contents in the plant. Electronic supplementary material The online version of this article (10.1186/s12870-017-1214-0) contains supplementary material, which is available to authorized users.

Due to pharmaceutical importance of these phytoconstituents, research is now focusing on to enhance the saponins content in C. borivilianum. Increased saponin contents were reported when plantlets were inoculated with mycorrhizal fungi [20]. However, the molecular factors regulating the production of saponins contents are not identified yet. Since the discovery of miRNAs, these are known to regulate gene expression at transcriptional and post transcriptional level. In recent years, miRNAs were found to be major biological factor regulating the secondary metabolite production in many plants [21]. In Mentha spp., miR5021 was reported to check the essential oil biosynthesis by regulating expression of genes coding for enzymes such as geranyl di-phosphate synthase involved in 2-C-methyl-D-erythritol 4-phosphate/1-deoxy-D-xylulose 5-phosphate (DOXP) pathway [22]. miR5021 and miR5293 were found to regulate first enzymatic function of the MVA pathway in Panax notoginseng [23]. Based on the above studies, we hypothesized that there must be some miRNAs regulating the secondary metabolic pathway in C. borivilianum. Our aim was to profile miRNAs in C. borivilianum and to predict targets of the identified miRNAs. Further involvement of miRNAs targets in biochemical pathway were studied, with special emphasis on saponins biosynthetic pathway. To achieve this, we prepared small RNA library and sequenced it using Illumina platform. Known and novel miRNAs were predicted along with their gene targets. RT-qPCR analysis was performed to establish a correlation between miRNAs and their corresponding targets predicted by computational analysis.

Plant material
C. borivilianum plants were raised from vegetative buds of old plants in the soil:soilrite (1:1.5) mix under regulated environment in plant growth chamber with 27°C day/night temperature, illuminated with white light (flux density 200 μmolm −2 s −1 ) at Panjab University, Chandigarh (India) [latitude: 30°44′14 N; longitude 76°47′14E; altitude 350 m above mean sea level]. Young leaves were collected from the 2 month old plant for isolation of total RNA enriched with small RNA.

RNA isolation & small RNA sequencing
Extraction of total RNA including small RNA was carried out by combining the Ghawana et al. 2011 [24] protocol with miRNeasy kit (Qiagen, Germany). Fine powder of young leaf tissue (100 mg) was homogenized with solution-I and solution-II. Mixture was then treated with 200 μl chloroform to separate the organic phase. Upper aqueous phase was then separated, and 1.5 volume of 100% ethanol was added. The solution was loaded onto miRNeasy column for washing and DNase digestion. After on-column DNase digestion and washing, total RNA enriched with small RNA was eluted using nuclease free water. Further, RNA quantity and quality were examined using Nanodrop spectrophotometer and Agilent 2100 bioanalyzer system respectively. Small RNA sequencing was performed using Illumina HiSeq 2000 Platform at AgriGenome Labs Private limited, Smart City Kochi, Infopark Road, Kakanad, Kerala, India (http://www.aggenome.com/).
Raw data of small RNA sequencing was submitted to SRA (Sequence Read Archive) database of NCBI (Natioanl Centre for Biotechnology Information). Small RNA reads (50 bp length) were generated using Illumina platform. Low quality reads were removed on the basis of GC content, average base quality, and phred score. Raw reads were then processed and 5′ and 3′ adapter sequences were eliminated using Cutadapt tool (v-1.3) [25]. Remaining reads were then aligned against a number of databases like GtRNAdb, Rfam, piRNABank, siRNAdb, NCBI Genbank, deepBase to eliminate other non-coding RNAs like tRNA, rRNA, piRNA, siRNA, snRNA, and snoRNA, respectively using Bowtie2 program (version 2.1.0) [26]. The remaining reads were used to predict known and novel miRNAs.

Identification of known and novel miRNAs
The filtered small RNA reads (17-35 bp) were used for known miRNA identification. These reads were aligned with mature miRNAs and precursor sequences of viridiplantae dataset in miRBase-21 [27] using Bowtie program (version 0.12.9) allowing only 2 mismatches with mature miRNAs [28].
Unaligned reads left after known miRNAs prediction were taken further for novel miRNAs prediction using miRDeep2 (version 2.0.0.7) [29]. Reads were aligned to leaf and root transcriptome sequences of C. borivilianum. These reads were checked for hairpin loop formation using RNA folding form (version 2.3 energies) application from Mfold webserver [30]. miRNA and miRNA* sequences were found in the arm of the stem-loop. Ultimately, these reads were considered as the precursor of novel miRNAs. Sequences with the following characteristics can be considered as potential precursor of novel miRNAs: (a) sequence forming stem-loop hairpin structure (b) sequence of mature miRNA within one arm of hairpin structure (c) less than 6 mismatches are allowed between miRNA and miRNA* (d) secondary structure with minimum fold free energy lower than or equal to -18 kcal/mol (e) AU% within the range of 30-70% [31].
Validation of miRNA using stem-loop RT-PCR Randomly chosen 5 known miRNAs and all 5 novel miRNAs were selected for validation. miRNA-specific stem-loop RT-primers for reverse transcription were designed according to Kramer, 2011 [32] (Additional file 1). cDNA was prepared for each miRNA (iScript™ Select cDNA Synthesis Kit, Biorad) and PCR amplification of each miRNA was carried out using miRNA specific forward primer and universal reverse primer. PCR product was eluted, ligated into pGEM®-T easy vector and ligated product was transformed to competent E. coli DH5α cells. Recombinant plasmids were isolated, sequenced and analyzed.

Target prediction for identified miRNAs
In silico target gene prediction of known and novel miRNAs in C. borivilianum was done using target gene prediction softwares, psRNATarget [33] and miRanda [34]. Identified miRNAs were targeted against root and leaf transcriptome library of C. borivilianum, Oryza sativa (rice, transcript, TIGR (The Institute for Genomic Research) genome cDNA OSA1 Release 5 (OSA1R5), version 5) and Arabidopsis thaliana (transcripts, removed miRNA gene, TAIR (The Arabidopsis Information Resource) version 10, released on 2010_12_14) using psRNATargetV2 (2017 release) with default parameters [maximum expectation: 3; length for complementarities scoring (hspsize): 19; target accessibility-maximum energy to unpair the target site (UPE): 25; flanking length around target site for target accessibility analysis: 17 bp (upstream)/13 bp (downstream); range of central mismatch leading to translational inhibition: 10-11 nt]. It is likely that psRNATarget allocates more than one target for each miRNA sequence and multiple target sites of a particular mRNA molecule because it takes into account both complementarities and site accessibility of the targets. Once potential target mRNA sequences were obtained, target transcripts were subjected to BLASTX (Basic Local Alignment Search Tool) search against nr (Non Redundant) database from NCBI to predict function of potential targets. To further see the role of miRNAs in saponin biosynthetic pathway, mature known and novel miRNAs were targeted specifically against the transcripts coding for enzymes involved in MVA and MEP pathway using miRanda. The prediction was made using miRanda with detailed statistical study of minimum free energies ( Finally for all identified targets, Gene Ontology (GO) annotations were retrieved using Blast2GO 4.1 software at Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) levels [35]. Finally, identified targets were located in biochemical pathways by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

RT-qPCR based expression analysis of miRNAs and their predicted targets
To confirm relationship of miRNAs and their targets, leaf tissue of C. borivilianum at two stages were collected. One, at the young stage (2 month old) and another at dormancy stage (7 month old). Total 11 conserved miRNAs were selected that were computationally predicted to regulate the genes coding for enzymes involved in saponins biosynthesis. Primers for stem-loop RT-qPCR of conserved miRNAs were designed as described previously [32], and the primers for target genes were designed using Primer 3.0 input software (Additional file 1).
Total RNA along with small RNA isolation [24] and cDNA preparation was done according to protocol described earlier. Reaction cocktail for RT-qPCR of miRNAs and their targets were prepared using SsoFast™ Evagreen® Supermix and the experiment was performed using Bio-Rad CFX96™ Real-Time PCR system. β-actin was used as an internal control. 1 μl of each cDNA sample was used for analysis along with 10 μl of SYBER green, 7 μl of nuclease free water, 1 μl of each miRNA specific forward primer and stem loop complimentary universal reverse primer in cycling condition of hot start at 95°C for 2 min, 45 cycles of denaturation at 95°C for 15 s, annealing at 55°C for 30 s and extension at 72°C for 30 s. All experiments were run in triplicates. To finally calculate the fold change expression of miRNAs and their targets ΔΔC T method was used [36]. 2^(− ΔΔC T ) values were changed to (log 2 ) to generate fold change expression.

Known miRNA identification
Total 442 known miRNAs belonging to 47 miRNA families were identified from leaf's small RNA population of C. borivilianum (Additional file 2). Nomenclature of distinct mature miRNAs originated from same precursor were differentiated by adding a number suffix [37,38]. By this method, miRNA species diced from a common precursor can be differentiated. During our research, a numbers of homologs were found for each miRNA, miRNAs having more than one homolog but similar sequence were considered in same family. miR166 family with 56 member followed by miR159 (52 members) family were found to have maximum number of members in the library, although 14 miRNA families such as miR393, miR444, miR473, miR531, miR1425, miR1862, miR1873, miR3623, miR3634, miR5072, miR5077, miR7486, miR9662, and miR9674 were found to have only one member. Members of same miRNA family do not reflect the degree of nucleotide (nt) sequence similarity, but the functional equivalence i.e. to have shared functions. Number of members in each family are mentioned in Fig. 1. Total read count for each miRNA was observed, which provided an idea about apparent expression level of miRNA in C. borivilianum young leaf tissue. Out of all, on the basis of bioinformatic analysis maximum expression was observed for miR159 family i.e. 315,441 reads followed by miR166 and miR167 family with 56,445 and 25,592 reads respectively and 17 miRNA families were reported to have less than 10 reads. This wide range of expression of different miRNAs predicts that different miRNAs play distinct role in plant growth and development. Graphical representation of abundance of each miRNA family is mentioned in Fig. 2. In C. borivilianum, miRNA length was found in the range of 18-24 nt, miRNAs containing 21 nt were the most abundant (31.45%), followed by miRNAs of length 20 nt (Additional file 3).

Novel miRNA identification
Due to unavailability of whole genome sequence of C. borivilianum, root and leaf transcriptome data was used as reference for prediction of novel miRNAs. In total, 5 novel miRNAs were scrutinized (Table 2), 2 were aligned with leaf transcriptome and 3 with root transcriptome. Interestingly one miRNA i.e. cbo-miR1 was observed to have different precursor in leaf and root with same level of expression. This suggests that cbo-miR1 functions in  both root and leaf but have discrete origin. Identified novel miRNAs were ranged from 21 to 24 nt. miRNA cbo-miR3 was found to be most abundant. Variation in the abundance of plant specific miRNAs was observed in a wide range i.e. from 28 reads for cbo-miR2 to 23,252 reads for cbo-miR3. Contrary to many previous studies, sequence analysis of novel miRNAs indicated that none of them started with 5'U. Ability of miRNA precursor to form a stem-loop can confirm its biosynthesis of mature miRNA. For that precursors of predicted novel miRNAs were used for secondary structure formation using RNA folding form (version 2.3 energies) from Mfold web server. This gives dG (free energy) for each folding structure. Secondary structures of all 5 novel miRNAs are described in Fig. 3. cbo-miR3 (dG = −63.60 kcal/mol at 37°C) was found to be the most stable novel miRNA on the basis of free energy and cbo-miR5 (dG = −35.70 kcal/mol at 37°C) with minimum stability in C. borivilianum. Stem-loop RT-PCR was carried out to confirm the in silico identified known and novel miRNAs. Cloning of miRNAs was done and sequencing of recombinant plasmids verified the computationally identified miRNAs (Additional file 4). In silico analysis followed by miRNA cloning and sequencing confirmed that novel miRNAs were indeed part of miRNAs family of C. borivilanum.

Target prediction and gene ontology
To better understand the function of identified miRNAs in C. borivilianum, target prediction of known and novel miRNAs was carried out. To annotate the potential targets, BLASTX search against nr (Non Redundant) database was performed. These targets were annotated as transcription factors like GAMYB (Gibberellin and Abscisic acid-regulated MYB) transcription factor, Auxin Response Factor (ARF) family proteins, Growth regulating factors, D-tyrosyl-tRNA(Tyr)deacylase, Squamosa promoterbinding-like (SPL) protein, UDP-glucuronate decarboxylase,  Additionally, novel miRNAs were also found to regulate important biological processes. Total 7 targets for 4 novel miRNAs from leaf and root transcriptome of C. borivilianum were found. Cbo-miR3 regulates the expression of genes coding for transketolase, cytochrome P450 family proteins, and epoxide hydrolase A. And novel miRNAs cbo-miR4 and cbo-miR5 target the transcripts coding for cytochrome P450 family protein, glutathione S-transferase-1 respectively, and cbo-miR1 targets betaglucosidase 12-like isoforms. This suggests, miRNAs with their targets directly involved in the saponin biosynthesis are described in Additional file 8.
Gene functional annotation was performed by Blast2GO 4.1 software. Here we categorized the targets on the basis of their involvement in cellular process, metabolic process, and biological role. Total 135 targets predicted by psRNA-Target and miRanda were subjected to BLASTX search, mapping, and annotation to analyze and categorize using GO analysis (Fig. 4).

Relative expression profiling of miRNAs and their targets
To study the association of miRNAs and their corresponding target genes, expression pattern of 11 selected target genes of conserved miRNAs were validated by RT-qPCR. The list of genes along with their targeting miRNAs is mentioned in Table 3. Relative expression study in leaf tissue of young plant and the mature plant reaching dormancy stage was carried out. Young plant was considered as reference to calculate fold change expression. It was observed that the genes involved in saponin biosynthesis were having lower expression during dormancy than active growth period. Simultaneously, the miRNAs expression pattern was observed just opposite to that of targets. The expression of Cytochrome p450 90b1-like (CYP) was down-regulated by 7.9 log fold while its corresponding miRNA (miR396e-5p.4) was up-regulated by 1.5 log fold. On the same pattern log fold down-regulation of gene coding for a key enzyme of MVA pathway that is Squalene synthase (SQS) was by 2.8 and the comparative up-regulation of targeting miRNA (miR166i-3p) was by 1.9 log fold. This converse pattern of fold change expression of all targets and their targeting miRNAs is well explained in Fig. 7. This link of miRNAs and their targets suggest the role of miRNAs in regulation of genes controlling the saponin synthesis in plant. Relative change in expression of miRNAs and their targets was not exactly same, this means that each target can be regulated by many miRNAs [39].
Discussion miRNAs are small, non-coding, single stranded regulatory elements, subjected to deep research since last two decades. Recently, identified plant miRNAs were found to help in plant adaptation during stress conditions. These findings concluded that altered expressions of miRNAs regulate plant growth and development in several  plant species subjected to abiotic stress conditions such as drought, salinity, extreme temperatures, nutrient deprivation, and heavy metals. So, miRNAs can be used further for genetic manipulations to make crop plants more stress tolerant. In the past few years, scientific endeavors have been directed towards understanding the post-transcriptional regulation of secondary metabolites involving miRNAs. Effective role of miRNAs to directly alter the plant biochemicals are already reported, understanding of this mechanism will help to further improve secondary metabolite concentration in plants. First such report was published in 2011, which disclosed the role of miR156 in regulating the amount of anthocyanin by targeting the SPL genes [40]. After this, a number of studies reported the role of miRNAs in regulation and biosynthesis of many secondary metabolties such as flavonoids, terpenoid and alkaloids [41]. Zhang et al. 2012 have shown in their landmark study that miRNA from food products can travel safely to mammalian gut and enter the bloodstream. This suggests the plant miRNAs can be used to regulate expression of target genes even for the treatment of many human diseases. Hence, miRNAs are reported to be used as bioengineering tool to alter gene expression in plants and animals [42].
Till now, many plants have been explored e.g. in rice, miR168 was found to bound human/mouse low-density lipoprotein receptor adapter protein 1 (LDLRAP1) mRNA and reduce its protein level in liver, and consequently decrease LDL removal from mouse plasma [43]. In Curcuma longa, ath-miR167d homolog was found to hybridize with EIF2AK2 (Tyrosine-protein kinase) and ZFYVE16 (Zinc Finger FYVE-Type) genes responsible for blocking the pathway of protein processing in endoplasmic reticulum and regulation of membrane trafficking in the endosomal pathway. This can be helpful in treating diseases like arteriosclerosis and hyperglycemia [44]. In Gmelina arborea, 6 putative miRNAs were found to play a significant role in preventing diseases like cancer, blood borne disease, and other urinary infections [45]. Apart from this, 36 medicinal plants belonging to families like Fabaceae, Asteraceae, Brassicaceae, Theaceae, Caricaceae, Apocyanaceae, Rutaceae, Rubiaceae, Zingiberaceae, Scrophulariaceae, Myrtaceae, Verbenaceae, Linaceae, Euphorbiaceae, Solanaceae, Araliaceae, Coniferophyta, Salicaceae, Lamiaceae, Malvaceae, Ericaceae, Vitaceae, and Gramineae have been studied to explore their miRNA pool [46]. This is the first attempt to study miRNA of a medicinal monocot herb from Liliaceae family. Extensive therapeutic and medicinal  properties of C. borivilianum made it a choice for miRNA study. C. borivilianum has been focused for the study of secondary metabolites biosynthetic pathways since a decade. Till now, miRNAs and their function in C. borivilianum are not known. Knowledge of miRNA can help in understanding the loop hole in other molecular and metabolic studies in this plant. In this study, we applied deep sequencing combined with bioinformatic analysis to identify and characterize miRNAs and their targets in C. borivilianum.
There is no information present worldwide about sRNA of any liliaceae family plant. So this study may serve as a foundation for further exploration and application of complex metabolic mechanism of bioactive substances found in C. borivilianum. Larger amount of data produced by deep sequencing helped to identify even less abundant and small sized miRNAs, because analyzing larger amount of reads increases the odds of recovering rare transcripts [47]. A total of 442 known miRNAs belonging to 47 families and 5 novel miRNAs were identified from C. borivilianum leaf sRNA library. It signifies that, diverse and highly complex small RNA population exists in C. borivilianum. An almost similar trend has been observed in other monocot plants e.g. in Aegilops tauschii, Brachypodium distachyon, O. sativa, Sorghum bicolor, Triticum aestivum, and Zea mays, number of mature miRNAs reported are 173, 525, 713, 241, 119, and 321, respectively [27]. Although the miRNA gene sequences may vary, the seed regions of the miRNAs belonging to a same family can be identical [48]. Another difficulty regarding miRNA annotation is the repetitive presence of paralogous MIRNA loci in genome producing identical or nearly identical mature miRNAs. Based on this rationale, these miRNAs can be categorized into same families [49]. Highly conserved miRNAs like miR156, miR157, miR159, miR160, miR162, miR164, miR165, miR166, miR167, miR168, miR169, miR171, miR172, miR319, miR390, miR393, miR394, miR395, miR396, miR397, miR398, miR399, miR408, and miR444 in other species [50] were also observed in C. borivilianum in high abundance. Novel miRNAs or species specific miRNAs considered to play specific functions in plants whereas conserved miRNAs are thought to be involved in generalized functions like signal transduction, development of leaf and flower etc. [51]. Conserved miRNAs have unique property of being represented by multiple loci in sequenced genomes; most of them were generated through genome duplication events, giving some indication of their antiquity [52]. Conservation pattern of miRNAs vary from species to species. Most often conservation of miRNAs is associated with its sequence features such as base content and cleavage sites. Plant miRNAs show a negative/positive correlation between conservation and AU/GC content. At the 5′ end, conserved miRNAs usually starts with base U, while less-conserved miRNAs have a non-U base at start position in mammals. But this is not true in case of insects and plants [48]. On the basis of data analysis, it can be predicted that miR159 and miR166 have maximum expression in leaf during the period of its active growth. Highest abundance of same miRNAs is also reported in other plants like Panax ginseng, Stevia rebaudiana, and A. thaliana [53][54][55]. miR159 represents one of the most ancient miRNAs in the plant kingdom [52]. miR159 family members were reported to regulate ABA stress response and seed germination in plants by regulating the level of MYB transcription factor in A. thaliana [56,57]. In tomato, miR159 regulates leaf and flower development by targeting the SGN-U567133 [58]. As the sample was collected during the active plant growth period, high abundance of miR159 suggests its active regulation of leaf and root development.
Monocot specific miRNAs such as miR437, miR444, miR396 were reported in monocot plant species like rice, maize, sorghum, and sugarcane [50]. In our study, we found the presence of only miR396 and miR444 but not miR437.
Complete reference genome should be the primary limiting factor for sRNA-seq. Unfortunately, the genome of C. borivilianum has not been published, but transcriptomic study of root and leaf tissue of C. borivilianum has already been conducted. Transcriptomic data from previous studies was used as reference for novel miRNA prediction. But many more miRNAs from C. borivilianum can be annotated in future when complete genome information will be available. Further, our study revealed 5 putative novel miRNAs in C. borivilianum using the available sequence data of this plant as reference. MFE (Minimum Fold Energy), sequence length and base composition are important features in predicting plant miRNAs. So, predicted novel miRNAs were confirmed by forming secondary structure of their precursor. All novel miRNA precursors have stem loop hairpin structure, and this foldback hairpin structure has a low free energy as predicted by Mfold software. Novel miRNAs in C. borivilianum found to originate from precursor with AU percentage ranged from 40 to 53%. Seed region of miRNA can bind to the 3'-UTR, 5'-UTR, and ORF region of target mRNA but out of all these 3'-UTR targeting is much more frequent than other two [59].
In the present study, along with transcriptomic data of leaf and root of C. borivilianum, O. sativa and A. thaliana transcriptome have also been used to find maximum miRNAs targets. Target inhibition by cleavage was found to be more frequent. This result supports the fact that cleavage by Argonaute2 mediation maybe the main mode of gene suppression for many known plant miRNAs [60]. Moreover, we found that in A. thaliana ARF6 and ARF8 were targeted by miR167 and ARF10, ARF16, and ARF17 by miR160 and in O. sativa ARF16 by miR160; ARF6, ARF8 were targeted by miR167. In C. borivilianum ARF18 and ARF17 were found to be targeted by miR160; and ARF12 by miR167. Down-regulation of ARF6 and ARF8 by miR167 in tomato, leads to affects the floral development and female sterility [61]. This suggests the functional similarity of miRNAs in C. borivilianum, O. sativa, and A. thaliana.
Target prediction results were further analyzed using Blast2GO. Biological processes like signaling, response to stimulus, developmental process, reproductive process, were found to be regulated by miRNAs. This suggested the wide regulation by miRNA in C. borivilianum. The two most important pathways identified were terpenoid backbone biosynthesis and sesquiterpenoid and triterpenoid biosynthesis. This suggested the specific role of miRNAs in secondary metabolism. Predicted target genes and the functions of their products may provide valuable clues for research into essential biological processes and metabolism in C. borivilianum.
RT-qPCR was carried out to calculate the expression levels of miRNAs and their target mRNAs. The expression was measured in leaves at two growth stages of plant. All 11 miRNAs were found with increased expression during dormancy. Maximum variation was observed in the expression of miR172c-5p and miR398a-3p.5 which are targeting 4-Hydroxy-3-methylbut-2-enyl diphosphate reductase. This concluded that multiple miRNAs can regulate same gene collectively. Respective up and down regulation of miRNAs and their targets during dormancy indicate vital role played by miRNAs in regulation of secondary metabolite accumulation. It ensures that miRNAs are directly involved in negative regulation. This present study follows the same scenario as discussed previously by some other investigators [41].

Conclusion
Total 442 known and 5 novel miRNAs were identified with the help of bioinformatic tools. The plant specific miRNA and known miRNAs were found to be associated with MVA/MEP pathway. This helped to establish a correlation between the miRNA and secondary metabolism. miRNAs targeting genes involved in secondary metabolite biosynthetic pathway can be used to further elaborate the regulation of genes by miRNAs to moderate their expression. Our study also found some mRNA coding transcription factors as target of miRNAs which suggested the role of miRNAs in plant growth and development. GO and KEGG analysis provide foundation for further research. RT-qPCR results imply that during dormancy, synthesis of secondary metabolites slows down and helped to understand the role of miRNAs in secondary metabolism. This will facilitate to develop database of C. borivilianum useful for designing RNAi experiments to regulate secondary metabolite content.