Metabolomic analysis of methyl jasmonate treatment on phytocannabinoid production in Cannabis sativa

2.2. Reagents and standards

Reagent and analytical standard preparation followed (Welling et al., 2021). For identification purposes, olivetolic acid (OA) was purchased from Novachem Pty Ltd (Heidelberg West, Vic., Australia), while JA, 12-oxo-Phytodienoic acid (OPDA), and Jasmone were purchased from Sapphire Bioscience Pty Ltd (Redfern, NSW, Australia). MeJA was purchased from Sigma-Aldrich (St. Louis, MO, USA). Δ9-THC-d3 (2 µg/mL) purchased from Novachem Pty Ltd (Heidelberg West, Vic., Australia) was used as an internal standard.

2.3. Methyl jasmonate (MeJA) and neomycin treatments

Fourteen d after transfer to a short-d photoperiod and on the onset of terminal flowering, defined here as meristem termination with an apical female flower, plant clones were subject to either one of five treatments: 0 mM MeJA (control), 1 mM MeJA, 7.5 mM MeJA, 15 mM MeJA, or 0.5 mM neomycin (negative control) dissolved in RO water supplemented with 0.8% (v/v) EtOH and 0.1% (v/v) Tween 20. Solutions were applied as a whole-plant foliar spray. Three biological replicates (plant clones) were used per treatment, and each biological replicate received 25 mL of the spray solution. Treatments were applied in a separate growing area to prevent cross contamination, and these were administered on three occasions over a 14-d period, with plant material harvested 24 h after the final treatment or 29-d post-exposure to a short-d photoperiod. Plant material taken from the apical (top 30 cm) inflorescence was snap-frozen in liquid nitrogen and stored at -80 °C. Foliage leaves of the 1^st – 3^rd order were manually removed prior to collection to enrich for trichome-containing vegetative tissues and reproductive organs: bracts and flowers, respectively. The complexity and compactness of C. sativa inflorescence architecture complicated the development of a reproducible wounding treatment and was therefore omitted from our experiment.

2.4. Sample preparation for mass spectrometry analysis

Inflorescence material enriched for bracts and flowers (~6 g) was freeze-dried for 7 d. After drying, the material was manually passed through a fine-mesh (250 μm) sieve to remove stems. Floral material was homogenized using a Geno/Grinder 2010 at 1500 rpm for 2 × 30 s intervals within 15 mL polycarbonate vials containing 2 x stainless steel grinding balls. The finely ground floral tissue (50 mg) was weighed into 5 mL Eppendorf Safe-Lock microcentrifuge tubes and extracted in 5 mL LCMS-grade MeOH containing Δ9-THC-d3 (2 µg/mL; internal standard) by vortex mixing for 1 min, followed by sonication for 20 min at room temperature. Particulate material was allowed to sit for a minimum of 3 min and then filtered using a 0.45 μm polytetrafluoroethylene (PTFE) syringe filter. Sample extracts were diluted 40 x with MeOH containing Δ9-THC-d3 (2 µg/mL) and aliquots of the neat and diluted filtrate were transferred to amber HPLC vials (2 mL) and stored at –80°C prior to analysis.

2.5. Liquid chromatography-high-resolution mass spectrometry analysis

UHPLC-HESI/HRMS analysis was as described (Welling et al., 2021), with minor modifications. Briefly, sample runs were carried out on a Thermo Fisher Vanquish Flex UHPLC system with solvent degasser, quaternary pump, temperature-controlled sampler/auto injector and column compartment, and photodiode array detector (DAD) coupled to an Orbitrap ID-X Tribrid high resolution mass spectrometer (Thermo Fisher Scientific Inc., MA, USA). Chromatographic separation was performed using a Phenomenex Kinetex C18 column, 1.7 μm, 150 mm × 2.1 mm (Phenomenex Australia Pty Ltd, NSW, Australia). The mobile phase A was water with 0.1 % (v/v) formic acid and the mobile phase B was acetonitrile with 0.1 % (v/v) formic acid. The following gradient program was used: 10 % B, 0–2 min; 10 to 40 % B, 2–3 min; 40 % B, 3–5 min; 40 to 80 % B, 5–6 min; 80 % B, 6–9 min; 80 to 90 % B, 9–11 min; 90 to 100 % B, 11–12 min; 100 % B, 12–15 min; 100 to 10 % B, 15–16 min; and 10 % B from 16–20 min.

