The dry weight of the CMCs and DDCs showed variation within days. Taking same dry weight as initial material, the CMCs produced maximum dry mass on 15th day, whereas DDC produced maximum dry weight of cells on 12th day (Figure 1I). Therefore, we chose 12 days old cells to analyze the alkaloid contents. Total alkaloids, vindoline, catharanthine, ajmalicine, and vinblastine were measured in CMCs and DDCs by HPLC. Interestingly, we observed that the total alkaloid content was significantly higher in CMCs (43.72 mg/G) as compared to DDCs (2.59 mg/L) (Figure 1ii). Other indole alkaloids also showed a similar trend i.e., higher concentrations in CMCs as compared to DDCs. Vindoline and vinblastine were only detected in CMCs. CMCs had 18.77-fold and 182.57-fold higher content of catharanthine and ajmalicine as compared to DDCs, respectively (Figure 1iii–vi). These results indicate that the total alkaloid content as well as individual indole alkaloid content is higher in CMCs, thus these cells can be used for commercial alkaloid biosynthesis.

The Illumina HiSeq transcriptome sequencing of four samples of C. roseus cells i.e., CMCs and DDCs resulted in 50,296,033 and 46,706,898 raw and clean reads, respectively. On an average 91.33% clean reads could be mapped 42,641,912 of 46,706,898 clean reads. Overall, this project resulted in 28.02 Gb clean data; the clean data of each sample reached 6 Gb, and the Q30 base percentage was 93% and above (Supplementary Table S1). A total of 120,808 unigenes could be annotated in the seven databases; 96,004 unigenes could be annotated in at least one database (Figure 2A). The FPKM for CMCs was higher than that of DDCs (Figure 2B). Pearson Correlation Coefficient between the CMCs replicates was 0.98% whereas between the DDC replicates it was 0.96% (Figure 2C).

On the basis of the screening conditions i.e., FDR (<0.05) and log2 foldchange (≥1) 9,564 genes were differentially expressed between CMCs and DDCs (Supplementary Table S2; Figure 3A). Of these, 1,184 and 759 genes were exclusively expressed in DDCs and CMCs, respectively. Overall, 3,713 genes had higher expression while 5,851 had lower expression in CMCs as compared to DDCs. The DEGs were enriched in 137 different KEGG pathways (Supplementary Table S3; Figure 3B). The KEGG classification showed that the DEGs were classified in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. Most DEGs were classified in pathways related to metabolism with the highest number of DEGs (1706) classified in metabolic pathways followed by biosynthesis of secondary metabolites (1,008). In the case of environmental information processing the DEGs were classified in plant hormone signal transduction (245) and MAPK-signaling pathway-plant (235) (Supplementary Figure S1). These observations indicate that the studied cells are experiencing large-scale metabolic changes and hormones and MAPK signaling play essential role in their characteristics.

The exclusively expressed genes in the DDCs could be annotated as 522 KEGG terms. The GO annotation classified the genes into three categories i.e., BP (biological process), CC (cellular component), and molecular function (MF) (Figure 3C). Major BP in which the DEGs were classified included cell prolification, reproduction, cellular component or biogenesis, biological regulation, and regulation of biological processes. The DEGs were classified in GO for 10 MFs; the highest number of DEGs were classified as transcription regulator activity and transporter activities.

Since we observed higher content of alkaloids in CMCs as compared to DDCs, we explored the expression of genes that were enriched in isoquinoline alkaloid biosynthesis, Tropane, piperidine and pyridine alkaloid biosynthesis, and Indole alkaloid biosynthesis (as well as the upstream pathways i.e., tryptophan biosynthesis and monoterpenoid biosynthesis). In the indole alkaloid biosynthesis pathway, only six transcripts were differentially regulated which were annotated as strictosidine synthase (Figure 4). Of these, two transcripts i.e., Cluster-11739.14165 and Cluster-11739.372 showed reduced expression in DDCs as compared to CMCs suggesting a possible role of these transcripts in the downstream alkaloids’ biosynthesis (Figure 1).

