Comparative Transcriptome Analysis Reveals Coordinated Transcriptional Regulation of Central and Secondary Metabolism in the Trichomes of Cannabis Cultivars

2.2. Transcriptional Regulation of the Central Metabolism in Trichomes for Enhancing Carbon Feedstock Production

Differential gene expression statistics were then mapped to the metabolic pathways to examine the likely effects of their transcriptional regulation on cannabinoid and terpenoid production. Most gene expression encoding enzymes in the central metabolism were up-regulated, as exemplified by the disproportionately larger number of positive M-values compared to negative ones (Figure 2a). Upon further analysis, the glycolytic pathway in the trichomes was found to be re-routed in a consequential manner compared to the stem, due to a significant increase in the expression of specific genes (Figure 2b). For example, the genome-wide expression of genes encoding the glucose-6-phosphate isomerase (GPI) enzyme increased by 73% overall (D = 2.09 × 10⁴, Table S6), thus allowing more glucose to be channeled into the pathway. In addition, a 154% net increase in the gene expression for fructose-bisphosphate aldolase (ALDO, D = 3.71 × 10⁵, Table S6) implied a markedly enhanced availability of G3P metabolites to the MEP pathway for cannabinoid and terpene production. G3P can also be used to produce more NADPH via higher expression of the gene encoding the glyceraldehyde-3-phosphate dehydrogenase [NADP+] (GAPN) enzyme (log₂FC = 2, D = 1.23 × 10⁴) to meet the increased demand of the cofactors during secondary metabolite production (Figure 2b). The enhanced expression of GAPN in the trichomes has been verified by qPCR (Figure S3). On the other hand, the alternative production of NADH and ATP using G3P, through the action of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and phosphoglycerate kinase (PGK), was found to be suppressed by the down-regulation of the PGK-encoding gene in the trichomes (Figure 2b). This suggests a shift from NADH and ATP production to that of NADPH in the trichomes to support secondary metabolism.

Besides G3P and NADPH, pyruvate formation was also found to be promoted in the trichomes. There was a significantly higher expression of genes for pyruvate kinase (PK) enzymes (Figure 2b) that produce pyruvate, with an overall 78% increase in expression (D = 7.69 × 10⁴, Table S6). In further coordination, the phosphoenolpyruvate carboxylase (PEPC) enzyme was also down-regulated (log₂FC = −1.9, D = 1.99 × 10⁴) to redirect phosphoenolpyruvate to pyruvate conversion, instead of oxaloacetate production. Pyruvate is also used as a substrate for forming acetyl coenzyme A (AcCoA) via pyruvate dehydrogenase (PDH) reactions [24]. Consistently, gene expression encoding a component of PDH (dihydrolipoyllysine-residue acetyltransferase component 5) was found to be up-regulated in the trichomes, with a corresponding log₂ FC value of 2.2, representing an increase of 2.71 × 10⁴ in the ‘Trimmed Means of M-values’ (TMM) expression (Figure 2b). In turn, the AcCoA was in high stoichiometric demand by the MVA pathway and for the biosynthesis of malonyl coenzyme A and hexanoyl coenzyme A, all of which are needed to produce cannabinoids and terpenes. Further, AcCoA was found to be essential for replenishing the tricarboxylic acid (TCA) cycle in the trichomes, as evidenced by the up-regulation of the ATP-citrate synthase (CS) enzyme in comparison to stem tissue (log₂FC = 2.9, D = 6.64 × 10⁴, Figure 2b). Altogether, the importance of enhancing and tuning the balance between four intracellular feedstocks—NADPH, G3P, pyruvate, and AcCoA—in the central metabolism for promoting the production of secondary metabolites in the trichomes was uncovered.

2.3. An Uncovered Pathway for Meeting NAD(P)H Demand in Trichomes