MS was operated using a HESI interface in positive and negative ion modes. For MS and MS², the orbitrap was the mass analyzer used. Pierce FlexMix calibration solution was used to calibrate the MS prior to data acquisition, and the internal mass calibrant fluoranthene (Easy-IC) was activated for real-time mass calibration during data acquisition. Xcalibur (v4.4) software (Thermo Fisher Scientific Inc., MA, USA) was used for data acquisition. For quantitative analysis, Trace Finder (v4.1) software (Thermo Fisher Scientific Inc., MA, USA) was used on full MS data. Integration of extracted ion peaks and comparison of retention time with reference standards was used for peak identity. Δ9-THC-d3 (2 µg/mL) was used as an internal standard, and concentration of each analyte was determined by interpolation from standard calibration curves generated in MS Excel. Results were expressed as average concentration ± standard deviation (s.d.) from duplicate extraction replicates.

2.6. Compound annotation

Small-molecule identification was carried out using Thermo Scientific Compound Discoverer (v3.3) software. The following parameters were used for compound detection and alignment of retention times: mass tolerance [ppm] = 5; min. Peak Intensity: 50 000; Chromatographic S/N Threshold: 1.5; Ions (positive mode): [2M+H]⁺¹, [2M+K]⁺¹, 2M+Na]⁺¹, [M+H]⁺¹, [M+H-H2O]⁺¹, M+K]⁺¹, [M+Na]⁺¹; Base Ions (positive mode): [M+H]⁺¹; [M+H-H2O]⁺¹. Ions (negative mode): [2M-H]^-1, [M-2H]^-2, M-2H+K]^-1, M-H]^-1, [M-H-H2O]^-1; Base Ions (negative mode): [M-H]^-1; alignment model: adaptive curve; maximum shift [min]: 2; mass tolerance: 5 ppm; min. element counts: C H; max. element counts = C₁₀₀ H₃₀₀ Br₅ Cl₁₀ N₂₀ Na₅ O₅₀ P₁₀ S₁₀. To correct for batch effects and reduce drift, peak areas were quality control (QC)-corrected using a QC sample developed by pooling inflorescence samples across all treatments. These were measured after every 15 injections.

Annotation of putative compounds was conducted using calculated elemental composition and searches with high-resolution accurate mass (HRAM) MSⁿ spectral library mzCloud, as well as compound databases (ChemSpider and Metabolika) and ranking (mzLogic) algorithmic tools, along with fragment ion searching (FISh) to predict in silico fragmentation. The following parameters were used for mzCloud searches: compound classes = all; precursor mass tolerance: 10 ppm; FT fragment mass tolerance: 10 ppm; IT fragment mass tolerance [Da] = 0.4; library = auto processed, reference; post processing = recalibrated; max # results = 10; annotate matching fragments = true. Search parameters for DDA were as follows: identity search: HighChem HighRes; match activation type: true; match activation energy: match with tolerance; activation energy tolerance: 20; apply intensity threshold: true; similarity search: confidence forward; match factor threshold: 60. The following databases were included in ChemSpider searches: Cambridge Structural Database, Carotenoids Database, ChEBI, ChEMBL, DrugBank, KEGG, LipidMAPS, NIST Chemistry WebBook, Phenol-Explorer, PlantCyc, PubMed. Database searches were performed by mass and formula using the following parameters: mass tolerance: 5 ppm, max. # of results per compound: 100; max. # of predicted compositions to be searched per compound: 3; The following parameters were used for mzLogic: FT fragment mass tolerance: 5 ppm; IT fragment mass tolerance: 0.4 Da; max. # compounds: 0; max. # mzCloud similarity results to consider per compound: 10; match factor threshold: 30.