Nineteen DEGs were enriched in tropane, piperidine, and pyridine alkaloid biosynthesis (Figure 4). Four of these were annotated as aspartate aminotransferase mitochondrial (GOT2); all of these had higher expression in CMCs as compared to DDCs. Two transcripts annotated as tyrosine aminotransferase (TAT) showed increased expression in DDCs as compared to CMCs. Rest of the genes i.e., primary-amine oxidases (PAO) and tropinone reductase I (TR1) showed both up- and down-regulation in both types of cells. However, notably only one tropinone reductase I had higher expression in CMCs while rest of the transcripts had lower expression in CMCs as compared to DDCs. Two transcripts (Cluster-11739.15071 and Cluster-11739.71246) annotated as tyrosine aminotransferase showed lower expression in CMCs as compared to DDCs. A bifunctional aspartate aminotransferase and glutamate/aspartate-prephenate aminotransferase (PAT, Cluster-11739.18216) was exclusively expressed in CMCs as compared to DDCs. The GOT2, PAT, PAO, and TAT were also enriched in isoquinoline alkaloid biosynthesis pathway. Additionally, two transcripts annotated as codeine 3-O-demethylase (CODM, Cluster-2412.0 and Cluster-9877.0) had lower expression in CMCs as compared to DDCs. The reduced expression of GOT2 in DDCs as compared to CMCs indicates that interconversion of l-phenylalanine and phenylpyruvate is lower in DDCs. This same reaction is also catalyzed by TAT; however, these genes showed lower expression in CMCs as compared to DDCs. These observations indicate that this reaction might be differently regulated in both types of cells (Figure 4).

Since tryptophan biosynthesis and monoterpenoid biosynthesis pathways provide the key molecules that are further processed in indole alkaloid biosynthesis, therefore, we explored the DEGs enriched in these pathways. Forty-six DEGs were enriched in tryptophan metabolism pathway. Of these, 10 DEGs annotated as acetyl-CoA C-acetyltransferase (ACAT1), amidase, aromatic aminotransferase (AroA), and dihydrolipoamide dehydrogenase (DLD) had higher expression in CMCs as compared to DDCs, while the rest of the genes had higher expression in DDCs (Figure 4). Finally, nine DEGs were enriched in monoterpenoid biosynthesis pathway (Figure 4).

There are two types of transporters i.e., H+ antiport carrier system and ABC transporters, involved in intra-cellular alkaloid transport mechanisms (N. Shitan and K. Yazaki, 2007). In this regard, we noted the differential expression of ATP-binding cassette subfamilies A (ABC1), B (MDR/TAP), C (CFTR/MRP), E, and G (WHITE). All ABC1s were highly (solely) expressed in DDCs, whereas, the 1 E member and maximum number of B and G family members were highly expressed in CMCs. The C family had four members up and four downregulated in CMCs as compared to DDCs. For H+ antiport system, we found differential expression of 24 Ca²⁺/H⁺ antiporters (18 of which belonged to TMEM165/GDT1 family) (see orange highlighted gene IDs in Supplementary Table S2). These transporters showed mixed expression i.e., some transcripts were upregulated while others were downregulated in CMCs as compared to DDCs.

Five genes related to auxin biosynthesis were differentially expressed; a tryptophan synthase beta chain (Trp b), a tryptophan---pyruvate aminotransferase (TAA), and three anthranilate synthases (AS). Trp b was exclusively expressed in CMCs, whereas TAA had fractional expression (0.055) in DDCs as compared to CMCs (2.45). The three AS’s were highly expressed in DDCs. These expression changes indicate that anthranilate synthesis is higher in DDCs leading to higher auxin biosynthesis in DDCs, which in turn may negatively regulate TIA biosynthesis. As far as cytokinin biosynthesis is concerned, we noted that 35 of 46 cytokinin biosynthesis related DEGs were highly expressed in CMCs, indicating higher cytokinin biosynthesis in CMCs. The DEGs with higher expression were beta-glucosidases (BGLs), cytokinin dehydrogenase (CDH), cytokinin riboside 5′-monophosphate phosphoribohydrolase (LOG8-like), a CS, cytokinin-N-glucosyltransferases (UGT76C2), cis-zeatin O-glucosyltransferase, cytokinin riboside 5′-monophosphate phosphoribohydrolase, and Arabidopsis histidine kinase 2/3/4. Whereas those with lower expression in CMCs were eight BGLs, a cytokinin synthase (CS), and two Arabidopsis histidine kinase 2/3/4s. Two 1-aminocyclopropane-1-carboxylate synthase-like (ACSs) and three aminocyclopropanecarboxylate oxidases (ACO1s) showed higher expressions in CMCs. However, other ACSs (two ACS3s) and ACOs (ACO1 and ACO3) had contrasting expression patterns. Indicating, ethylene biosynthesis is being regulated by different transcripts in the 2 cell types. We also noted that allene oxide synthases (AOC), and phospholipase A1s (PLA1) had higher expression in CMCs as compared to DDCs, implying higher JA biosynthesis in CMCs. Similarly, seven genes related to SA biosynthesis i.e., isochorismate synthases (ICSs), phenylalanine ammonia-lyases (PALs), and enhanced disease susceptibility one protein (EDS1) showed increased expression in CMCs. However, one of the ICSs (Cluster-23336.0) showed DDCs specific expression. Overall, these expression trends indicate higher SA biosynthesis in CMCs (see green highlighted gene IDs in Supplementary Table S2).

