Charting the Cannabis plant chemical space with computational metabolomics

Metabolite extraction and standard preparation

Fifty milligrams (50 mg) of the plant materials were weighed and extracted in 1 mL of 80% methanol (Chemlab, UK). The samples (leaves and flowers per respective cultivar) were spun overnight in a digital rotisserie tube rotator at 70 rpm. The crude extracts were then centrifuged at 5678× g in a benchtop fixed-angle centrifuge (Thermo Fisher, Johannesburg, South Africa). After centrifugation, the supernatants were filtered using 0.22 μm nylon filters into glass vials with 500 µL inserts. In this study, seven independent replicates for each sample group were weighed and prepared. The prepared samples were then stored at 4 °C until analysis. For chemical baiting, the vitexin standard was prepared in 50% methanol to a concentration of 5 ppm. The standard was treated in the same manner as the samples for experimental procedures.

Data acquisition – LC-MS/MS analyses of crude extracts

Samples (methanol extracts) were analyzed on a liquid chromatography–quadrupole time-of-flight mass spectrometry (LC-qToF-MS) instrument (LCMS-9030, Shimadzu Corporation, Kyoto, Japan). For chromatographic separation, a sample volume of 3 µL was injected on a Shim-pack C18 column (100 mm × 2.1 mm, 2.7 μm) (Shimadzu Corporation, Kyoto, Japan) thermostatted at 55 °C. In addition to the stationary reverse phase, the chromatography was carried out with a binary mobile phase, applying a gradient elution method. The solvent system comprised solvent A consisting of 0.1% formic acid in Milli-Q water (both HPLC grade, Merck, Darmstadt, Germany) and solvent B being methanol (UHPLC grade, Romil SpS, Cambridge, UK) with 0.1% formic acid, with a flow rate of 0.4 mL/min. The gradient elution was performed as follows, B referring to organic composition (i.e., solvent B): 5% B for 3 min, 5–40% B over 3–5 min, 40–95% B over 5–12 min, then 95% B for from 12 min to 18 min. The gradient was changed back to initial 5% B at 18 min and kept to 20 min. The column was re-equilibrated for 3 min.

The effluents from chromatographic separation were further analyzed with the high-definition mass spectrometer, equipped with electrospray ionization (ESI), acquiring both negative and positive spectral data. The mass spectrometer parameters used were the following: 4.5 kV interface voltage, with interface temperature of 300 °C; 3 L/min flow rate for nebulization and dry gas; DL temperature of 250 °C, and 400 °C for heat block; detector was operated at 1.8 kV voltage. For monitoring the accuracy of acquired mass-to-charge ratio (m/z), sodium iodide (NaI) was used as a calibration solution. Both non-fragmented (MS1) and fragmented (MS2) spectral data were acquired with m/z range of 100–1200 Da. For MS/MS experiments, data-dependent acquisition method was applied, with an intensity threshold of 5000 counts. The collision-induced dissociation (CID) method was applied for the fragmentation of ions, using argon as a collision gas at a collision energy of 30 eV with a spread of 5 eV.

Data mining: processing and chemometrics analyses

The acquired raw datasets (ESI negative and positive centroid data) were processed using MetaboAnalyst version 5.0 (Pang et al., 2021). The data was processed with UPLC-Q/TOF default parameters using the centWave method. From the default setting, the max peak width was changed to 15. The resulting feature tables with 3560 and 3570 features for (ESI) negative and positive datasets respectively were imported into soft independent modelling of class analogy (SIMCA) version 17.0 software (Sartorius, South Africa) for chemometrics modelling. Principal component analysis (PCA) and hierarchical cluster analysis (HCA) models were computed for data exploration, and group classification and discriminant analysis.

Molecular networking and metabolite annotation