The annotation process corresponded to level 2 of confidence as set out by the COordination of Standards in MetabOlomicS (COSMOS) and the Chemical Analysis Working Group (CAWG) Metabolomics Standards Initiative (MSI) (Sumner et al., 2007). The compound list was restricted to compounds with chromatographic peak areas > 50 000 as well as those with MS² fragmentation spectra using the following filtering parameters: background is false; annot.source has status full match in source predicted compositions; group area is greater than 50 000 in any sample group; norm. area has any value in any file; MS² is equal to DDA for preferred ion. Putatively annotated compounds from this list were used for co-expression network construction and metabolomic comparisons between control and treated plants.

2.7. Gene expression analyses

Snap frozen inflorescence material (~2 g) enriched for bracts and flowers was ground using a Geno/Grinder 2010 at 1500 rpm for 2 × 30 s intervals within 15 mL polycarbonate vials containing 2 x stainless steel grinding balls. The Spectrum Plant Total RNA kit (Sigma-Aldrich) was used to isolate total RNA from homogenized inflorescence material (~100 mg fresh weight), with DNase I on-column digestion performed in accordance with the manufacturer’s instructions. DNA-free RNA was eluted in 30 µL of RNase-free water. RNA concentrations were determined using a Qubit 4 Fluorometer and Qubit RNA BR Assay Kit (Thermo Fisher Scientific, USA), with RNA integrity assessed by gel electrophoresis. RNA was extracted from three biological replicates (plant clones) per treatment (0 mM MeJA (control), 1 mM MeJA, 7.5 mM MeJA, 15 mM MeJA, or 0.5 mM neomycin). Single strand cDNA was synthesized using SuperScript IV VILO Master Mix with ezDNase Enzyme (Thermo Fisher Scientific, USA) and diluted 1:5 to a total RNA equivariant of 10 ng/µL.

Quantitative PCR with threshold cycle (Cq) determination was performed using a fluorescence baseline setting of 0.3 (QuantStudio 5 Real-Time PCR System, Applied Biosystems). A total of 2.5 µL of 8 × diluted cDNA was used for quantitative PCR in a total volume of 10 µL with SYBR Select Master Mix (Applied Biosystems). Data were normalized by using a previously validated reference gene EF1α (Riedel et al., 2014; Guo et al., 2018). PCR efficiencies as well as the sensitivity and robustness of the assay were determined by means of calibration curves, with PCR efficiencies calculated using the formula: PCR efficiency = 10^−1/slope – 1 (Bustin et al., 2009). Calibration curves were generated for each target gene primer pair (Supplementary Table S2). PCR products were amplified using Phusion High-Fidelity PCR Master Mix with HF Buffer (Thermo Fisher Scientific, USA) and amplicons sequence-verified at the Australian Genome Research Facility (AGRF). Data were expressed as 40-ΔCq, where ΔCq is the difference between the Cq of the target gene and reference gene. As such, a 40-ΔCq value equal to 40 represents a transcript amount equal to the transcript abundance for the reference gene.

3.2. MeJA treatments increase PC content in a dose-dependent manner

To determine the impact of MeJA treatment on PC content and composition, 14-d-old inflorescences with 1^st – 3^rd order foliage leaves removed were subject to targeted quantitative PC profiling. The industrial hemp genotype used for this analysis exhibited a chemotype high in the C₅ alkyl side chain PC Δ9-tetrahydrocannabinolic acid A (Δ9-THCA) and C₃ alkyl side chain PC Δ9-tetrahydrocannabivarinic acid (THCVA), as can be seen from the total ion chromatogram and extracted ion chromatograms m/z 331.1904 and m/z 359.2217 (Figure 2A). PCs were predominantly in their acidic forms, with neutral PC dry weight concentrations (w/w) orders of magnitude lower than their carboxylic acid precursors (Supplementary Table S4). Other acidic PC species could also be quantified using the extracted ion chromatogram data in conjunction with known retention times of standards, including cannabichromenic acid (CBCA) and cannabigerolic acid (CBGA) as well as cannabidiolic acid (CBDA) and cannabidivarinic acid (CBDVA), although these were at much lower concentrations than Δ9-THCA or THCVA (Supplementary Table S4).