A total of 245 and 235 DEGs were enriched in plant hormone signal transduction and MAPK signaling-plant pathways, respectively. Seventy-seven DEGs were common in both pathways. All hormone signaling pathways were differentially regulated in the 2 cell types of C. roseus. Most prominently, we observed the higher expression of AUX1 (auxin transporter protein 1), AHP (Arabidopsis histidine phosphotransfer protein), TF (phytochrome interacting factor 4), SnRK2 (serine/threonine-protein kinase SRK2), EIN3 (ethylene-insensitive protein 3), ERF1/2 (ethylene-responsive transcription factor 1), BKI1 (BRI1 kinase inhibitor 1) and BZR1/2 (brassinosteroid resistant 1/2), and JAZ (jasmonate ZIM domain-containing protein) in DDCs as compared to CMCs. On the other hand, TIR1 (transport inhibitor response 1), PYR/PYL (abscisic acid receptor PYR/PYL family), ETR (ethylene receptor) and EBF1/2 (EIN3-binding F-box protein), and CYCE3 (cyclin D3, plant) had higher expression in CMCs as compared to DDCs. The other DEGs showed varied expression in both types of cells i.e., multiple transcripts were expressed. Some had higher expression while others had lower expression in the same type of cells (Figure 5A). Considering that the Jasmonic acid (JA) signaling controls multiple steps upstream of the monoterpenoid biosynthesis and indole alkaloid biosynthesis, the differential expression of JAZ is important. The lower expression of JAZ in CMCs indicates the higher activity of MYC2s as compared to DDCs. However, four MYC2s were exclusively expressed in DDCs, while four others were upregulated in DDCs. Only four MYC2s had higher expression in CMCs as compared to DDCs (Figure 5A). Thus, this step can’t be solely responsible for the observed changes in the indole alkaloid biosynthesis.

Large-scale changes were observed in the MAPK signaling plant pathway (Figure 5B). Specific to JA, we observed that three transcripts annotated as MMK3 (Cluster-11739.74927, Cluster-11739.87079, and Cluster-11739.10429) were highly expressed in CMCs as compared to DDCs. MMK3 is also present upstream the MYC2, thus its importance can’t be negated in terms of the upstream regulation of indole alkaloid biosynthesis. Together with JAZ, the MMK3 genes’ role should be explored in detail in future studies for their effect on the biosynthesis of indole alkaloid biosynthesis.

In addition to the above-described pathways that were related to alkaloid biosynthesis, the KEGG pathway analysis showed the significant enrichment of DEGs in terpenoid backbone biosynthesis, metabolic pathways, biosynthesis of secondary metabolites, flavonoid biosynthesis, phenylpropanoid biosynthesis, and isoflavonoid biosynthesis (Figure 3B). Two of these pathways i.e., metabolic pathways (1,689) and biosynthesis of secondary metabolites (1,008) had a higher number of DEGs enriched in them. These observations indicate that the two types of C. roseus cells differ in their metabolic states and secondary metabolites contents. In the case of the phenylpropanoid biosynthesis pathway, the gene [(trans-cinnamate 4-monooxygenase (C4H), Cluster-11739.25135)] controlling the conversion of cinnamic acid to p-coumaric acid as well as cinnamoyl-coA to p-Coumaroyl-coA were had higher expression in CMCs as compared to DDCs. Whereas genes annotated as phenylalanine ammonia-lyase (PAL), ferulate-5-hydroxylase (F5H), cinnamyl-alcohol dehydrogenase (CAD), coniferyl-aldehyde dehydrogenase (CALDH), and scopoletin glucosyltransferase (TOGT1) had lower expression in CMCs as compared to DDCs. All other genes showed variable gene expression in both types of cells i.e., some transcripts of the same annotation were expressed higher in one type of cells while other had higher expression on the second type of cells (Supplementary Table S4). The C4H gene was also enriched in flavonoid biosynthesis. Most of the other genes in this pathway showed increased expression in DDCs as compared to CMCs. These observations indicate that CMCs have lower phenylpropanoid and flavonoid biosynthesis as compared to DDCs.