Prior to computing the feature-based molecular networks (FBMN) through the Global Natural Product Social (GNPS), the acquired raw datasets obtained from the Shimadzu LCMS-9030-qTOF-MS were converted to open-source format (.mzML) files. The mzML files were then uploaded to Mass Spectrometry-Data Independent AnaLysis (MS-DIAL) platform for data processing. The parameters used for MS-DIAL data processing parameters included mass accuracy MS1 and MS2 tolerance of 0.25 Da and 0.1 Da respectively, with the MS/MS range of 50–800 Da; the minimum peak height was 2000 amplitude, mass slice width of 0.1 Da for peak detection, a sigma window value of 0.5 was used, and retention time tolerance was 0.1 min and MS1 tolerance set at 0.015.

Post MS-DIAL data processing, the GNPS export files, both GNPS MGF files and feature quantification tables, were exported into the GNPS ecosystem using the WinSCP server for molecular networking. FBMNs were then computed and generated using the GNPS molecular networking workflow. The parameters used for computing FBMN for the leaves and flower spectral datasets included precursor ion mass tolerance of 0.25 Da with fragment ion mass tolerance of 0.25 Da. For spectral similarity, a cosine score cut off was 0.65, with a minimum of 4 matched fragment ions. For annotation, various spectral libraries were searched and these included MassBank, ReSpect and NIST. To visualize and analyze the computed networks, the Cytoscape network visualization software (version 3.8.2) (Shannon et al., 2003; Smoot et al., 2011) was used.