To understand the overall effect of MeJA on PC content, PC values were combined to give a total dry weight (w/w) in the inflorescences of MeJA and neomycin-treated plants. Significant increases in PC content were found in the 1 mM MeJA (p < 0.01) as well as in the 7.5- and 15-mM MeJA (p < 0.001) treatments when compared with the 0 mM MeJA control (Figure 2B). PC content was found to increase in a dose-dependent manner, with an almost two-fold increase in PC content measured in the 15 mM MeJA treated plants. The dry weight (w/w) content of individual PCs followed a similar pattern for both major and minor acidic PCs, with significant increases in PC content in the 7.5- and 15-mM MeJA treatments (Figure 2C). The exception was CBGA, whose content was significantly higher in both the 1- and 7.5-mM treatments, but not the 15 mM treatment (Figure 2C). No significant differences in PC content were observed for the neomycin-treated group (Figures 2B, C).

3.3. Composition of alkyl side chain analogues is impacted by MeJA

The use of a C. sativa chemotype capable of producing high levels of both C₃ and C₅ alkyl side chain PCs allowed for the examination of the impact of MeJA on alkyl side chain analogue composition. The length of the alkyl side chain of PCs is directed by the availability of acyl-CoA starter units (Luo et al., 2019). The use of hexanoyl-CoA during polyketide assembly results in the formation of C₅ alkyl side chain analogues CBDA and Δ9-THCA, while assembly with butanoyl-CoA would result in the formation of C₃ alkyl side chain analogues CBDVA and THCVA (Figure 3A). Significant increases in the CBDVA : CBDA ratio were observed at the 7.5- and 15-mM MeJA treatments (one-way ANOVA, Dunnett’s test, p < 0.05), and a significant linear trend of increasing C₃:C₅ alkyl ratios with increasing MeJA molar concentrations was also found among the Δ9-THC(V)A-type subclasses (one-way ANOVA, trend test, p < 0.05) (Figure 3B). For the control and neomycin treatments, CBDVA and CBDA levels were close to parity, with increases in CBDVA exceeding those of CBDA at the 1-, 7.5- and 15-mM MeJA treatments (Figure 3C). Levels of THCVA also increased at a faster rate than for Δ9-THCA, with THCVA and Δ9-THCA becoming close to parity at the 7.5- and 15-mM MeJA treatments (Figure 3C).

3.4. Metabolomic profiles differ among MeJA treated plants

To further investigate the impact of MeJA treatment on PC biosynthesis, we compared the metabolite profile of inflorescence samples across treatments using an untargeted metabolomic approach by UHPLC-HESI/HRMS with data-dependent MS² acquisition in both positive and negative ion modes. In the positive ion mode, analysis across all treatments identified 137,196 features that accounted for 6,398 putative compounds distinguishable by retention time and molecular weight. In the negative ion mode 78,852 features were detected, and these accounted for 3,778 putative compounds. After filtering by eliminating low abundant chromatographic peak areas < 50 000 counts and those without MS² fragmentation spectra, 1,224 (positive ion mode) and 769 (negative ion mode) compounds remained in the list (Supplementary Table S5).

Results from the multivariate analysis of the methanol-soluble metabolomes of MeJA and neomycin treatments were comparable in both positive and negative ion modes (Figures 4A, B). Principal component analysis (PCA) showed a clear separation between the 1-, 7.5- and 15-mM MeJA treatments, with the first two principal components accounting for ~ 90% (first principal component) and ~ 3% (second principal component) of the variability in the dataset, respectively (Figure 4A). Biological replicates of the 0 mM MeJA control and neomycin treatments clustered together, indicating that the neomycin treatment had negligible impact on the metabolite profile of these plants. Hierarchical clustering analysis (HCA) of inflorescence samples formed two major clades (Figure 4B), one comprising the 0 mM MeJA control, 1 mM MeJA and neomycin treatment samples and the other the 7.5- and 15-mM MeJA treatments (Figure 4B). This suggests that the inflorescences treated with the high molar MeJA concentrations have distinct methanol-soluble metabolomes compared to the neomycin, 0 mM and 1 mM MeJA samples.