The terpenoid backbone biosynthesis pathway, which is present upstream of the monoterpenoid biosynthesis pathway, was also differentially regulated between CMCs and DDCs. Particularly, 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR), (E)-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase (ISPG), and diphosphomevalonate decarboxylase (MVD) had higher expression in CMCs as compared to DDCs. Whereas several other intermediate steps of the terpenoid backbone biosynthesis i.e., between the glycolysis and monoterpenoid biosynthesis were either upregulated in DDCs as compared to CMCs and/or multiple transcripts showed variable expression in both types of C. roseus cells (Supplementary Table S4). Since these steps are present upstream of the indole alkaloid biosynthesis pathway, therefore, the genes can be important targets for desired alkaloid production in CMCs.

The comparative transcriptome analysis of CMCs and DDCs indicated differential expression of 529 and 162 transcription factors (TFs) and transcriptional regulators (TRs), respectively. The TFs were classified in 50 different families; 148 and 381 TFs showed increased and decreased expression in CMCs as compared to DDCs, respectively. Of the TRs, 88 and 74 had higher and lower expression in CMCs as compared to DDCs, respectively; the TRs were classified in 19 families of which the most prevalent were SNF2s. The TFs that were exclusively expressed in CMCs were classified as Alfin-like, AP2/ERF, B3, ARF, bHLH, C2H2, C3H, FAR1, HD-ZIP, LOB, MADS-M-type, MYB-related, and NAC; overall 28 TFs were specific to CMCs. The most important observation was the higher expression of CPP, E2F-DP, HB-WOX, OFP, SRS, trihelix, TUB, and whirly TFs in CMCs as compared to DDCs. On the contrary, most other TFs found in our transcriptome sequencing were highly expressed in DDCs as compared to CMCs (Supplementary Table S5). These expression trends indicate that these TFs and TRs can be associated with the changed TIA concentrations in both cell types.

We studied the expression of 12 key genes involved in the biosynthesis of indole alkaloid biosynthesis. The qRT-PCR results showed that the expression of all the studied genes was higher in CMCs as compared to DDCs, except for STR gene (Figure 6). These results are consistent with the higher concentrations of total alkaloids as well as the individual alkaloid i.e., vindoline, catharanthine, ajmalicine, and vinblastine (Figure 2).