All semi-annotated and some unmatched/unidentified nodes, visualized in Cytoscape, were verified using their empirical formulae calculated from accurate mass and fragmentation patterns. In addition to the manual inspection of annotations, some natural products dereplication databases such PubChem (https://pubchem.ncbi.nlm.nih.gov) and Dictionary of Natural Products (http://dnp.chemnetbase.com/faces/chemical/ChemicalSearch.xhtml) were also used. The annotations were also verified and confirmed against available literature. All metabolite annotations were carried out at level 2 and 3 of the Metabolomics Standards Initiative (MSI) (Sumner et al., 2007). All validated annotations (manual and computational) are listed in supplementary Table S1. The networks generated were further explored using network annotation propagation (NAP), where first in silico fragmentation was applied to the generated FBMN where structural searches were performed in databases e.g., GNPS, CHEBI, SUPNAT, and DRUGBANK.

To perform network annotation propagation (NAP) (da Silva et al., 2018) jobs, both the fusion and consensus scores were used based on the first 10 structural candidates. Using scoring methods, NAP re-ranks the candidate structure list based on the network topology (Ernst et al., 2019; Kang et al., 2019; Nephali et al., 2022). The two scoring methods, as previously mentioned, utilized by NAP include (a) fusion scoring, which uses MetFrag in silico prediction with MetFusion based on spectral library matches within a molecular family; and (b) consensus scoring, which uses the structural similarity from in silico candidates across the spectral nodes of a molecular family. The FBMNs were also explored using DEREPLICATOR (Mohimani et al., 2018) for peptidic structural annotations.

Substructure annotation and initial exploration was carried out using the MS2LDA (MS2 latent Dirichlet allocation) tool (van der Hooft et al., 2016; Rogers et al., 2019) in GNPS and substructure annotations from MotifDB were included in the analyses. From the default settings, some of the parameters were changed as follows for the ESI negative dataset: Bin width was set at 0.01 (Tof data) and LDA free motifs set at 150. The MotifDB included was Rhamnaceae Plant and the rest were excluded. For the ESI positive dataset, the the parameters were changed as follows: LDA free motifs set at 200. The MotifDB included were GNPS, Massbank and Euphorbia.

Finally, MolNetEnhancer (Ernst et al., 2019) which incorporated outputs from both FBMN and in silico tools such as MS2LDA, DEREPLICATOR, NAP, and the automated chemical classification through ClassyFire, was applied, providing thus a holistic chemical overview of measured metabolomics spectral data, with enhanced structural details for each fragmentation spectrum. The GNPS job links are provided in the Supplementary materials.

Pathways analysis and relative quantification

The metabolomic pathway analysis (MetPA) tool in MetaboAnalyst, version 5 (Pang et al., 2021) was used to perform functional analysis, particularly pathway analysis. The latter was computed using identifiers of the annotated metabolites, KEGG IDs, as input data. For overrepresentation analysis, the enrichment method, hypergeometric test was used; and for the node importance measure, i.e., topological analysis, the relative betweenness-centrality was employed. For visualization was done using scatterplots and the pathway library used was Arabidopsis thaliana (thale cress) (KEGG). Using integrated peak areas, (relative) quantitative analysis was done, and colour-coded heatmap generated in MetaboAnalyst. Pareto-scaling and log-transformation were applied as data pre-treatment methods.

Cell culture

The MIA PaCa-2 pancreatic cancer and HEK 293 human embryonic kidney were cultured in Dulbecco’s Modified Eagle Medium (DMEM) (Sigma Aldrich, USA) supplemented with 10% fetal bovine serum (FBS) (Biogen, UK) and 1% Pen-Strep (penicillin-streptomycin) (Biowhitakker, Germany) at 37 ºC in 95% humidity and 5% CO₂. About 75% of confluent cells were harvested and rinsed with 6 ml Dulbecco’s Phosphate-buffered Saline (DBPS) (Sigma Aldrich, USA). The cells were then trypsinized with 2 ml of 1 X trypsin and were incubated for 45 s (HEK 293) and 90 s (PaCa-2) at 37 ºC in 95% humidity and 5% CO₂. Four millilitres of fresh DMEM were added to both cell lines respectively to stop trypsinization. Into a new 50 ml centrifuge tube, 6 ml of cell solution was transferred into the centrifuge tube, and the cell solution was centrifuged for 4 min at 2200 rpm (779 × g). After centrifugation, the pellet was resuspended in 3 ml of DMEM for cell quantification and further subculturing.

After cell culturing, cell quantification was done using trypan blue and TC20 automated cell counter (Bio-Rad). In a 1.5 ml microcentrifuge tube, 10 µl of cells in DMEM was gently mixed with 10 µl of trypan blue and was followed by aliquoting 10 µl of the mixture into the cell counter slide compatible with the TC20. Cell viability and cell concentration were read and cells that exhibited 95% viability were used for cytotoxicity and caspase activity assays.

Extract preparation

The 80% methanol flower extracts of AMNH and RDC were prepared using the protocol reported in Sect. 2.2. The 80% methanol plant extracts were then dried into solid using the rotary evaporator and the two extracts (AMNHF and RDCF) were stored at room temperature until further analysis. The dried extracts were made up to a stock concentration of 100 mg/ml w/v in 100% dimethyl sulphoxide (DMSO) (PanReac AppliChem, Germany) respectively. The extract stock solutions were stored in the dark at 4 ºC awaiting further analysis.

Alamar blue cellular viability assay

The HEK 293 and MIA PaCa-2 were plated in 96 well plates (at 5000 cells per well). The cells were treated with AMNHF and RDCF and were incubated for 24 h with an untreated control, the negative control set as 0.1% DMSO and the positive control set as 100 µM etoposide (Sigma Aldrich, Germany), and a media blank column (plated per column). The AMNHF and RDCF extracts were used to treat the cells in a concentration-based series starting from 100 µg/ml and serially diluting it 2X to 3.125 µg/ml. A volume of 10 µl of the Alamar blue (Thermofisher, USA) reagent was added to each well following the 24 h treatment and the plate was incubated for 2 h at 37 ºC. The fluorescence was measured in each well using a plate reader and the Gen 5 program 530/25 nm excitation and 590/35 nm emission. The analysis was repeated three times.

10 caspase Glo® 3/7 activity assay

MIA PaCa-2 and HEK 293 cells were plated and treated with AMNHF and RDCF in 33 mm petri dishes in a similar manner as described in Sect. 2.9 with the IC50s and controls. The MIA PaCa-2 cell line was subjected to a caspase detection assay to determine if caspase-3 and -7 were activated due to the compound treatment. The cells were plated for 24 h in 33 mm petri dishes at 1 × 10⁵ cells/ml. The cells were treated 24 h after plating with 2 ml treatments consisting of an untreated control, the positive and negative controls, and the AMNHF and RDCF treatments with the IC50s of each extract. The assay was repeated three times and was conducted using the Caspase Glo® 3/7 assay kit (Promega, USA) based on the manufacturer’s recommendations.

Global chemical profiling – key constituents of Cannabis metabolome

The chemical profiling of plant extracts such as Cannabis, achieved through GNPS-based molecular networking (MN) tools, can assist in bettering and laying the foundation for the differentiation of Cannabis cultivars into chemovars based on their chemical class content. In doing so, metabolite biomarkers can be identified for each chemovar, and this would aid in understanding the non-generic medicinal or bio-active potential elicited by each Cannabis chemovar or plant-tissue thus narrowing down the medicinal applications of the chemovar of interest. Thus, in this study, the metabolic inventory of the two cultivars was investigated. The obtained spectral data were submitted to mass spectral molecular networking (MN) through the GNPS ecosystem, a scalable workflow that digitalizes the diversity and distribution of metabolites in plants (Kang et al., 2019).

Feature-based molecular networking (FBMN) enabled the semi-automated putative annotation of metabolites with confidence level 2 (or level 3) annotations as defined in the proposed minimum reporting standards of the metabolomics standards initiatives (MSI) (Sumner et al., 2007) (Fig. 1). FBMN represents a computational strategy that facilitates the visualization and compound annotation of complex, high-resolution untargeted LC-MS/MS metabolite data from natural extracts (Hammerle et al., 2021). This type of MN is designed to distinguish structural isomers by incorporating features such as chromatographic retention times which enhances metabolite annotation and thus the dereplication of metabolites whilst also retaining semi-quantitative information to perform statistical analyses (Quinn et al., 2017; Nothias et al., 2018, 2020). The computed FBMN of spectral data from leaf samples (from both cultivars) comprised 2914 nodes (metabolite features with an MS//MS spectrum) where 997 nodes clustered into 128 molecular families (MF) based on spectral similarity while 1917 nodes formed singletons (stand-alone nodes that are not incorporated into a MF) (Fig. 1A). Searching spectral libraries in GNPS, 166 nodes (from the 2914 nodes) were putatively annotated, with varying relative distributions across the cultivars, which point to diverse chemistries in the leaves of these Cannabis cultivars. For the samples from the flower tissues, the computed FBMN was of 2829 nodes of which 119 nodes were putatively annotated and 998 nodes clustered to form 172 molecular families while 1831 formed singletons (Fig. 1B). The singletons (selfloop nodes) positioned at the bottom of both networks (Fig. 1A-B) represented spectra that were not clustered into molecular families, i.e., low spectral similarity scores.

The molecular families (MFs) in the networks provided a global visualization and insights on the metabolome of the two cultivars, revealing diverse metabolite classes including phospholipids, cannabinoids and flavonoids, which are differentially distributed across the cultivars and across the tissues (Fig. 1 and Figure S2). From the 128 MFs found in the leaf data (where two or more connected nodes are considered as an MF), 94 MFs were comprised of unknown metabolites. However, 34 MFs had at least one metabolite automatically matched to the GNPS libraries, a similar trend that was also observed in the MFs found in the flower MN. Considering that MN is rationally constructed based on spectral similarities (Aron et al., 2020; Vincenti et al., 2020; Yu et al., 2022; Zhang et al., 2023), automated annotation of one node in a cluster can aid in decoding and annotating other structurally similar metabolites or features in the same molecular family. As such, MN improves on metabolite annotation, decoding ‘dark matter’ in spectral data, which subsequently provides an improved coverage on the annotated metabolome. In this study, this is illustrated by Fig. 2A where one spectral node in the MF was putatively annotated through GNPS spectral library matching as CBD. Based on the concept of molecular networking, theoretically this cluster posed as a cannabinoid cluster formed by metabolites with similar structures and fragmentation patterns. This then propagated the manual annotation of unidentified metabolites such as cannabidiolic and cannabidivarinic acids in the cluster (Fig. 2A). The unannotated metabolites in this cluster also suggest that there may be other unknown CBD related compounds or analogues present.

Furthermore, as may be expected from plant samples, flavonoids were also one of the major metabolite classes found in the Cannabis tissues. In this flavonoid cluster, there are spectral nodes (such as m/z 489.10, 505.28, 479.15, 639.23 and 669.31) that could not be annotated through GNPS spectral library matching, yet they are postulated to be structurally related to glycosylated flavonoids (Table S1) due to their grouping with flavonoids such as isoquercitrin (m/z 433.07) and isorhamnetin-3-glucoside (m/z 477.09) (Fig. 2B) (Table S1).

As above mentioned, it is worth noting that applying these computational tools improves annotation as it provides insights into possible molecular families in the measured metabolome. However, there are some limitations, related to similarity scoring algorithms and the complexity of untargeted metabolomics spectral data. To improve on this, the “chemical baiting” concept was explored and applied. The baiting approach is a concept that is predominantly used in protein biochemistry where chemical baits or probes are used to bind or concentrate a diverse range of biomarkers (proteins & peptides, metabolites, lipids & fatty acids, nucleic acids, and post translationally modified peptides) that are either occurring in low concentrations or masked by other dominant resident proteins (Luchini et al., 2010). Nonetheless, within metabolomics approaches based on molecular networking (MN), a known exogenous compound (a bait) can be employed to screen for related endogenous metabolites residing within the same molecular network. Thus, in this study, the vitexin standard was used as a probe to determine the effect of chemical baiting on FBMN (Fig. 3).

The overall observed effect of chemical baiting on FBMN (Fig. 3) showed that using a known refererence standard, in this case the vitexin standard, can aid in probing or fishing out known, “unknown knowns” and “unknown unknown” that are structurally similar and thus have similar fragmentation patterns depending on one or more standards used. As such, the use of molecular baits would speed up the annotation process, increasing the accuracy in forming molecular families or clusters (based on structural similarities) in the network. This was shown by the clustering of vitexin and vitexin-2-O-rhamnoside into one molecular family when the vitexin standard was added as a chemical bait.

We also note that metabolite features annotated with cannflavins A and B, flavonoids unique to Cannabis cultivars and structurally related to vitexin, were also initially detected in the same MF as vitexin-2-O-rhamnoside before the addition of the vitexin standard (a molecular bait). However, once the vitexin standard was added, the cannflavins formed their own MF (Fig. 4). It was previously demonstrated that cannflavins (canflavin A, B and C) are geranylated flavones postulated to be unique to Cannabis cultivars (Rea et al., 2019; Bautista et al., 2021; Tomko et al., 2022).

In our study, cannflavins A and B were detected in both the leaves and flowers of AMNH and RDC and reported in Table S1. Cannflavin A (m/z 435.1825) was distributed evenly across the plant tissues of both cultivars, and cannflavin B (m/z 367.1193) was recognisably present in the leaves and flowers of RDC but significantly abundant in the flowers of AMNH (Fig. 4). Using the GNPS analogues library search, the unknown metabolites in the cannflavin cluster were matched to flavone analogs such as afrormosin and 3’,3’,4’-trimethoxyflavone which are structurally similar to the detected cannflavins.

Thus, the computed FBMN (Fig. 3) provided spectral clustering and putative annotations of the measured metabolomes from the leaves of the cultivars, revealing actionable insights and key characteristic features of the chemical space of Cannabis such as 289 spectral library matches of flavonoids, cannabinoids, phospholipids and more (Table S1), 1265 connected nodes and 192 molecular families. Furthermore, chemical baiting assisted, for instance, in more confident annotation of vitexin-related structures. Moreover, the generated network showed differential distribution of metabolites across the cultivars. For instance, clusters of lipid-like molecules were more abundant in the plant tissues of RDC and flavonoids such as isoquercitrin were more abundant in the tissues of AMNH, while a cluster of hydroxy acids such as citric acid occured in similar levels across both the cultivars (Figure S2).

To further explore and characterize the metabolome of the Cannabis cultivars, in silico annotation tools such as substructure recognition topic modelling through the MS2 Latent Dirichlet Allocation (MS2LDA), and the network annotation propagation (NAP) were performed for all spectral datasets from both cultivars, AMNH and RDC. To illustrate the contribution of these tools to the Cannabis metabolite annotation and identification at a scaffold diversity level, the outputs obtained from both AMNH and RDC dataset are reported herein (Fig. 5 and Table S1). The machine learning (ML)-based tool, MS2LDA, allows an unsupervised decomposition of fragment spectra, discovering mass fragmental patterns of co-occurring mass fragments as well as neutral losses (termed ‘Mass2Motifs’, or shortly ‘m2m’) from different MS/MS spectra, which enables extracting information on substructural diversity within each class of metabolites (van der Hooft et al., 2016). During MS2LDA analyses, sets of previously annotated Mass2Motifs can be added from MotifDB (Rogers et al., 2019), whilst unknown chemistry is captured in Mass2Motifs that are numbered to enable highlighting and further studying them during subsequent analysis (van der Hooft et al., 2020). As such, the generated substructural information highlights the (bio)chemical relationships of the connected compounds. Moreover, these (shared) substructures or scaffolds can point to functional groups and/or core structures of the compounds, potentially revealing common biosynthetic routes in a molecular family (van der Hooft et al., 2016; Beniddir et al., 2021; Ramabulana et al., 2021).

In total, one hundred and eighty-one (181) Mass2Motifs were discovered, and most were common in both the Cannabis leaves and flowers datasets (ESI negative) of the two cultivars, respectively. The Mass2Motifs, annotated and unannotated were visualized through the MS2LDA website (ms2lda.org) (Wandy et al., 2018). As an illustration, some of the characterized Mass2Motifs found in the Cannabis leaves are highlighted and shown in Fig. 5. These discovered substructural Mass2Motifs are significant in the sense that they can aid in the structural elucidation of “unknown knowns’’ (metabolites or spectral nodes whose reference spectra are not found in the GNPS spectral libraries and therefore represented as unknowns) such as m/z 609.1476 identified as luteolin-glucoside-arabinoside. When evaluating the chemical makeup of the glucosylated flavonoids, Fig. 5A illustrates that m/z 609.1476 shares a common substructure (unknown motif_61) with the manually annotated metabolites luteolin-glucoside and luteolin-rutinoside. This resulted in the structural elucidation of luteolin-glucoside-arabinoside (m/z 609.1476) and led to the shared unknown motif_61 being annotated as a substructure that is related to luteolin (m/z 285). Additionally, this also suggests that motif_107, motif_48, and motif_45 are all related to some sugar loss, respectively.

Moreover, it can be postulated that motif_ 76 and motif_42 are both related to the apigenin scaffold as they both form the core structure of the putative annotations of vitexin (m/z 431.0996) and vicenin II (m/z 593.1529) where motif_76 is found in compounds related to vitexin while motif_42 is found in compounds related to vicenin II. This implies that the mass fragmental pattern captured in motif_76 is related to apigenin with mass features specific for vitexin, whereas the pattern captured in motif_42 is related to apigenin with mass features specific for vicenin. Furthermore, MS2LDA also highlighted some of the chemical moieties that make up cannabinoids as shown in Fig. 5B. It is worth noting that the fragmentation patterns (in mass spectrometry) of C-linked flavonoids, such as vitexin display key structural information. The latter can aid in differentiating flavonoids by examining first-order ion spectra or low-energy CID spectra. However, a detailed account on this (C-linked and O-linked glucosides) is beyond the scope of this study, and the reader is referred to the literature (Ma et al., 2000; Cuyckens & Claeys, 2004; Cao et al., 2014; Geng et al., 2016). It suffices here to mention that O-glycosides, linked through an oxygen atom, frequently undergo glycosidic cleavages, producing a variety of sugar moiety fragments and prominent Y ions formed by rearrangement reactions at the interglycosidic bonds. On the other hand, C-glycosides (such as vitexin), with a more stable carbon-carbon linkage, mostly molecular ions are observed together with cross-ring cleavages of the saccharidic residue and the loss of molecules of water, i.e., mostly X ions, in the nomenclature model by Domon and Costello (Cuyckens & Claeys, 2004).

Furthermore, in addition to flavonoids, compounds such as cannabidiol (CBD) were annotated through GNPS spectral library matching and with the use of FBMN, it was shown to be structurally similar to cannabidiolic acid (CBDA – m/z 357.2074) and cannabidivarinic acid (CBDVA – m/z 329.1765) (Fig. 2). The MS2LDA discovered substructures of these clustered metabolite features clearly describe the chemical relationship between the compounds through rhamn_motif_172.m2m which describes the carboxylic acid that is attached to CBDA and CBDVA (Fig. 5B). Furthermore, this led to the shared unknown motif_86 being annotated as a substructure that is related to cannabidiol (m/z 313.2175). Combining Mass2Motif annotations and library matches, we can observe how the known compounds and unknown cannabidiol analogues are associated with the cannabidiol substructure. This makes sense as the cannabidiol derivatives would have a cannabidiol core structure. Altogether, the here highlighted Mass2Motifs are indicative of flavonoids and cannabinoids, respectively, and their decorations (i.e., glycosylations). Consequently, this confirms and gives confidence to the FBMN annotations such as those previously made on the highlighted molecular family (Fig. 2) and listed in Table S1.

In addition to MS2LDA, Network Annotation Propagation (NAP) was also applied to facilitate metabolite feature annotation and to serve as an input for chemical compound classicifation performed by MolNetEnhancer. NAP performs in silico fragmentation-based metabolite predictions or annotations of candidate structures through searching several databases including GNPS, ChEBI (Chemical Entities of Biological Interest), SUPER NATURAL (II) and DNP. Therefore, using these two scoring methods, NAP is a beneficial tool when experimental spectra are matched to a few spectral libraries matches as it also allows the propagation of annotations in the absence and presence of library spectral matches (da Silva et al., 2018; Nephali et al., 2022; Kim et al., 2023). The outputs of NAP obtained in this study are reported in the supplementary material. In brief, when looking at the annotation of individual features, we could verify that NAP was able to correctly predict 18 (out of 997 spectral features) metabolites for the leaves data and 25 (out of 998) for the flower data. In both leave and flower exctracts, the NAP predicted metabolites include, cannabidiol, vitexin, orientin, and luteolin. These metabolite annotations were verified through manual inspection of fragmentation patterns and literature. Whilst other NAP annotations could not be verified as being exact matches, many are likely representing structurally related molecules present in compound databases, and all structural annotations were used by MolNetEnhancer to infer chemical compound classes of the molecular families in the generated networks.

To get an improved chemical compound class annotation and a comprehensive chemical overview of the measured spectral data, the outputs from FBMN, MS2LDA, NAP and DEREPLICATOR + were combined in an enhanced molecular network workflow, MolNetEnhancer, which integrates also the automated chemical classification of molecular families with ClassyFire (Ernst et al., 2019). MolNetEnhancer enabled the chemical annotation, visualization, and discovery of the subtle substructural diversity within molecular families (Fig. 6).

The generated networks indicated the varying distribution of major superclasses such as phenylpropanoids & polyketides, lipids and lipid-like molecules, organic oxygen compounds, benzonoids and organoheterocyclic compounds across the leaves (Fig. 6A) and flowers (Fig. 6B). A major difference observed between the plant-tissues of the cultivars (Fig. 6) was the abundance of phenylpropanoids and polyketides superclass, which was the second most dominating superclass in the leaves of the cultivars (Fig. 6A). Phenylpropanoids and polyketides include flavonoids, lignins, phenolic acids, stilbenes and coumarins. A considerable number of flavonoids, a chemical class whose biosynthetic pathway was identified as one of the significant metabolic pathways in the two cultivars (Table S2 and Figure S3), were identified in the leaves. These metabolites, in addition to their plant functions, exhibit biological activities such anti-inflammatory and antioxidant properties.

Other polyphenolic compounds such as feruloyl quinic acid and caffeoylquinic acid reported in Table S1 were also detected in the leaves. These phenolic acids which form part of chlorogenic acids are known to possess anticarcinogenic properties due to their antioxidant activities (Makhafola et al., 2016; Bouyahya et al., 2022). When evaluating the superclasses present in the flowers of the two cutivars, the Phenylpropanoids and polyketide superclass (comprised of polyphenolic compounds) was less dominant. Contradictory to this observation, flowers are regarded as important reproductive organs thus high levels of polyphenols are expected (Piccolella et al., 2020). However, the clear lack of polyphenols in the flowers of the studied cultivars (Fig. 6B) poses questions on the presence of intact biosynthetic routes of these compounds or the effects that the breeding process could have had on the chemistry of these cultivars. Moreover, literature regarding polyphenols in the flowers of Cannabis cultivars remains scarce and limited (Izzo et al., 2020). It has been shown that Cannabis cultivars lose genetic variation due to domestication and excessive breeding for selective traits (Clarke & Merlin, 2016). Therefore, such losses could contribute to the observed differences in the distribution of the phenylpropanoids and polyketides (polyphenols) between the leaves and flowers of the studied cultivars.

Additionally, Cannabis genomics data lacks plant-tissue-specific data (Hussain et al., 2021) but the metabolite data generated in this study could lay the foundation for further exploration of the regulatory genetic networks of chemical classes in plant tissues of Cannabis cultivars. For example, when zooming in on these identified superclasses and detailing the relative abundance of the identified metabolites, Figure S2 highlights the cultivar-specific and plant-tissue-specific metabolite differences where lipid-like molecules such as Glc-Glc-octadecatrienoyl-sn-glycerol (isomer 2) are dominant in the tissues of RDC while flavonoids (phenylpropanoids and polyketide superclass) such as isoquercitrin are more dominant in the tissues of AMNH. These metabolic differences could be related to changes in regulatory networks underlying the specialized metabolome of these cultivars.

The chemical profiling achieved through the GNPS-molecular networking tools revealed a complex and diverse phytochemical space of the cultivars based on the leaf and flower metabolite profiles. The leaves and flowers were characterized by the presence of chemical classes such as flavonoids, cannabinoids and phospholipids which have various medicinal properties as discussed above. Possible anti-cancer properties are one of the medicinal properties that were postulated from the revealed metabolomes of the two cultivars. To further explore this postulation, in silico molecular docking studies were performed to computationally predict the anti-cancer or anti-proliferative properties of the revealed cultivar chemistries supported by biological assay studies.

University of Venda, Biochemistry Department are gratefully thanked for access to the SHIMADZU LCMS-9030. The South African national research fund (NRF) is highly thanked for bursary support to A.M.