Figure 4

Metabolome effects of MeJA and neomycin treatment on inflorescence. (A) PCA of quality control-corrected metabolite chromatographic peak areas in positive and negative modes. (B) HCA and heatmap of quality control-corrected metabolite chromatographic peak areas in positive and negative modes. G1-G3 represent the major metabolite groups by abundance profile across samples. Red boxes indicate the three main metabolite groups. Median values from duplicate extraction replicates used per biological replicate. Hierarchical clustering, heatmap and PCA performed on 1,224 (positive mode) and 769 (negative mode) putative metabolites. For multivariate analysis, mean-centred normalization and log2 transformation scaling were applied to quality control-corrected metabolite abundances.

Three main groups can be distinguished from the heatmap of putative compounds. Group 1 (top, G1), which accounts for approximately two-thirds of all metabolites evaluated, has high chromatographic peak area values in the 7.5- and 15-mM MeJA treatments, followed by intermediate values in the 1 mM MeJA treatment and low peak area values in the 0 mM MeJA control and neomycin treatments (Figure 4B). Group 2 (middle, G2) which is a smaller group and has the opposite chromatographic peak area values (high values in the 0 mM MeJA control and neomycin treatments and low values in the 7.5-and 15-mM MeJA treatments), and Group 3 (bottom, G3) which has peak area values unrelated to the treatment groups (Figure 4B).

3.5. Metabolite co-expression modules correspond to PC content

To gain an understanding of the cellular responses to exogenous MeJA treatments and its impact on PC biosynthesis, weighted gene co-expression network analysis (WGCNA) was used. Originally developed for high-throughput microarray and RNA-seq analyses, WGCNA can also be applied to metabolomic datasets for meaningful biological interpretation (Pei et al., 2017). Using chromatographic peak area abundances from MS data in positive and negative ion modes, WGCNA was employed to identify correlative relationships between metabolites across MeJA-treated samples by clustering highly correlated metabolites into modules and examining their correlations to chemotypes scored through targeted PC profiling.

Weighted co-expression network construction identified 13 modules (six from positive mode, seven from negative mode), as can be visualized from the hierarchical clustering dendrograms together with assigned module colors in Figure 5A. To identify modules that are significantly associated with PC biosynthesis, module eigenmetabolites (ME) were correlated with PC content (Figure 5B). Module–chemotype correlations identified that the eigenmetabolites of the turquoise modules from positive and negative ion modes showed strong positive correlations with PC content (r = 0.99), while eigenmetabolites of the blue modules from positive (r = -0.97) and negative (r = -0.93) ion modes showed strong negative correlations with PC content (Figure 5B; Supplementary Tables S6, S7). These modules corresponded to Group 1 (turquoise modules) and Group 2 (blue modules) from the HCA analysis (Figures 4B, 5A), with the eigenmetabolites of the turquoise modules having the highest abundance at 7.5-and 15-mM MeJA treated samples while the eigenmetabolites of the blue modules having the lowest abundances at these concentrations. There were significant differences in the enrichment of KEGG pathways between the turquoise and blue modules (Table 1). In the blue negatively correlated modules, “Alanine, aspartate and glutamate metabolism (map00250)”, “D-Glutamine and D-glutamate metabolism (map00471)” and “Arginine and proline metabolism (map00330)” were the most significantly enriched pathways, while in the turquoise positively correlated modules the KEGG pathways “Biosynthesis of plant hormones (map01070)”, “Phenylalanine, tyrosine and tryptophan biosynthesis (map00400)” and “Biosynthesis of secondary metabolites (map01110)” were the most significantly enriched (Table 1).