Besides GAPN, the current analysis further implicated the transcriptional amplification of four consecutive reactions catalyzed by malate synthase (MS), malic enzyme [NADP+ dependent] (ME [NADP+]), pyruvate decarboxylase 1 (PDC1), and aldehyde dehydrogenases (ALDH) to further meet the demand of NAD(P)H cofactors for secondary metabolite production in the trichomes (Figure 2c). The pathway is as follows: the MS enzyme first transfers the acetyl moiety of AcCoA to glyoxylate for the formation of malate. In turn, malate, which is also well replenished by the TCA cycle due to enhanced CS expression (Figure S3), is oxidized by ME (NADP+) in the second reaction to form pyruvate, as well as the first NADPH molecule produced by the pathway. In the third reaction, pyruvate, which is also well provided by glycolysis, is then oxidized by PDC1 to form acetaldehyde. The resulting aldehyde molecule can then be oxidized by ALDHs (ALDH2C4 (ALDH family 2 member C4) and ALDH3F1 (ALDH family 3 member F1)) to generate the second NADH molecule as well as acetate (Figure 2c). Besides being major cellular sources of reducing equivalents through NAD(P)H production [25], ALDHs are also required for generating hexanoic acid, a key precursor for cannabinoid production [4]. Importantly, the pathway prioritized NADPH production over NADH, which was similarly observed with the relative regulation of GAPN, GAPDH, and PGK, as described previously. For example, while the fold-change value of the pathway enzymes was generally large (log₂FC = 3.8 to 8.3) compared to stem tissue, the ME (NADP+) enzyme stood out in having an increased expression (D) that was at least 50 times those of the three individual pathway enzymes (1.2 × 10⁶ versus 1.7 × 10⁴–2.3 × 10⁴) and close to a hundred times that of GAPN (1.23 × 10⁴) (Figure S3). This observation thus implicates ME (NADP+) as another major contributor to NADPH production in trichome cells, besides GAPN.

2.4. Modular and Calibrated Regulation of MEP and MVA Pathways

There are also notable features regarding the regulation of MEP and MVA pathways. Firstly, the regulation of the MEP pathway appeared to be modular and calibrated with the up-regulation of four successive end reactions in an increasing manner, i.e., both log₂FC and D-values increased in the order of 4-diphosphocytidyl-2-C-methyl-d-erythritol kinase (CMK), 2-C-methyl-d-erythritol 2,4-cyclodiphosphate synthase (MCS), 4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthase [ferredoxin] (HDS), and 4-hydroxy-3-methylbut-2-enyl diphosphate reductase (HDR) (Figure 3a and Figure S4). Among them, HDS and HDR are known to be rate-limiting [26], and, as such, their regulation is expected to enhance the sensitivity of tuning the pathway. Figure S3 shows the up-regulated expression of MCS, HDS and HDR in the trichomes of all varieties compared to the stem. In this regard, the current finding provides fresh evidence that they serve as control points in the pathway [27]. On the other hand, in the MVA pathway, only the initial fast reaction catalyzed by acetyl-coenzyme A acetyltransferase 1 (ACAT1) was found to be significantly increased in terms of gene expression (Figure 2b, Figure 3a, and Figure S3), thus suggesting a limited increase in the utilization of the pathway. Overall, the modular and calibrated mode of regulation for both pathways might be due to their broader role in integrating environmental and developmental cues, thus necessitating the careful and effective adjustment of pathway activity by specialized modules [27].

2.5. En Bloc and Multi-Loci Up-Regulation of the Hexanoate Pathway and Cannabinoid Biosynthesis

The regulation of the hexanoate pathway was distinct from the MEP and MVA pathways in two ways: (i) the en bloc transcriptional up-regulation of all seven steps in the pathway, which is frequently achieved via (ii) numerous gene loci (Figure 3a). Such multi-loci regulated enzymes included linoleate lipoxygenase (LOX), ALDH, AAE, and OAC, with their log₂ FC values of gene expression ranging from 2.5 to 7.6 (5 loci), 4.8 to 12 (3 loci), 2.1 to 8.1 (3 loci), and 9.1 to 9.2 (2 loci), respectively. This mode of regulation was also found with enzymes directly involved in cannabinoid biosynthesis (Figure 3a). For example, the aromatic prenyltransferases (APT) and THCA synthase (THCAS) enzyme genes showed respective log₂ FC values ranging between 5.1 and 9.9 (5 loci), and 7.7 and 91 (2 loci). The en bloc and multigene loci mode of regulation shared by both pathways can be ascribed to their more straightforward role in cannabinoid production in the trichomes compared to the MEP and MVA pathways.

2.6. Metabolic Profiles of Cultivars

