Modular Design of Picroside-II Biosynthesis Deciphered through NGS Transcriptomes and Metabolic Intermediates Analysis in Naturally Variant Chemotypes of a Medicinal Herb, Picrorhiza kurroa

Picroside-II (P-II), an iridoid glycoside, is used as an active ingredient of various commercial herbal formulations available for the treatment of liver ailments. Despite this, the knowledge of P-II biosynthesis remains scarce owing to its negligence in Picrorhiza kurroa shoots which sets constant barrier for function validation experiments. In this study, we utilized natural variation for P-II content in stolon tissues of different P. kurroa accessions and deciphered its metabolic route by integrating metabolomics of intermediates with differential NGS transcriptomes. Upon navigating through high vs. low P-II content accessions (1.3–2.6%), we have established that P-II is biosynthesized via degradation of ferulic acid (FA) to produce vanillic acid (VA) which acts as its immediate biosynthetic precursor. Moreover, the FA treatment in vitro at 150 μM concentration provided further confirmation with 2-fold rise in VA content. Interestingly, the cross-talk between different compartments of P. kurroa, i.e., shoots and stolons, resolved spatial complexity of P-II biosynthesis and consequently speculated the burgeoning necessity to bridge gap between VA and P-II production in P. kurroa shoots. This work thus, offers a forward looking strategy to produce both P-I and P-II in shoot cultures, a step toward providing a sustainable production platform for these medicinal compounds via-à-vis relieving pressure from natural habitat of P. kurroa.