The dedifferentiated cells [(DDCs) (CAT-1)] of C. roseus were established according to the modified method of Zhou et al., (Zhou et al., 2015). Whereas, the cambial meristematic cells [(CMCs) (CAT-2] are re-induced establishments of vinblastine containing CMCs reported in Zhou et al., (Zhou et al., 2015) (Figure 7). Briefly, twigs were collected from C. roseus growing in wild, rinsed in the running tap water for 30 min, disinfection with 75% ethanol for 30 s, and rinsing five times with sterilized distilled water (dH₂O). The disinfected twigs were then put in 0.05% HgCl₂ for 6 min, rinsing with dH₂O five times, followed by rinsing with 150 mg/L citric acid (dissolved in dH₂O). Cambium, phloem, cortex, and epidermal tissues were peeled from the xylem, and the epidermal tissue was placed on pH-adjusted (5.75) Murashige and Skoog (MS) medium (Murasnige and Skoog, 1962) for 24 days. The MS medium was supplemented with 1.0 mg/L a-naphthalene acetic acid (NAA), 30 g/L sucrose, and 4 g/L gelrite. After 24 days, the cambial cells formed a flat plate while phloem, cortex, and epidermis originated DDCs proliferated in irregular fashion. Cambial cells were separated, transferred to 250 ml Erlenmeyer flasks containing 100 ml of MS medium supplemented with 2.0 mg/L NAA, 10 g/L sucrose, and 4 g/L gelrite. The CMCs were cultured at 25°C under continuous dark and were maintained by serial subculture every 12 days. Two months prior to the experiments, CMCs were transferred to 100 ml of MS solid medium in a 250-ml Erlenmeyer flask and grown at 25°C with a 12/12-h light/dark photoperiod. Since cells were sub-cultured for every 12th day, therefore after 6 months, the cells used for the alkaloid analyses were 13th generation.

For determining growth rates, CMCs and DDCs were subcultured every 12 days. CMCs or DDCs (5.0 g fresh weight) were transferred to 250-ml Erlenmeyer flasks respectively, containing 100-ml MS liquid medium supplemented with 2.0 mg/L NAA and 30 g/L sucrose. Suspension cultures of C. roseus CMCs and DDCs were incubated on a HZT-2 gyrotory shaker (Donglian Electronic and Technol. Dev. Co., Beijing, China) at 25°C and 120 rpm under continuous light. Cultures were harvested after 3, 6, 9, 12, 15, 18 days. Cell growth was determined by grams of dry weight (DW) per liter as reported earlier (Zhou et al., 2015).

The total alkaloid content was estimated using HPLC by following the method described by Harborne (Harborne, 1998). We also measure the concentrations of vindoline, catharanthine, ajmalicine, and vinblastine in the CMCs and DDCs. For this, first we isolated these compounds following the method described by Schripsema and Verpoorte (Schripsema and Verpoorte, 1992). The extracts were filtered and analyzed by HPLC. The separation of the four compounds was achieved by using Diamonsil C18 column (250 mm × 4.6 mm, 5 μm), by using 380 ml solution A [water: diethylamide = 986:14 (v/v), pH 7.2] combined with 620 ml solution B [methanol: acetonitrile = 4:1 (v/v)] as the mobile phase. The detection was performed at a flow rate of 1.5 ml min-1 and a wavelength of 220 nm within 40 min as described earlier (Hashemabadi et al., 2018).

Before analyses, the raw sequencing data was checked for quality control by removing reads (having adaptors, paired reads having >10% N content, and if the low-quality basis Q ≤ 20 was >50%). We then determined the error distribution and GC content in the sequencing reads. We then used BLAST (Johnson et al., 2008) to compare the unigene sequences with KEGG (Kanehisa and Goto, 2000), NR (Deng et al., 2006), Swiss-Prot (Apweiler, 2001), GO (Ashburner et al., 2000), COG/KOG (Koonin et al., 2004), Trembl databases (Apweiler et al., 2004). Furthermore, we predicted unigenes’ amino acid sequences and used HMMER software to compare the sequences with Pfam (Bateman et al., 2002). The gene expression was quantified by splicing the transcripts by Trinity (ref. sequence) and then map the clean reads to each ref. sequence by using bowtie2 (Langmead and Salzberg, 2012) in RSEM (Li and Dewey, 2011). The Fragments Per Kilobase of transcripts per Million fragments mapped (FPKM) was calculated and visualized in R. Pearson Correlation Coefficient (PCC) and Principal Component Analysis (PCA) were done on the gene expression data in R. The differentially gene screening between the CMC and DDC was done in DESeq2 (Varet et al., 2016). Then Benjamini–Hochberg method (Benjamini and Hochberg, 1995) was used to perform hypothesis test correction on p-value to obtain false discovery rate (FDR) and screen the differentially expressed genes (DEGs) with the conditions that if FDR (<0.05) and log2 foldchange (≥1), the gene was differentially expressed. After identifying DEGs, they were enriched in KEGG pathways in KOBAS2.0; care was taken to reduce false positives for the prediction of KEGG pathways by considering FDR <0.05. Finally, we used iTAK software to predict plant transcription factors (TF); it integrates PlnTFDB and PlantTFDB and uses TF family and identifies TF through HMM-HMM scan comparison (Zheng et al., 2016).

We studied the expression of 12 genes that are associated with the TIA biosynthesis pathway (Table 1). First strand cDNAs were synthesized from 100 ng of total RNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystem, United States). The primers for each gene are given in Table 1. The reactions were carried out in 20 μl each on a Rotor-Gene 6,000 machine (Qiagen, Shanghai, China). 40S Ribosomal protein S9 (RPS9) was used as a single reference gene (Menke et al., 1999). The reaction mixture contained the following amounts of the reagents. SYBR ^® Premix Ex Taq (Tli RNaseH Plus) (2x) 10. µL, forward primer (5 Mm) 0.8 µl, reverse primer (5 Mm) 0.8 µl, DNA template 2.0 µl, RNase Free dH2O 6.4 µl. The reactions conditions were as follows. Pre-denaturation at 95°C for 30 s, 20°C/s (1 cycle), 40 cycles of 95°C for 5s, 60°C for 20s and 20°C/s, followed by the melting curve analysis at 95°C for 0 s, 20°C/s, 65°C for 15s, 20°C/s, 95°C for 0 s, and 0.1°C/s. For each sample, two biological replicates were analyzed in independent runs.