Furthermore, the gene expressions of the metabolic enzymes were compared among the cultivars to elucidate their relative production capacities. Most genes in the core central metabolism, as well as MEP and MVA pathways, were highly expressed in the WS cultivar, unlike others, suggesting that it has the highest NADPH and carbon feedstock productivities (Figure 4a). Consistently, the WS cultivar also had one of the lowest gene expressions for the PEPC enzyme, thus enabling it to redirect phosphoenolpyruvate to pyruvate feedstock production. Additionally, in the WS cultivar, the low gene expression for PDC and ALDH enzymes that drive NADH production (Figure 2c) may be compensated by high CS expression replenishing the TCA cycle to make the cofactors. Similarly, the high expression of CS and ME (NADP+) in WS may make up for its relatively low GAPN expression in producing NADPH. On the other hand, the HB cultivar seemed to excel in NADPH production via the up-regulation of both GAPN and ME (NADP+) genes, but was partly offset by its lower CS expression level compared to the WS cultivar. The WS cultivar also had broadly higher gene expressions in the hexanoate pathway (Figure 4a), suggesting its larger upstream capacity for cannabinoid production. In this regard, its lower ALDH level could already be sufficient for generating the required hexanoic acid for synthesizing cannabinoids (Figure 3a). Both the WS and CD cultivars also had some of the highest transcript abundances for APT, THCAS, and CBDAS enzymes but were still somewhat matched by the HB cultivar for THCAS and by the GT cultivar for CBDAS. Notably, the TH cultivar showed relatively low gene expression levels across all metabolic pathways (Figure 4a).

Expectedly, the gene expression for five major TPSs, namely, myrcene synthase (MYS), limonene synthase (LS), α-humulene synthase (αHS), nerolidol synthase (NES), and germacrene A synthase (GAS), were significantly up-regulated in the trichomes compared to the stem (Figure 3b and Figure S4). The CD cultivar had the highest total transcript levels for LS (TMM = 1.87 × 10⁶) and NES (6.58 × 10⁴), which were about twice those of the TH cultivar with the lowest expression (TMM = 9.66 × 10⁵ and 3.32 × 10⁴, respectively) (Figure 4b). On the other hand, the WS cultivar had the highest gene expression for MYS, αHS, and GAS.

As no chemical profiling data were available for rationalizing the metabolic outputs based on gene expression, publicly available information was used for the few comparisons that were possible (Table S7). Although THC levels were indeed higher in both WS and CD cultivars compared to others, as would be expected from their THCAS expression (Figure 4a), the levels of cannabinoids were, however, not in line with the relative gene expression of THCAS in WS and CD. CBDA and β-myrcene were also expectedly higher in WS compared to CD, based on respective CBDAS and MYS expressions (Figure 4b). On the other hand, the relative levels of limonene in WS and CD cultivars were inconsistent with the relative expressions of LS in the two cultivars.

2.7. Functional Characterization of New TPSs

Besides previously reported TPSs from Cannabis [28,29,30,31], three new TPSs were identified from the RNA-seq data containing full-length open reading frames (ORFs). They were designated as CsTPS3GT, CsTPS4WS, and CsTPS5TH, to reflect their respective origins in the GT, WS, and TH cultivars. CsTPS4WS and CsTPS5TH had very high expressions in the trichomes of all cultivars, whereas CsTPS3GT showed lower expression in the trichomes of the HB, CD, and WS cultivars compared to the stem (Figure S5). Based on amino acid sequence similarity, CsTPS3GT was 96% similar to CsTPS31 (accession number: QLC36839.1), which was previously characterized as producing an unknown sesquiterpene with an RI of 1916 and a base peak of 93 [29]. Furthermore, CsTPS4WS was identical to CsTPS29 (accession number: QLC36833.1), which was previously identified as a linalool synthase [29], while CsTPS5TH was 98.7% similar to CsTPS35 (accession number: QLC36840.1), a linalool/nerolidol synthase [29].

As expected for TPSs, all three peptides contain the ‘DDXXD’ and ‘NSE/DTE’ motifs, which allow them to participate in divalent metal ion-assisted binding of substrates and cofactors [32]. Furthermore, both CsTPS3GT and CsTPS4WS contained a tandem arginine/tryptophan motif, ‘RR(X8)W’, at the N-terminal region (Figure S6), which could be used for monoterpene cyclization [33]. However, the plastid-targeting peptide, which is unique to mono-TPSs, was absent in CsTPS3GT (Figure S6). Phylogenetic analysis placed CsTPS3GT and CsTPS4WS in the TPS-b sub-family, whereas CsTPS5TH fell under the TPS-g sub-family (Figure 5).