INTRODUCTION
Picrorhiza kurroa, locally named as Kutki, is a high altitude (3000-5000 m) medicinal herb distributed in the North-Western Himalayas of India. It is used for the treatment of various ailments due to pharmacological properties like hepatoprotective, antiallergic, antiasthmatic, antioxidant, anticancerous, and immunomodulatory and therefore, consequently provide livelihood and health security to a large segment of high altitude populations (Kumar et al., 2015a(Kumar et al., , 2017. It is listed among top 15 most traded medicinal herbs in India with respect to revenues generated by traded material (Uniyal et al., 2011). P. kurroa is a self-propagating medicinal herb. It originates as a young bud on stolons which develops to produce mature stolon and subsequently rhizomes with self-governing roots and shoots compartments (Pandit et al., 2013a). Thus, during the collection of mature rhizomes of P. kurroa, all the young stolons and buds are also get uprooted which subsequently disrupt it's proliferation in natural habitat. The extensive harvesting has thus put this herb in the category of critically endangered medicinal herbs and as a result legal restrictions have been levied on its collection from the wild (Kumar et al., 2017). The scarcity of herbal raw material not only caused economic constraints on local communities but also resulted in adulteration of P. kurroa, thereby affecting quality and efficacy of drugs due to lower levels of picroside-I (P-I) and picroside-II (P-II), the signature bioactive compounds of P. kurroa (Shitiz et al., 2013).
The picrosides, P-I and P-II, possess vast therapeutic potential and their bio-efficacy has also been tested via both in vitro as well as in vivo studies (Kumar et al., 2017, Table 1). Moreover, the herbal drug formulations required both P-I and P-II in a ratio of 1:1.5 (Kumar et al., 2017). Therefore, to provide a constant supply of raw material without sacrificing the levels of major chemical constituents, mass propagation of quality plant material is a pressing need. Mass propagation of P. kurroa shoots through in vitro culture has been optimized but with lower yields of P-I (Sood and Chauhan, 2010). Moreover, the content of P-I is also ∼5-fold lesser in P. kurroa shoots grown in vitro compared to those available in natural habitats (Sood and Chauhan, 2010;Shitiz et al., 2015;Sharma et al., 2016).
In P. kurroa, it has been hypothesized that the P-I and P-II are biosynthesized in a tissue specific manner, i.e., shoots synthesize  Comu et al., 2016 Anti-allergic asthma Choi et al., 2016 Renal ischemia/reperfusion injury  Neuroprotective Guo et al., 2010;Li et al., 2010;Pei et al., 2012;Liu et al., 2014;Zhao et al., 2014; Hepatoprotective Gao and Zhou, 2005a,b Hypoxia/reoxygenation induced cardiomyocyte injury Meng et al., 2012a,b Myocardial ischemia reperfusion injury Wu et al., 2014 P-I Hepatoprotective Dwivedi et al., 1992 Anti-cancer Rathee et al., 2013 P-I and P-II Neuritogenesis Li et al., 2000Li et al., , 2002 P-I, roots contained only P-II and stolons/rhizomes produce both P-I and P-II (Pandit et al., 2013a). The shoots of P. kurroa have been established as biosynthetic tissues of only P-I as P-II was not detected in P. kurroa shoots and the P-I content also showed progressive increase in different growth and developmental stages of P. kurroa shoots (Sood and Chauhan, 2010;Pandit et al., 2013a;Kumar et al., 2015b). In contrast, both P-I and P-II contents have been detected and showed progressive increase in different growth and developmental stages of P. kurroa stolons/rhizomes. In previous reports, we had suggested that rhizomes are a storage tissue for P-I and P-II not the biosynthetic tissue, whereas P-I was postulated as biosynthesized in shoots and possibly transported to rhizomes (Sood and Chauhan, 2010;Pandit et al., 2013a). The root tissues of P. kurroa contained only P-II which was detected at the later stage compared to rhizome and thus, it was hypothesized that P-II might be transported from rhizomes to roots (Pandit et al., 2013a). Therefore, it would be of paramount importance if P-II can be produced in shoots, which will not only relieve pressure from natural habitat of P. kurroa for uprooting rhizomes/roots but will also provide a continuous supply of herbal raw material through shoot cultures with desired amounts of P-I and P-II. However, lack of knowledge on P-II biosynthesis has precluded implementing any genetic intervention strategy toward enriching P-II in P. kurroa shoot cultures.
In previous studies, maximum efforts have been focused on P-I production in P. kurroa as P-II was not produced under in vitro shoot cultures. The biosynthetic pathways for P-I and P-II were, therefore, proposed including; mevalonate (MVA), non-mevalonate (MEP), shikimate/phenylpropanoid, and iridoid pathways so as to pinpoint key genetic intervention points to maximize picrosides content using metabolic engineering approaches (Kumar et al., 2013;Shitiz et al., 2015). It has been hypothesized that the catalpol (CAT) backbone derived from iridoid pathway undergoes esterification reactions with transcinnamic acid (CA) and vanillic acid (VA) to produce P-I and P-II, respectively (Kumar et al., 2013). The metabolic route of P-I biosynthesis has been ascertained through precursor feeding and enzyme inhibitor studies along with the assortment of key candidate genes through expression analysis using qRT-PCR visà-vis P-I content (Pandit et al., 2013b;Kumar et al., 2015bKumar et al., , 2016bShitiz et al., 2015). However, the complexity of P-II biosynthesis is still unresolved as both protocatechuate and ferulic acid (FA); the metabolites of shikimate/phenylpropanoid pathway, are plausible to produce VA, the immediate precursor of P-II (Figure 1) (Kawoosa et al., 2010;Gahlan et al., 2012;Kumar et al., 2013). In this study, we have selected mature stolons to elucidate the mechanism controlling P-II biosynthesis in P. kurroa.
To elucidate the metabolic fluxes, radioactive/stable isotopes labeled flux analysis has been used in various plant species (Ishihara et al., 2015;Xiong et al., 2015) but concurrently accompanied with cost intensive and time consuming experiments (Stephanopoulos, 2012). Therefore, in this study, we have utilized natural variation in P-II content among different accessions (genotypes) of P. kurroa so as to guide us in discerning the route and flux of intermediates leading to the biosynthesis of P-II. Moreover, the integration of metabolites contents data with FIGURE 1 | Modular design of P-I and P-II biosynthetic pathways. The structures of P-I and P-II linked to cinnamic acid/vanillic acid (purple/blue color) moieties and catalpol (orange color). Question marks indicate plausible routes for vanillic acid production. the corresponding transcriptomes would enable us to understand the flux of intermediates vis-à-vis gene transcripts controlling P-II biosynthesis in P. kurroa.

Collection of Plant Material
A total of 10 accessions of 1 year old P. kurroa were collected from different locations of North-Western Himalaya, India. The details of accessions are provided in Supplementary Table 1. All the accessions were kept in the greenhouse of Jaypee University of Information Technology, Waknaghat, Himachal Pradesh, India (31 • 0 58.55 N and 77 • 4 12.63 E; 1700 m, altitude). The plants were washed with tap water and segregated into shoots, mature stolons, and roots. The shoots were coded with PKS, while mature stolons were labeled as PKST. In contrast, roots were labeled as PKR. We have used PKS and PKST tissues in this study. All tissue samples were kept at −80 • C for further analyses.

Extraction and Quantification of P-I and P-II Contents
The P-I and P-II were extracted in triplicates from shoots and stolons tissues of all 10 P. kurroa accessions as per the protocol described in Kumar et al. (2016c). The quantification of P-I and P-II in all the samples was performed in triplicates by using HPLC method as described in Pandit et al. (2013a).

Extraction and Quantification of Intermediate Metabolites
Four accessions with variation for P-I content in shoots (PKS-1, PKS-4, PKS-5, and PKS-21) and with variation for P-I and P-II contents in stolons (PKST-3, PKST-5,  were subjected to the extraction of intermediate metabolites viz. CA, p-coumaric acid (p-CA), FA, protocatechuic acid (PA), and VA. The extraction was performed as per the method described in Mattila and Kumpulainen (2002) with some modifications. The samples (200 mg) were weighed and homogenized with liquid N 2 in prechilled pestle and mortar. To each homogenized sample, 7 mL of combination of methanol [containing 2 g/L butylated hydroxytoluene (BHT)] and 10% acetic acid in a ratio of 85:15 was added and vortexed for 10 min at room temperature (25 • C). The volume of extracts was made up to 10 mL with distilled water and mixed properly. Each extract was then transferred to 100 mL reagent bottle and suspended with 12 mL of distilled water and 5 mL 10 M sodium hydroxide. The contents of the extracts were then bubbled with N 2 gas, sealed and kept on a magnetic stirrer for overnight at room temperature (25 • C).
On following day, the pH of extracts was adjusted to 2 with concentrated hydrochloric acid. The metabolites released in the solution were then extracted three times with 15 mL of a mixture of cold diethyl ether and ethyl acetate in a ratio of 1:1 by shaking and centrifuging for 10 min at 7000 rpm and 4 • C temperature condition. The diethyl ether/ethyl acetate layers thus obtained, were evaporated to dryness and reconstituted in 1.5 mL of absolute methanol. The obtained samples were filtered through 0.22 µm filter (Millipore) and subjected to quantification of selected metabolites by using HPLC as per the method described in Nour et al. (2013). The selected compounds viz. CA, p-CA, and FA were detected at a wavelength of 290 ± 4 nm, 310 ± 4 nm, and 330 ± 4 nm, respectively while PA and VA were detected at an absorbance of 260 ± 4 nm. All the selected compounds were identified upon comparison of retention time and UV spectra with standards (Sigma-Aldrich, USA). The extraction and quantification experiments were performed in triplicates.

Feeding of In vitro Grown P. kurroa Shoots with Different Precursors
The solutions of different precursors including, p-CA, PA, and FA (Sigma-Aldrich, USA) were prepared as neutral aqueous stocks at 150 µM concentrations. Conversely, 150 µM aqueous solution of VA (Sigma-Aldrich, USA) was mixed with 70 µM CAT (Sigma-Aldrich, USA) and used for the feeding of in vitro grown P. kurroa shoots. These concentrations were selected based on their effects on P-I production in P. kurroa when applied alone and in combination under tissue culture conditions (Kumar et al., 2016b). Following this, the solutions of different precursors were added into the culture tubes containing optimized MS medium [containing 3 mg/L indole-3-butyric acid (IBA) and 1 mg/L kinetin (KN)] after filter sterilization through 0.22 µm syringe filter (Millipore) . In each culture flask, the shoot explant of P. kurroa grown at 25 ± 2 • C (having negligible P-I content) was aseptically transferred and all the cultures were kept in a plant tissue culture chamber maintained at 15 ± 2 • C (temperature favors P-I production) (Sood and Chauhan, 2010). The P. kurroa shoots without precursor's treatment were referred to as controls. The experiment was performed in triplicates. The shoot samples were collected after 30 days and immediately stored at −80 • C for further analyses. All the samples were subjected to the extraction and quantification of CA, p-CA, FA, VA, P-I, and P-II metabolites in triplicates as per the protocol mentioned above.

Culturing of P. kurroa Shoots Grown
In vitro with Different Concentrations of VA+CAT Mixture and P-II The different concentrations of VA+CAT mixture viz. 25 µM each VA+CAT, 70 µM each VA+CAT, 70 µM VA with 150 µM CAT, 150 µM each VA+CAT and 230 µM each VA+CAT, were also tested after determination of P-II content in shoot cultures fed with different precursors. Moreover, P-II at a concentration of 25 µM was also introduced separately in the MS medium to check the accumulation of P-II in P. kurroa shoots. The culturing of P. kurroa shoots was performed according to the same protocol as mentioned above. The untreated shoots were labeled as controls. The experiment was performed in triplicates. The shoot samples were harvested after 30 days and immediately stored at −80 • C for the analysis of P-I and P-II contents. The analysis of picrosides was carried out in triplicates.

Treatment of In vitro Grown P. kurroa Shoots with VA+CAT Mixture and P-II under Liquid Culture Conditions
The solutions of VA+CAT mixture and P-II were prepared at concentrations of 70 and 25 µM, respectively to check the intake of precursors by the plant. These solutions were then incorporated into the culture tubes having optimized liquid MS medium (containing 3 mg/L IBA and 1 mg/L KN) after filter sterilization through 0.22 µm syringe filter (Millipore). The culturing of P. kurroa shoots was performed under the same conditions as mentioned above. The untreated shoots were labeled as controls. The experiment was performed in triplicates. The shoot samples were harvested after 30 days and immediately stored at −80 • C for the further analyses. The analysis of VA and P-II contents was performed in both shoots and liquid media left after the harvesting and carried out in triplicates.

Isolation of Total RNA, Preparation of Libraries, and Sequencing for RNA-seq Analysis
Total RNA was isolated from four accessions of shoots (PKS-1, PKS-4, PKS-5, and PKS-21) and stolons (PKST-3, PKST-5, PKST-16, and PKST-18) by using total RNA isolation kit (RaFlex TM ) according to the instructions mentioned in the user manual. The construction of libraries and their sequencing for transcriptome analysis was performed as per the protocol described in Pal et al. (2015).

Transcriptome De novo Assembly and Validation
For de novo assembly, we chose go with De brijn graph based Trinity Assembler based on the criteria of; (a) default K-mer, (b) less memory foot print, (c) optimized for Illumina paired end data, (d) reproducibility, and (e) configurable for all computing capacities (Henschel et al., 2012). To pool assembly, 50 core threads of processing with 2.4 GHz speed and a maximum Heap Space of 50 GB was allotted. Nevertheless, de novo transcriptome assemblers are capable of producing in fragmented/missassembly, the validation of assembled transcriptome was carried out by mapping back the high quality (HQ) filtered reads to the ESTs. The HQ filtered reads from each library thus obtained were mapped to the assembled transcriptome by using Bowtie software. This analysis generated.bam (Binary Sequence Alignment/Map) files which were processed through bedtools and samtools for quantitation (read count estimation) of each transcript in a library and also to calculate the total coverage and average depth of the transcriptome in each library (Li et al., 2009;Quinlan and Hall, 2010).

Transcripts Annotation
The annotation and statistical analysis of each transcript was performed by homology based method against NCBI nrdb protein database using Blast2GO software through Blastx module (Conesa et al., 2005). The parameters used for this analysis were; (a) e-Value b = 10-e4, (b) Similarity ≥35%, (c) Annotation cutoff ≥55, (d) GO weight cutoff ≥5. Moreover, domain level annotation was also performed by using the Online InterProScan tool RunIprScan-1.1.0 1 (Henschel et al., 2012).

Normalization of Transcripts and their Expression Profiling
The benefit of using NGS based transcriptome profiling is to recognize sample/condition specific expressed transcripts which was not easy with earlier hybridization methods. Transcripts with a read count of ≥10 in any one of the libraries were considered to be expressed. To perform normalization and expression profiling, a sub bam file was created from the master bam file using RSEM software (Li and Dewey, 2011). In RSEM, the default parameters/commands were used to normalize each library which provides an output with expected normalized read count, TPM (tags per million) and FPKM (fragments per kilobase per million). Log to the base 2 of FPKM was considered as absolute expression or Delta CT equivalent value.

Transcripts Mapping, Differentially Expressed Gene (DEG), and Correlation Analysis
The annotated transcripts obtained from BLAST2GO were mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016). Transcripts which were involved in secondary metabolic pathways were shortlisted for DEG analysis using sequence similarity threshold of >98%. Only those transcripts were selected and used for the analyses which were present in all the eight different samples, i.e., four accessions of each shoots and stolons. Hierarchical clustering of selected DEGs were performed using heatmaply, gplot2, plotly packages in R. Further, the correlation analysis was carried out by using R-package through "Corrgram" function (Friendly, 2002).
genes were selected and subjected to qRT-PCR analysis. The primers of the selected genes were designed from transcriptomic sequences of P. kurroa through Primer3 software and the details with annealing temperature are provided in Supplementary  Table 2. The gene expression analysis was performed in triplicates according to the method described in Kumar et al. (2015b). The housekeeping gene, 26S was used as reference in this study for normalization of gene expression.

Statistical Analyses
The one and two-way analysis of variance (ANNOVA) followed by a Bonferroni test were performed from data in triplicates (mean ± SD) using GraphPad prism software version 6.0. To examine the correlations between gene expression profiles obtained from RNA-seq and qRT-pCR, log2 fold change values were estimated among different shoots and stolons accessions. The scatterplots were then created by comparing log2 fold change values determined by RNA-seq and qRT-PCR using GraphPad prism software version 6.0. The correlations were represented in terms of coefficient of determination, i.e., R 2 .

Alterations in the Levels of Intermediate
Metabolites of Shikimate/Phenylpropanoid Pathway among Selected P. kurroa Accessions The shoots and stolons of P. kurroa accessions selected for high and low P-I/P-II contents were examined for the levels of intermediate metabolites in shikimate/phenylpropanoid pathway viz. CA, p-CA, FA, PA, and VA (Figures 3A-F). The investigation of shoot tissues revealed that PKS-1 contained maximum CA content of 0.5% which showed significant progressive decrease in PKS-5, PKS-4, and PKS-21 with 1.5-(p < 0.0001), 2.5-(p < 0.0001), and 5.3-fold (p < 0.0001), respectively ( Figure 3E). These results are akin to P-I accumulation patterns in selected shoot tissues of P. kurroa. In contrast, the lowest p-CA content (0.04%) was observed in PKS-1 and showed significant increase in PKS-5 and PKS-21 with 2.1-(p < 0.05) and 5.9-fold (p < 0.0001), respectively. The p-CA content of PKS-4 showed non-significant modulation compared to PKS-1 shoots ( Figure 3E). The examination of FA content in shoot tissues revealed 2-fold (p < 0.05) significant elevation in PKS-21 compared to PKS-1, PKS-5, and PKS-4 accessions ( Figure 3E). Moreover, the VA content showed non-significant modulation among all the four selected shoots of P. kurroa. The PA content was not detected in the selected shoots of P. kurroa. It was evident from these results that CA is positively correlated with P-I content in shoots with Pearson correlation coefficient (PCC) of 0.99 ( Figure 3F).
The investigation of four different stolons revealed that PKST-18 contained maximum CA content of 0.2% which showed significant increase with 4.7-fold (p < 0.0001) compared to PKST-3 ( Figure 3B). The CA content of PKST-3 showed non-significant modulation with PKST-5 and PKST-16 stolons. The analysis of p-CA content showed non-significant alterations in all the four selected stolons of P. kurroa. Upon examination of FA content, we observed highest level in PKST-5 (0.2%) which exhibited 1.7-fold (p < 0.001) significant increase compared to PKST-3, PKST-16, and PKST-18 ( Figure 3B). In contrast, highest VA content was observed in PKST-3 and PKST-5 followed by significant decrease in PKST-16 and PKST-18 with 2.2-(p < 0.0001) and 1.3-fold (p < 0.0001), respectively ( Figure 3B). Moreover, PA content was not observed in the stolons of P. kurroa. These observations infer positive correlations of p-CA and VA with P-II through PCC of 0.45 and 0.79, respectively ( Figure 3C).

Effect of Treatment with Different Precursors on Picrosides and Intermediate Metabolites Contents
The in vitro grown shoots of P. kurroa were treated with 150 µM concentration of PA, FA, and p-CA along with combination of 150 µM VA and 70 µM CAT to observe their effects on the CA, p-CA, FA, VA, P-I, and P-II contents (Figures 4A-C). The analysis of CA content revealed 2.2-(p < 0.0001) and 1.7-fold (p < 0.0001) significant increase in shoots treated with FA and p-CA, respectively, whereas 1.3-fold (p < 0.001) significant decrease was observed in VA + CAT fed shoots compared to P. kurroa shoots without treatment, i.e., control shoots ( Figure 4A). The p-CA content, on the other hand showed non-significant modulations in all the fed shoots compared to control. The investigation of FA content revealed significant increase with 2.9-fold (p < 0.0001) in shoots treated with FA compared to untreated shoots. In contrast, 1.9-(p < 0.05) and 2-fold (p < 0.01) significant elevation in VA content was observed in shoots fed with FA and VA + CAT, respectively compared to control ( Figure 4A). Moreover, the examination of P-I content revealed 1.8-fold (p < 0.0001) significant enhancement in shoots fed with p-CA compared to control. Conversely, small level of P-II, i.e., 0.04% was only detected in P. kurroa shoots treated with VA + CAT ( Figure 4A). Thus, it was evident from these results that VA and P-II shared high positive correlation with PCC of 0.64 and the former was also positively correlated with CA and FA with PCC of 0.19 and 0.22, respectively ( Figure 4C). showing correlations between tested metabolites and picrosides contents among stolons and shoots, respectively. Correlations are presented in the form of graphs filled in proportion to the Pearson's correlation coefficient values. Clock-wise occupied with blue color depict positive correlations while anti-clockwise graphs filled with red color indicate negative correlations. The data presented as means ± SD (n = 3). Significance was evaluated within metabolites contents between different accessions ( * p < 0.05, * * * p < 0.001, * * * * p < 0.0001).
Frontiers in Plant Science | www.frontiersin.org FIGURE 5 | Differential expression analysis of genes involved in secondary metabolism among different tissues of P. kurroa accessions; (A,B) cluster analysis of differentially expressed genes patterns between different P. kurroa accessions selected for shoots and stolons, respectively; (C,D) correlations in terms of Pearson's correlation coefficients depicted using correlogram between selected genes and metabolites contents among different stolons and shoots, respectively. The data is represented in the form of graphs filled in proportion to the Pearson's correlation coefficient values. Clock-wise occupied with blue color depict positive correlations while anti-clockwise graphs filled with red color indicate negative correlations.
Since, P-II was only detected in VA + CAT fed shoots therefore, different concentrations of VA + CAT were further tested to observe their influence on the P-II content. The data revealed non-significant alterations in P-II content among P. kurroa shoots treated with different concentrations of VA + CAT ( Figure 4B). Moreover, the P-II content also exhibited non-significant alteration in shoots fed with P-II while 0.2% P-II was observed in the media remained after the sampling of shoots treated with P-II ( Figure 4B).
The feeding of different precursors in vitro also increased shoot biomass compared to untreated control which might be due to the enhanced carbon source utilization.

Correlation Analysis between Picrosides Content and Transcripts Abundance Values (FPKM)
To get insight into flux dynamics of different metabolic modules connected to picrosides biosynthesis, correlation maps were constructed between the FPKM values of selected transcripts and picrosides content among different shoots and stolons accessions of P. kurroa 2 . The transcripts encoding enzymes catalyzing regulatory reactions in different metabolic pathways were selected based on their connection with picrosides biosynthesis. The details of selected transcripts with their role in picrosides production are provided in Table 2.
The   (Figures 5B,C). These observations implied that DAHPS, PAL, PEPCK, and 4CL which showed high positive correlations (>0.60) with P-I content might be considered as the key candidate genes for P-I biosynthesis. Conversely, HMGR, C4H, DAHPS, and GS might play an important role in the biosynthesis of P-II in P. kurroa.
Finally, the correlations between differential gene expression patterns of six selected genes obtained from transcriptomes data and qRT-PCR were determined and we observed coefficient of determination (R 2 ) ≥ 0.75 which indicates a good positive correlation (Figure 7). It was thus evident from the results FIGURE 6 | Expression profiling of selected genes in different tissues of P. kurroa accessions; (A) shoots and (B) stolons. Expression values were normalized with levels of 26S reference gene. The data presented as means ± SD (n = 3). Significance was evaluated for each gene between different accessions ( * p < 0.05, * * p < 0.01, * * * * p < 0.0001).
that gene expression patterns investigated by qRT-PCR were corroborated with those obtained by RNA-seq analysis, thereby supporting the reliability of our transcriptomes data.

DISCUSSION
The metabolic network of picrosides production is overwhelmingly complex as multiple pathways contribute to their biosynthesis in P. kurroa. In recent years, the major attention was paid to elucidate the metabolic flux dynamics in P-I production through the use of various inhibitors and precursors specific to the enzymes catalyzing critical steps in its biosynthesis (Shitiz et al., 2015;Kumar et al., 2016b). However, the biosynthesis of P-II remains untapped since it is not produced in P. kurroa plants grown in vitro (Sood and Chauhan, 2010). Moreover, the shoots of P. kurroa plants grown in natural habitats exhibited only P-I while stolons contained both P-I and P-II (Pandit et al., 2013a). Thus, to investigate the metabolic basis of P-II biosynthesis, we have employed a strategy utilizing natural genetic diversity existing for the production of P-II among P. kurroa accessions collected from different geographic locations of North-Western Himalayas, India. Moreover, the genes/enzymes information for most of the enzymatic steps has been obtained through generating and mining NGS transcriptomes from differential content phenotypes of P. kurroa. The differential transcriptomes provided greater insights not only into the genes involved in biosynthetic pathways but also molecular mechanisms that regulate the contents of picrosides in P. kurroa. In the present study, the analysis of stolon tissues from different P. kurroa accessions revealed prominent increase FIGURE 8 | Complete framework of metabolic network depicted P-I and P-II biosynthesis in different compartments of P. kurroa. The established routes of P-I and P-II production were presented with different colors to highlight the modulations in shoots and stolons tissues of P. kurroa. Up-regulated gene expressions were highlighted with bold and green color. The abbreviations are elaborated in Supplementary Table 3 in P-II content upon reduction in P-I level indicating that both P-I and P-II biosynthesis skewed from a common metabolic node. It is likely since VA, the immediate precursor of P-II; is supposed to be formed from either PA or FA, both are derived from shikimate/phenylpropanoid pathway which also produces CA, the immediate biosynthetic precursor of P-I (Funk and Brodelius, 1992;Kumar et al., 2013Kumar et al., , 2016b. The shikimate/phenylpropanoid pathway produces numerous branchpoints such as chorismate, arogenate, CA, p-CA, caffeic acid and FA, which act as precursors for the formation of phenylpropanoids, alkaloids, flavonoids, lignins, and picrosides, i.e., P-I and P-II (Tohge et al., 2013;Alvarez et al., 2016;Kumar et al., 2016b). The chorismate is involved in the formation of both tryptophan and arogenate metabolites (Kumar et al., 2016b), the later one on the other hand, can be converted into tyrosine and phenylalanine, the branchpoints for the formation of alkaloids and phenylpropanoids, respectively (Pascual et al., 2016;Schlager and Drager, 2016). The CA formed from phenylalanine acts as immediate precursor to P-I and also produces various phenylpropanoids via the formation of p-CA (Kumar et al., 2016b;Rigano et al., 2016). Moreover, the CoA pools generated from CA, p-CA, caffeic acid, and FA can deviate P-II biosynthesis by directing the flux of shikimate/phenylpropanoid pathway toward the formation of flavonoids and lignins (Alvarez et al., 2016; Figure 8).
Therefore, it is crucial to monitor dynamic profiles of intermediate metabolites in shikimate/phenylpropanoid pathway to get deep insight for metabolic switches of P-II biosynthesis. To address this, we have analyzed stolon tissues of four P. kurroa accessions based on their varying P-II content ranged from minimum to maximum (1.36-2.65%). The results observed in this study revealed statistically significant increase in VA content among high P-II content stolon tissues (PKST-3 and PKST-5) compared to stolons possessing low P-II content  which established that VA is associated with P-II biosynthesis. This result is consistent with the proposed hypothesis that VA is the immediate precursor of P-II biosynthesis (Kumar et al., 2013). It is noteworthy that analysis of CA and p-CA metabolites content in selected stolon tissues revealed the flux direction of shikimate/phenylpropanoid pathway through to the CA step leading to the biosynthesis of p-CA and P-I as CA is located at the branching point between P-I and p-CA production (Kumar et al., 2016b). Interestingly, 87% CA was converted to p-CA in PKST-5 (P-1, 0.08%; P-II, 2.19%), 44.4% in PKST-3 (P-1, 0.17%; P-II, 2.65%), 21.1% in PKST-16 (P-1, 1.05%; P-II, 1.52%) and 8.2% in PKST-18 (P-1, 1.32%; P-II, 1.36%). This implies that with increase in the P-I content, the flux of shikimate/phenylpropanoid pathway limits through to the CA step to p-CA while increase in the P-II content directs the flux of shikimate/phenylpropanoid pathway through to the CA to p-CA leading to the enhanced P-II production in P. kurroa (Figure 8). This statement is also supported by comparative transcriptome analysis of selected stolon tissues which showed high positive correlation of gene encoding cinnamic acid-4-hydroxylase (0.60) with P-II content compared to gene encoding phenylalanine ammonia lyase enzyme which showed PCC of 0.28 with P-II content ( Figure 5C). The cinnamic acid-4-hydroxylase is an enzyme that catalyzes the conversion of CA to p-CA while phenylalanine ammonia lyase catalyzes the conversion of phenylalanine to CA, the immediate precursor of the P-I biosynthesis (de Jong et al., 2015;Kumar et al., 2016b). Therefore, based on this analysis, we have hypothesized that CA is the common metabolic node between P-I and P-II biosynthesis in P. kurroa rather than 3-dehydroshikimate which serves as the branch point between PA and shikimate production; the metabolites also supposed to be the common nodes for P-I and P-II biosynthesis (Shitiz et al., 2015).
The analysis of FA, however, did not show statistically significant variations among selected stolon tissues. Albeit, increased VA content was observed in PKST-3 and PKST-5; stolons having high P-II content, compared to PKST-16 and PKST-18; low P-II content stolons. Keeping in view that FA is a canonical intermediate metabolite in the biosynthesis of guaiacyl-and syringyl lignin's in addition to its proposed role in VA production, the results observed in this study might indicate that flux of FA is directed toward both VA and P-II in PKST-3 and PKST-5, whereas maximum flux of FA might deviate to lignin biosynthesis in PKST-16 and PKST-18. This statement is in agreement with correlation map constructed between P-II content and transcripts abundance values (FPKM) which revealed negative correlation of gene encoding caffeic acid-3-Omethyltransferase (CMT) with P-II content (−0.63) followed by a strong positive correlation of CMT through a PCC value of 0.82 with ferulic acid-5-hydroxylase (F5H), an enzyme that catalyzes the conversion of ferulic acid to 5-hydroxy ferulic acid, which represents as a precursor for syringyl lignin biosynthesis (Wagner et al., 2015). Consequently, it is apparent that up-regulated transcript level of CMT; an enzyme catalyzing the conversion of caffeic acid to FA accompanied by F5H, might be linked to the activation of lignin biosynthesis in PKST-16 and PKST-18. Trabucco et al. (2013) also reported the reduced total lignin content upon down-regulation of gene encoding CMT enzyme in Brachypodium distachyon. As a whole, the analysis of stolon tissues with natural variations for P-II content hypothesized that P-II is produced via shikimate/phenylpropanoid pathway through degradation of FA to VA, which finally integrates with catalpol (CAT) to produce P-II in P. kurroa (Figure 8).
Previous studies have reported that shoots of P. kurroa are independent biosynthetic tissues for P-I, a secondary metabolite produced by the combination of CA and catalpol (Kumar et al., 2016b). Nevertheless, current study established stolons of P. kurroa as autonomous biosynthetic machinery for P-I and P-II, both share a common biosynthetic pathway, it is inspiring to comprehend the negligence of P-II biosynthesis in P. kurroa shoots. To address this, the shoot tissues of four P. kurroa accessions based on their varying P-I content ranged from minimum to maximum (0.13-2.31%) were selected for analysis of p-CA, CA, FA, and VA with the aim to investigate the fate of shikimate/phenylpropanoid pathway in P. kurroa shoots. The analysis of shoot tissues of four P. kurroa accessions for CA and p-CA revealed that increased P-I production limits flux of CA toward p-CA biosynthesis. It is likely since CA shares a common biosynthetic node for P-I and p-CA (Kumar et al., 2016b). In contrast, the analysis of FA revealed slight significant variation between the selected shoots and also showed a high positive correlation by PCC of 0.85 with p-CA. However, C4H and chalcone synthase (CHS) showed a strong positive correlation with PCC of 0.90 compared to PCC of 0.23 between C4H and CMT ( Figure 5D). CHS is an enzyme that catalyzes the formation of naringenin chalcone; the later one serves as a starting metabolite for flavonoids biosynthesis (Zuk et al., 2016). Therefore, this finding underscores that phenylpropanoid pathway might allocate the flux of p-CA to produce both FA and flavonoids in P. kurroa shoots. This statement is in agreement with previous study which revealed that silencing of gene encoding hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyltransferase enzyme (HCT) involved in lignin biosynthesis, re-routed the flux of phenylpropanoid pathway into flavonoids through CHS activity (Besseau et al., 2007). Interestingly, negligible levels of VA were observed in all the selected shoot tissues which implies that P-II biosynthesis arrests in the downstream steps of FA, which possibly limits the VA supply.
Therefore, p-CA, FA, and PA were exogenously applied along with mixture of VA and catalpol in tissue culture conditions not only to establish the course of P-II biosynthesis but also to produce P-II in P. kurroa shoots. The results revealed significant enhancement of VA content in FA fed shoots compared to shoots without treatment. This observation further ascertained that FA is associated with VA production rather than PA which is previously proposed to be the precursor for VA biosynthesis (Kumar et al., 2013;Shitiz et al., 2015). Moreover, significant enhancement of VA content in VA+CAT fed shoots over untreated control might indicate it's accumulation due to exogenous application. This was also supported by the absence of VA content observed in the MS media left after the collection of VA+CAT fed shoots. Upon analysis of P-II content, minor detection was only observed in VA+CAT fed shoots which indicate that exogenous application of VA+CAT, both immediate precursors of P-II biosynthesis, might stimulate its production in P. kurroa shoots. Therefore, different concentrations of VA+CAT were further tested with the aim to observe the progressive increase of P-II content in P. kurroa shoots. Unfortunately, we did not observe significant increment in P-II content among different VA+CAT fed shoots which is possibly due to low activity of a probable acyltransferase that catalyzes the esterification of VA and CAT to produce P-II. However, we have not observed the activity of probable acyltransferase as it is not discovered yet. It is striking that P-II was not taken up by the P. kurroa plants upon its exogenous introduction in purified form. This was clearly demonstrated from the negligible level of P-II observed in shoots fed with P-II whereas its substantial content was detected in the MS media left after the collection of P-II fed shoots. This finding suggests that P. kurroa plants are unable to take P-II in its ready form and subsequently they are dependent on its endogenous production.
Taken together, the findings of this study shed light on metabolic modulations underlying the biosynthesis of P-I and P-II in P. kurroa. It offers convincing evidences that independent mechanisms control the biosynthesis of P-I and P-II in both compartments of P. kurroa, i.e., shoots and stolons. The intermediates of shikimate/phenylpropanoid pathway produced in stolons guided the flux toward P-II biosynthesis via degradation of FA to produce VA, the immediate precursor of P-II. In contrast, flux of precursors deviate from the non-canonical view of P-II biosynthesis in shoot tissues of P. kurroa. Theoretically, we have finally illustrated the basis of P-II negligence in shoot compartment that provides impetus for the future investigation of an enzyme catalyzing the conversion of VA and CAT to produce P-II, a task that can rewrite the P-II production in shoots of P. kurroa.

AUTHOR CONTRIBUTIONS
VK conducted tissue culture experiments. VK conducted molecular experiments and HPLC analysis. AB performed transcriptomes data analysis and VK conducted correlation analysis. VK and RC conceived and designed research. VK, AB, and RC analyzed data. VK and RC wrote the manuscript. All authors read and approved the manuscript.

ACKNOWLEDGMENTS
We are also thankful to the Himalayan Forest Research Institute (HFRI), Panthaghati, Shimla, India for providing authenticated plant material of P. kurroa and Dr. Hemant Sood, Assistant Professor, JUIT, for providing in vitro grown P. kurroa plants for culturing. A special note of thanks to Mr. Kamlesh Kumar, Lab attendant, JUIT for his assistance in plant tissue culture experiments.