3.6. Identification of hub metabolites

Consistent with the enrichment of secondary metabolite and plant hormone biosynthesis pathways, metabolites involved in oxylipin, and PC biosynthetic pathways are overrepresented in the turquoise module. Among the metabolites within the positively correlated turquoise modules were the PCs Δ9-THCA, THCVA, CBCA, CBDA, CBDVA, and CBGA, as well as the PC resorcinyl intermediate OA, and these compounds matched the nominal mass (Δ ppm < 5) and retention times of certified PC reference standards (highlighted in black bold in Supplementary Tables S6, S7). Other putative PCs were also identified among these modules, but these were not verified using reference standards. As evidence of the robustness of the module assignments, CBGA was among the metabolites present in the yellow module positively correlated with CBGA content (r = 0.96, p-value 0.002) (Figure 5B). Several jasmonate (JA)-related compounds were also annotated within the turquoise modules by consensus evaluation based on elemental composition prediction, spectral library (mzCloud) and database (ChemSpider) searches. These include JA, MeJA, 12-oxophytodienoic acid (OPDA), and the JA derivative jasmone (highlighted in blue bold in Supplementary Tables S6, S7), and these compounds were confirmed by comparison with certified reference standards. In addition, metabolites that have a high significance for PC content also showed a high module membership (MM) in the turquoise modules for both the positive and negative ion modes (Figure 5C). Taken together, these data suggest that the eigenmetabolite of the turquoise modules demonstrates a strong relationship between the oxylipin and PC biosynthesis pathways.

Metabolites which display a high degree of connectivity in interaction networks are considered to be biologically important (Pei et al., 2017). Among the centrally located intramodular hub metabolites and in the top 10 most connected metabolites in the turquoise module was a putative ω-oxo acid 12-oxo-(10E)-dodecenoic acid, commonly known as traumatin (Figure 5D; Supplementary Table S6). This metabolite had a MM of 0.99 (p-value 7.57E-14) in the turquoise module and is also in the top 10 most significantly correlated compounds to PC content in this module (r 0.99, p-value 2.05E-11) (orange bold in Supplementary Table S6). Analysis of the MS spectra of this peak at 7.715 min supported its annotation as traumatin (orange bold in Supplementary Figure S3). Traumatin is a product of lipoxygenase (LOX) and hydroperoxide lyase (HPL) activity, which, if cleaved from a C₁₈ unsaturated FA, would result in the formation of C₆ aldehydes that are hypothesized to be the progenitors of hexanoyl-CoA in PC biosynthesis (Figure 6A; Supplementary Figure S1). Interestingly, we also identified other putative ω-oxo-acids and structurally related FAs in the turquoise module, including a C₁₂ ω-oxo-acid (5Z,8Z,10E)-12-oxo-5,8,10-dodecatrienoic acid (r 0.98, p-value 2.77E-11) and a C₁₀ oxo-acid 9-oxo-decenoic acid (r 0.98, p-value 1.07E-10) which both had a similar level of MM and correlation to PC content (orange bold in Supplementary Tables S6, S7and Supplementary Figures S4, S5), indicating that multiple putative LOX-HPL cleavage products have a strong relationship with PC content. We also searched for intermediates associated with other metabolic routes capable of producing short-medium chain acyl-CoA PC precursors, including those associated with branched-chain amino acid catabolism & α-ketoacid elongation (e.g., 2-oxo-valeric, 2-oxo-hexanoic acid and 2-oxo-heptanoic acid (Kroumova et al., 1994) (scenario 1 Supplementary Figure S1). However, these metabolites were not found among the modules (Supplementary Tables S6, S7), and we did not detect any compounds that matched the molecular formula and mass of these intermediates from the MS data (data not shown).

The authors would like to acknowledge the academic input of Dr Ricarda Jost and Dr Oliver Berkowitz who designed the OAC and THCAS primers used for qPCR analysis. We would also like to thank Professor Mathew Lewsey and Dr Muluneh Oli for their useful discussions on jasmonates.