The intracellular localization of each TPS was investigated in N. benthamiana leaf cells. Upon transient expression of YFP-tagged TPSs in the leaf cells, both CsTPS3GT and CsTPS5TH were found to localize to the cytosol (Figure 6), which is consistent with their lack of transit peptide sequences and thus suggests that they are sesqui-TPSs, rather than mono-TPSs (Figure 6). On the other hand, CsTPS4WS was localized to the chloroplasts, indicating it to be a mono-TPS.

4.2. RNA Library, RNA-seq, Pre-Processing, and Quality Control

RNA libraries were constructed according to Illumina’s TruSeq strand-specific protocol, with a median insert size of ~230 bp (~100–350 bp). They were then sequenced using HiSeq 2500@PE100 platform (Illumina, San Diego, CA). Adaptor sequences were removed from the reads using Trimmomatic [61], while the quality of the reads was maintained by trimming using the ‘adaptive quality trimming’ algorithm [62]. The latter strikes a balance between quality score and sufficient read length, which is necessary for unambiguous mapping to the genome. Resulting reads < 36 bp were deemed to cause unspecific mapping and were thus discarded. Furthermore, only reads with a mean Q-score > 25 were retained, which translated to <1% error on average for each base pair. Upon completion, quality-control checks using FastQC did not reveal any issues with the final pre-processed reads.

4.3. Reference Genome Mapping and GC-Bias Correction

SALMON [63] was used to quantify transcript expression from the pre-processed reads in a rapid and memory-efficient manner by coupling quasi-mapping onto a reference transcriptome with a read count inference procedure. In doing so, it considers the specific experiment attributes and corrects their biases, including GC contents, to provide a more accurate estimate of the read counts. The latest cs10 transcriptomic and genomic sequences (https://www.ncbi.nlm.nih.gov/genome/?term=txid3482[orgn], accessed on 18 October 2020) were used as references for transcriptome mapping and identification of genomic decoys to prevent spurious mapping, respectively. A read count table based on transcripts was produced as a result.

4.4. Gene-Level Expression, Normalization, and Quality Control

Tximport [64] was used to sum up transcript read counts at the gene level. Following this, read counts were normalized using the TMM method and further divided by gene length to enable comparison across samples [65]. The TMM method normalizes the total RNA amount across samples and assumes that most genes are not differentially expressed. As such, it allows for the comparison of diverse samples, such as from different batches and tissue types, in contrast to within-sample normalization methods, such as RPM, FPKM, and TPM. Principal component analysis was then conducted on the gene expression dataset to assess whether there were any trichome sample outliers.

4.5. Differential Expression Analysis

The NOISeq algorithm [66] was used for the identification of DEGs, as it allows for a good control of false positives among lowly expressed genes by considering both M- and D-values. This is because lowly expressed genes tend to have high absolute M-values due to noise and may be construed wrongly as differentially expressed; however, by considering their low D-values, the false positives can be picked out statistically in a robust manner. This is especially applicable in discerning DEGs in the secondary pathways of plants (e.g., Cannabis plants) that may be lowly expressed but are biologically relevant. An 80% probability of being a true positive is used as the threshold for flagging a gene as differentially expressed. For genes having multiple genomic loci contributing to their expression, the percentage increase in trichome expression over the stem tissue for each gene can be determined by first taking the expression value of each locus in the trichomes as the average of all 5 cultivars. Then, the sum of all loci values can be taken to be the expression value of the gene. Gene expression in the stem tissue was similarly evaluated and then subtracted away from the corresponding expression in the trichome to determine the increase in expression. Gene expression increase was then reported as a percentage of its expression in the stem.

4.6. Hierarchical Clustering and Venn Diagram Analysis of DEG

Hierarchical clustering was carried out in R using complete linkage, and 1 − ρ² as a measure of dissimilarity between all pairwise genes, where ρ is the Pearson’s correlation between the log₂ fold-change (FC) values of pairwise genes. For representation of the clustering in heatmaps, the log₂ FC values of each gene were centered and scaled based on the cultivars. The Venn diagram analysis was conducted online (http://www.interactivenn.net/) [67].

4.7. Gene Mapping to Metabolic Pathways

To analyze gene expressions in the metabolic pathways, they were first identified by a BLASTX search of cs10 gene sequences against the NCBI database to identify their Viridiplantae protein homologs. The KEGG IDs of the latter were then mapped to the pathways using the Blast2GO [68] module in OmicsBox 1.3.11 software (https://www.biobam.com/omicsbox [BioBam Bioinformatics, Valencia, Spain], accessed on 12 November 2020). In this way, cs10 genes with membership in metabolic pathways were identified (Table S5). More were manually picked out by inspecting their cs10 descriptions.

4.8. Quantitative Real-Time PCR (qRT-PCR)

The quantity and quality of RNA were measured using a Nanodrop spectrophotometer (ND-1000, Thermo Fisher Scientific, Waltham, MA, USA), and approximately 500 ng of RNA was reverse transcribed to cDNA using an iScript^TM cDNA Synthesis kit obtained from Bio-Rad, Singapore. Expression levels of enzyme genes along the trichomes of five cultivars and the stem of the HB cultivar were analyzed using qRT-PCR. Primers for qRT-PCR were designed by exploiting the cDNA sequences obtained from the RNA-seq data. The qRT-PCR reactions were performed in a 384-well PCR plate using KAPA SYBR fast master mix (Roche, Singapore) and an ABI PRISM 900HT real-time PCR system. For a total PCR reaction of 5 µL, 0.3 µL of cDNA was used and the cycling profile was set at 50 °C for 2 min, 95 °C for 10 min, 40 cycles of 95 °C for 15 s, and 60 °C for 60 s. After thermal cycles, the dissociation analysis (melting curve) was carried out to confirm the specific amplifications of the PCR reaction by adding a profile of 95 °C for 15 s, 60 °C for 15 s, and 95 °C for 15 s. In the current study, CsGAPDH was used as an internal reference, due to its similar expression across the five cultivars. A non-template reference was included for each gene to eliminate the possibility of random genomic DNA contamination and primer dimer formation. SDS 2.4 software (Applied Biosystems) was used to analyze the obtained results. The threshold cycle (Ct) value of a gene is the cycle number required for the SYBR Green fluorescence signal to reach the threshold level during the exponential phase for detecting the amount of accumulated nucleic acid [69]. Comparative delta Ct values of target genes to CsGAPDH were taken as relative expressions among different tissues. The amounts of target genes, normalized to the CsGAPDH gene, were calculated by 2^{−(Ct[target gene]-Ct[GAPDH])}. Error bars represent means ± SDs. All primers used in this study were designed manually and are listed in Table S1.

4.9. Phylogenetic Tree and Clustal Analysis

A phylogenetic tree was constructed using MEGA7 software (Version 7.0, Pennsylvania State University, PA, USA) [70] by the neighbour-joining method with bootstrap values of 1000 replicates. The required sequences were obtained from the NCBI database, with their accession numbers listed in Table S2. The deduced amino acid sequences of CsTPSs were aligned with other functionally characterized CsTPSs using Clustal W, with the following parameters: gap open: 10; gap extension: 0.1; protein weight matrix: Gonnet; penalties: on; gap separation: 4; cut off: 30%.

4.10. Subcellular Localization of TPSs

Full-length ORFs of CsTPS3GT, CsTPS4WS, and CsTPS5TH were cloned into pENTR vectors using the pENTR™/D-TOPO^® Cloning Kit (Thermo Fisher Scientific, Singapore) and transformed into XL1-blue competent cells. Plasmids from the positive clones were isolated and cloned into the destination vector pBA-DC-YFP [71], which contains the YFP in frame at the C-terminal and cauliflower mosaic virus (CaMV) 35S promoter, to generate CsTPS3GT-YFP, CsTPS4WS-YFP, and CsTPS5TT-YFP, respectively. The constructs were transformed into the Agrobacterium tumefaciens EHA105 strain by a heat shock method [72]. The transformed EHA105 cells were cultured at 28 °C and resuspended in a solution containing 100 µM acetosyringone, 10 mM MES (pH 5.6), and 10 mM MgCl₂. The above mixture was incubated at room temperature for 3 h and later infiltrated into N. benthamiana leaves using a 1 mL syringe. After 2 d, the infiltrated leaves were excised, and the fluorescence signals were observed under a confocal scanning laser microscope (LSM 5 Exciter, ZEISS, Jena, Germany). All constructs were verified by DNA sequencing.

We thank the Temasek Life Sciences Laboratory’s central facility for support in confocal microscopy and Xiujing He for her technical assistance with regard to RNA-seq.