To get a general idea of the metabolite profile of leaves and cotyledons from B. napus, C. sativa and T. arvense (field pennycress), we carried out untargeted metabolomics through LC-MS analysis on a Chemical Isotope Labeling (CIL) Metabolomics Platform (Zhao S. et al., 2019). In total, of the thousands of metabolites that were detected 718 could be classified with confidence. The metabolites were classified into the following categories: short and medium-chain fatty acids, polyamine, neurotransmitter, phytohormone, phenol and quinone, flavonoid, vitamins and derivatives, phenylpropanoid, polyphenol, lipid, amino acids and derivatives, dipeptides and tripeptides, alkaloid, terpene, and others ( Table 1 ). The most abundant two categories are amino acids and derivatives, dipeptides and tripeptides, which contributed 131 (18.2%) and 122 (16.9%) compounds respectively.

To have an overview of the metabolites’ distribution across all examined samples, we performed a hierarchical clustering analysis, the clustering results indicated that the 718 metabolites could be classified into 10 subclasses according to their relative abundance ( Figure 1 ). The abundance level of metabolites in each subclass varied to a great extent across the experimental samples. Notably, subclasses 2, 3, and 7 exhibited distinctive patterns, which effectively distinguished B. napus, C. sativa, and pennycress from one another. Metabolites in subclass 2 displayed significantly higher levels of accumulation in both Camelina’s leaf and cotyledon tissues; metabolites in subclass 3 demonstrated higher levels in both B. napus lines compared to C. sativa and pennycress. In contrast, metabolites in subclass 7 showed a much greater accumulation in pennycress ( Figure 2 ). Analyses of each of the three differential subclasses showed no significant enrichment for metabolites in any one biochemical pathway. The specific details regarding the 10 subclasses as well as their included metabolites can be found in Supplementary Data Sheet 1 .

To determine the variability between different groups (species/tissue), as well as the metabolic differences among the three biological replicates for each group, a principal component analysis (PCA) were performed with all detected 718 compounds. Figure 2A showed that the first principal component (PC1) accounted for 30.1% of the variation, and largely separated pennycress from the other two species. Meanwhile the second principal component (PC2) accounted for 16.8% of the variation and separated B. napus from C. sativa ( Figure 2A ). The PCA plot shows that for each group, three biological replicates are highly clustered together, suggesting there is a high cohesion within each group. Meanwhile, the samples from the three oilseed crops separated into three distinct areas in the plot, indicating each species possessed a distinct metabolic profile overall. Further, two B. napus lines, DH12075 and N99, were closely grouped together, implying there was limited differences between these two genotypes with regards to the metabolic profile ( Figure 2A ). To further confirm the reproducibility, a sample-sample correlation analysis was performed, where the resultant heatmap showed very high correlation between the three biological replicates for each group. Although, for each species, the metabolic profiles of their cotyledons and leaves were very similar, it was clear from the B. napus data that the individual tissues from the two genotypes were more similar to each other than the differences between the genotypes ( Figure 2B ). Yet, as the main goal was to identify differentially accumulated metabolites at the species-level, results from leaves and cotyledons were combined for each species during subsequent analysis. Additionally, since the two B. napus genotypes (i.e. DH12075 and N99) possessed very similar metabolites overall, DH12075 was used to represent the oilseed B. napus during the following analysis.

To check the overall metabolic differences between these three oilseed crops, we performed pairwise comparisons among them. Each compound was compared between two species by a Student’s t-test, the p-value of which was further corrected through the Benjamini-Hochberg procedure. Significantly different metabolites were selected based on their adjusted p-values with 0.05 as a cut-off value. During the three pairwise comparisons, a total of 267 compounds showed differential abundance when comparing B. napus vs. C. sativa, 383 for B. napus vs. pennycress, and 376 for C. sativa vs. pennycress ( Figure 3A ; Table 2 ). The comparisons revealed very similar results between B. napus vs. pennycress and C. sativa vs. pennycress, with the number of their differential metabolites reaching 383 and 376 respectively, among which they shared 268 common metabolites ( Table 2 ; Figure 3B ). Meanwhile, the smallest difference was detected in comparisons of B. napus vs. C. sativa, where only 267 compounds showed differential abundance (170 had higher abundance in B. napus while 97 were higher in C. sativa) ( Figures 3A, B ; Table 2 ).

After the explorative analysis above, our main question concerned those metabolites that showed a discriminant pattern among these three oilseed crops. The partial least squares-discriminant analysis (PLS-DA) is an effective method for separating differentially abundant metabolites (DAMs) between different groups in metabolomics data because of its ability to handle highly collinear and noisy data, and maximize the differences between groups (Gromski et al., 2015; Qi et al., 2022). To find the most discriminant metabolites between B. napus, C. sativa and pennycress, we carried out a PLS-DA with all 718 identified metabolites. According to the variable importance in projection (VIP) scores from PLS-DA, DAMs were initially selected if their VIP scores were higher than 1. Further, all initial DAMs were filtered by an ANOVA to only keep the ones with p-value smaller than 0.05. Thus, DAMs were defined as metabolites with VIP score > 1 and ANOVA p-value < 0.05, which resulted in a total of 149 final DAMs ( Supplementary Data Sheet 2 ). Based on the relative metabolite abundance, these 149 DAMs were divided into three distinct sub-groups ( Figure 3C ). Among the three sub-groups, there were distinct metabolite distributions: sub-group 1 comprised 47 metabolites, exhibiting the highest abundance in C. sativa; sub-group 2 encompassed 52 metabolites, demonstrating highest concentration in pennycress; 50 metabolites in the 3rd sub-group, displayed the highest concentration in B. napus instead ( Figure 3C ).

To check the top DAMs among these three crops, we selected the 20 metabolites with the highest VIP scores, and plotted these metabolites along with their relative abundance across all three oilseed crops for visualization ( Figure 4A ). Agmatine, a precursor for polyamine biosynthesis, has the highest VIP score, and its abundance was significantly higher in C. sativa compared with the other two oilseed crops ( Figure 4A ). To further gain knowledge of what biological pathways these 149 compounds were involved in, we performed a functional enrichment analysis according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database through the pathway enrichment analysis of MetaboAnalyst6 (https://www.metaboanalyst.ca/). As shown in Figure 4B , several amino acid metabolism pathways, including “Tyrosine metabolism”, “Arginine and proline metabolism”, and “Cysteine and methionine metabolism” have very small p-values (p-value < 0.01), which is in line with the large number of overall “Amino Acids & Derivatives” identified in the untargeted metabolomics results ( Table 1 ). Meanwhile, “Isoquinoline alkaloid biosynthesis” also had a small p-value (p-value < 0.05). More noticeably, “Phenylpropanoid biosynthesis” pathway has the smallest p-value and the ensuing “Flavonoid biosynthesis” and “Flavone and flavonol biosynthesis” also have small p-values (p-value = 0.056 and 0.037 respectively).

To investigate the candidate genetic regulators related to the abovementioned DAMs, we performed a concomitant RNA sequencing along with the metabolomics. Since the three species vary in ploidy level, in order to compare the gene expression between these three oilseed crops, all their genes were first projected onto their corresponding orthologs in A. thaliana, and the sum of the normalized expression of all orthologous genes in B. napus and C. sativa was calculated for each A. thaliana gene (see details in Methods). Although it is recognized that there may be differential expression among the orthologues in B. napus and C. sativa, summing the expression across the duplicate copies allowed us to identify those genes which could play a role in differentiating the species. In total, there were 14,025 unique A. thaliana genes possessing corresponding orthologous genes in all three oilseed crops, which were selected for further comparison ( Supplementary Data Sheet 3 ). To identify the most differentially expressed genes (DEGs) among these three crop species, a likelihood ratio test (LRT) in DESeq2 and a PLS-DA were employed for these 14,025 genes. Under the definition of DEG as “LRT adjusted p-value < 0.05 and PLS-DA VIP score > 2”, there were 644 DEGs in total. Further, a hierarchical clustering was performed for these 644 DEGs based on their expression profiles. As shown in Figure 5A , these DEGs’ expression patterns varied to a great degree, and the three oilseed crops could be distinguished clearly according to the expression level of these DEGs ( Figure 5A ). To functionally characterize these identified DEGs, a KEGG enrichment analysis was performed, where pathways “2-Oxocarboxylic acid metabolism”, “Glucosinolate biosynthesis” and “Tropane, piperidine and pyridine alkaloid biosynthesis” etc. were significantly enriched ( Figure 5B ).

As the 149 DAMs above were mainly attributed to the “Phenylpropanoid biosynthesis” pathway ( Figure 4B ), to further examine this pathway in more detail, we selected all DAMs that were classified as phenylpropanoids or flavonoids by the KEGG database, and compared their relative abundances across these three oilseed crops. The compounds were involved in general phenylpropanoid pathway (e.g. p-Coumaric acid), flavonoid pathway (e.g. Afzelin, Vitexin, Kaempferol-3-O-galactoside), and phenolic acid pathway (e.g. Chlorogenic acid). The relative abundances of these phenylpropanoids varied to a great extent, for example, C. sativa contains the highest levels of Chlorogenic acid and cis-3,4-Leucopelargonidin, while B. napus contains the highest levels of Kaempferol-3-O-galactoside, trans-Ferulic acid, and p-Coumaroyl quinic acid, meanwhile Vitexin and Afzelin accumulated the most in pennycress ( Figure 6 ).

To identify genes that were responsible for the abovementioned phenylpropanoids, multiple putative regulatory genes involved in phenylpropanoid biosynthesis were selected to compare their expression levels. As might be expected, there was great variation in expression levels for phenylpropanoid biosynthetic genes between these three oilseed crops ( Figure 7 ); for example, homologs of flavonol 3’-hydroxylase (F3’H/TT7) and dihydroflavonol reductase (DFR/TT3) were predominantly expressed in C. sativa, meanwhile homologs of leucoanthocyanidin dioxygenase (LDOX/ANS/TT18) were mainly expressed in B. napus, these expression differences were in line with the high accumulation of cis-3,4-Leucopelargonidin in C. sativa ( Figures 6 , 7 ). However, given the limitation of untargeted metabolomics, there were only a small amount of identified phenylpropanoids, which makes it difficult to link the variation in gene expression directly to the corresponding change in phenylpropanoids.

To explore the 644 DEGs in a more general manner, co-expressed gene modules were identified through weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008). As shown in Figure 8A , these 644 DEGs were divided into 9 separate modules based on their expression levels. According to the eigengene expression of each module, module red, green and yellow were B. napus, pennycress and C. sativa specific, respectively ( Figure 8B ). Interestingly, C. sativa-specific module yellow was highly correlated with the 47 DAMs in sub-group 1, which accumulated the highest abundance in C. sativa ( Figure 3C ; Figure 8C ). Similarly, B. napus-specific module red and pennycress-specific module green were strongly correlated with 50 DAMs in sub-group 3 and 52 DAMs in sub-group 2 ( Figure 3C ; Figure 8C ), thus suggesting that the species-specific gene expression identified among these three oilseed crops might contribute to their distinct metabolic profiles. To mine potential candidate genes associated with variation of the abovementioned phenylpropanoids, we further examined the correlations between the identified phenylpropanoids and the 644 DEGs above. Pearson’s correlation coefficients showed that 171 DEGs were significantly related to 10 phenylpropanoids (R² > 0.9), whereat 104 pairs were positively correlated, and 76 pairs were negatively correlated ( Supplementary Data Sheet 4 ), these 171 DEGs were highly enriched in GO terms “secondary metabolite biosynthesis process” and “glycosyl compound metabolic process”.

In summary, to document the metabolic profiles of B. napus, C. sativa, and field pennycress, a total of 718 metabolites including flavonoids, phenylpropanoids, amino acids and derivatives, were characterized through untargeted metabolomics analysis. Moreover, 149 DAMs were identified through PLS-DA, indicating the different metabolic profiles of these three crops were mainly attributed to phenylpropanoids. Further integration with transcriptomics revealed a great variation of the candidate genes’ expression between these three species. Our data provides a novel resource for understanding the metabolism of these three oilseed crops at the molecular level.

Although, metabolomics has been widely used to study the metabolites of plants (Allwood et al., 2021; Wang et al., 2016; Yan et al., 2022), in oilseed crops, most attention has focused on dissecting the metabolic profile of their seeds, in order to further improve the oil quality. For example, Boutet et al. (2022) revealed a large diversity of specialized metabolites in Camelina seeds through untargeted metabolomics and lipidomics (Boutet et al., 2022). Li et al. (2023) utilized targeted metabolomics to study the oil content of Brassica napus seeds, which identified marker metabolites that were correlated with oil content (Li et al., 2023). However, there is very limited knowledge about the metabolites in the leaf and/or cotyledon tissues of these oilseed crops. Metabolites play key roles in plants’ physiological and biochemical responses to environmental stimuli, and the accumulation of metabolites also has been shown to vary to a great extent in different species (Qi et al., 2021; Wang et al., 2016). Understanding the metabolite differences between species in tissues such as cotyledons and young leaves, which are often the first to recognize new threats, could provide insights into why some species are more resilient than others. In this instance, two novel oilseed crops C. sativa and T. arvense that have often been associated with higher tolerance to a number of biotic and abiotic stresses compared to their more widely grown relative B. napus were used to test this hypothesis.

In this study, we detected more than 3,465 metabolic peak signals through LC-MS-based untargeted metabolomics ( Supplementary Data Sheet 5 ). Underlying the current limitations of these technologies, only 718 metabolites could be identified given that the majority of metabolites in the plant kingdom are still unknown or unannotated (Wang et al., 2019; Yonekura-Sakakibara and Saito, 2009). The analysis of these 718 identified metabolites revealed great variation between these three oilseed crops, indicating a species-specific accumulation pattern of metabolites ( Figure 1 ), which was further supported by the PCA plot that showed complete separation of these three species ( Figure 2A ). The 149 DAMs identified by PLS-DA were highly enriched in “Phenylpropanoid biosynthesis”, and amino acid metabolism pathways ( Figure 4B ), suggesting that both phenylpropanoids and amino acids play important roles in distinguishing these three oilseed crops.

Utilizing untargeted metabolomic analysis to these three oilseed crops provides us a comprehensive description of their metabolic variations. To mine the putative genes related to their distinct metabolic patterns, a cross-species transcriptomic analysis was employed between these three species to identify DEGs. As shown in Figure 5A , these 644 DEGs display distinctive expression patterns across three species. Interestingly, these DEGs are significantly enriched in KEGG pathways, “2-Oxocarboxylic acid metabolism” and “Glucosinolate biosynthesis”. Glucosinolates (GLSs), the unique class of secondary metabolites that are prevalent in the Brassicaceae family members. The type and concentration of GLSs vary to a great extent between different species of the Brassicaceae family and between different tissues within the same species (Bhandari et al., 2015; Essoh et al., 2020). Moreover, a number of environmental factors, including seasons, abiotic/biotic stresses, largely affect the accumulation of GLSs (Essoh et al., 2020; Mitreiter and Gigolashvili, 2021).

For the three oilseed crops in this study, all their seeds contain significant amount of GLSs (Bhandari et al., 2015; Chopra et al., 2020; Czerniawski et al., 2021). However, only B. napus and pennycress are reported to contain large amounts of GSLs in their leaves, with no detectable GLSs in the leaves of C. sativa (Bhandari et al., 2015; Chopra et al., 2020; Czerniawski et al., 2021), which coincided with the fact that the majority of GLSs biosynthetic genes were “off” in C. sativa leaves but highly expressed in both B. napus and field pennycress based on the expression level ( Supplementary Data Sheet 3 ). The expression level of these GLSs related genes in C. sativa leaves is verified by the publicly available transcriptome atlas of C. sativa (Kagale et al., 2016), which further suggests a tissue-specific expression pattern of these GLSs related genes in roots and seeds of C. sativa. Given the essential roles of GLSs in Brassicaceae plants’ response to biotic and abiotic stresses (Del Carmen Martinez-Ballesta et al., 2013; Variyar et al., 2014), it will be very interesting to investigate what’s the evolutionary benefits of no GLSs in C. sativa leaves?

Flavonoids are a subgroup of phenylpropanoids, which are an important class of plant secondary metabolites synthesized via the shikimate pathway from phenylalanine (Dixon et al., 2002; Fraser and Chapple, 2011; Vogt, 2010). In addition to flavonoids, phenylpropanoids contain other different subgroups, including lignins, phenolic acids, stilbenes, and coumarins (Deng and Lu, 2017; Vogt, 2010). Both phenylpropanoids and flavonoids have been shown to play important roles in plant defense against abiotic and/or biotic stresses (Ramaroson et al., 2022; Sharma et al., 2019). Our current results indicated that there was great variation in the phenylpropanoids, including flavonoids and phenolic acids, between these three oilseed crops, e.g. C. sativa contains high amount of chlorogenic acid ( Figure 4B , 6 ). However, given the limitation of using untargeted metabolomics, where a significant portion of the detected peaks are not characterized, only a small number of phenylpropanoids were identified, which hinders our understanding of where the different metabolic pathways diverge in the three species. A previous study by Onyilagha et al. (2003) using a targeted approach, presented a more detailed view of leaf flavonoid distribution among different crops from the Brassicaceae family through Thin-layer Chromatography (TLC) and high-performance liquid chromatography (HPLC), whereat quercetin was the only type of flavonoid detected in C. sativa; while B. napus accumulated at least two types of flavonols, with kaempferol as the major type and quercetin in a small concentration (Onyilagha et al., 2003). A further study confirmed that B. napus leaves mainly contained kaempferols, while C. sativa accumulated large amount of quercetins instead (Onyilagha et al., 2012). Interestingly, high kaempferols in B. napus and high quercetins in C. sativa is in line with our transcriptomic data, where the homologs of TT7, which encodes a flavonol 3’-hydroxylase to convert kaempferols to quercetins, were predominately expressed in C. sativa ( Figure 7 ). Genetic editing and kinetic analysis of TT7 homologs in C. sativa could be useful for further verification of their functions.

Compared with other Brassicaceae oilseed species, C. sativa has been shown to tolerate a broad range of abiotic stresses, including but not limited to drought, freezing, lodging, and salinity (Berti et al., 2016; Čanak et al., 2020; Gao et al., 2018; Matthees et al., 2018). However, the mechanisms underlying these tolerances of C. sativa are largely unknown. Our results here indicated that there were particularly high content of polyamines in C. sativa compared to other two crops. Noticeably there are two polyamines in the top 20 DAMs, agmatine and spermidine, both of which accumulate in C. sativa ( Figure 4A ). Interestingly, both agmatine and spermidine are involved in the “Arginine and proline metabolism” pathway ( Figure 4B ), where agmatine is an intermediate metabolite formed by arginine decarboxylase (ADC) from arginine, further agmatine is used to generate spermidine through sequential reactions (Chen et al., 2018).

More importantly, both agmatine and spermidine have been shown to play important roles during plant defense against a variety of abiotic stresses, including drought, freezing and salinity etc (Gupta et al., 2013; Kasukabe et al., 2004; Seifi and Shelp, 2019; Shao et al., 2022). For example, Kasukabe et al., 2004 showed that overexpression of spermidine synthase in A. thaliana significantly increased the spermidine content in leaves, and provided enhanced tolerance to various abiotic stresses (Kasukabe et al., 2004). Therefore, it’s reasonable to speculate that the high content of polyamines (e.g. agmatine and spermidine) in C. sativa is likely to contribute to its resistance to a broad range of abiotic stresses. In the future, quantifying and comparing other major polyamines (e.g. diamine putrescine and tetraamine spermine) between these oilseed crops, as well as genetic modification of related polyamine biosynthetic genes in C. sativa; for example, targeted knock-out of agmatine and spermidine biosynthetic genes, encoding arginine decarboxylase (ADC) and Spd synthase respectively, would help us better elucidate the function of polyamines in the broad abiotic resistances of C. sativa.

Brassica napus (DH12075 and N99), Camelina sativa (DH55), and pennycress (Thlaspi arvense; line collected from field at Saskatoon, SK, Canada) were grown in a growth chamber (22°, 16 h light/18°, 8 h dark cycles) for up to 2 weeks after germination. Three biological replicates of cotyledons and first true leaves for each species were collected and ground into powder in liquid nitrogen, then stored at -80° freezer for metabolomic and transcriptomic analyses.

Untargeted metabolomic analyses was carried out by the Metabolomics Innovation Centre (University of Alberta, CA). In short, metabolites were extracted from 40 mg of tissue powder for each of three biological replicates per species and tissue, where 6 ceramic beads were placed in the sample vials, and 500 µL LC-MS grade MeOH/water (4:1 v/v) was added before homogenizing at 4.5 m/s for 15 seconds. Then, the homogenates were incubated at -20°C for 10 minutes and centrifuged at 15,000 g for 10 minutes, after which the supernatants were carefully transferred into new vials and completely dried. Sample extracts were then re-suspended in 30 µL LC-MS grade water before chemical isotope labeling. The subsequent LC-MS analyses were carried out with a Thermo Scientific Vanquish LC linked to Bruker Impact II QTOF Mass Spectrometer (Bruker, Germany) using the eclipse plus reversed-phase C18 column (150 x 2.1 mm,1.8 µm particle size; Agilent) at 40°C with a flow rate of 400 μL/min.

The raw LC-MS data were first processed using DataAnalysis 4.4 (Bruker), which then exported data to IsoMS Pro 1.2.16 (Nova Medical Testing Inc, CA) for quality check and processing with the following parameters: minimum m/z: 220, maximum m/z: 1000, saturation intensity: 20,000,000, retention time tolerance: 9s, mass tolerance: 10 ppm. Further, data files were filtered out peak pairs present in less than 80% of samples, and normalized by ratio of total useful signal. Metabolite Identification were carried out against the NovaMT Metabolite Database v3.0 (The Metabolomics Innovation Centre, CA).

The hierarchical clustering analysis, principle component analysis (PCA), and partial least squares-discriminant analysis (PLS-DA) were carried out through scikit-learn (version 1.2.2) in Python (version 3.10.9). KEGG pathway functional enrichment analysis was performed via MetaboAnalyst6 (https://www.metaboanalyst.ca/) using the A. thaliana database.

Total RNA of each sample was extracted using RNAeasy plant mini Kit (Qiagen) according to the manufacturer’s instructions, with three biological replicates per sample (24 in total). The quality of total RNA were examined by BioAnalyzer with RNA 6000 Nano Kit (Agilent) to ensure RNA integrity number value > 7. Total RNA were used to prepare cDNA libraries following Illumina Stranded mRNA Prep guide, further 150 bp paired-end sequencing was performed on the NovaSeq 6000 platform (Illumina).

The raw RNA-seq data were first filtered using Trimmomatic (version 0.32) to remove adapter and low-quality sequences. Then clean reads were aligned to the corresponding reference transcriptomes of B. napus (DH12075 v3.1; cruciferseq.ca), C. sativa (Kagale et al., 2014) and pennycress (Nunn et al., 2022) respectively using Salmon (version 1.10.0). Subsequently, tximport (version 1.30.0) was used to obtain the gene expression levels as TPM (transcripts per million), which were further normalized through log2 transformation.

Syntelog tables of the A. thaliana orthologous genes were collected for B. napus (DH12075 v3.1; cruciferseq.ca) and C. sativa (Kagale et al., 2014). Homologous pairs between pennycress and A. thaliana were obtained through Reciprocal BLAST Hits (RBH). To make the expression comparable across three different species, all genes were projected to their corresponding A. thaliana orthologs for each species. Given the polyploidy of B. napus and C. sativa, they contain multiple copies of orthologs for each A. thaliana gene. The sum of the normalized expression for all orthologous genes in B. napus and C. sativa for each A. thaliana gene was used for subsequent comparison.

To identify the most differentially expressed genes (DEGs) between these three species, a partial least squares discriminant analysis (PLS-DA) through scikit-learn (version 1.2.2) in Python (version 3.10.9) and a likelihood ratio test (LRT) in DESeq2 (version 1.42.1) in R (version 4.3.3) were performed. The KEGG enrichment analysis of DEGs was performed using clusterProfiler (version 4.11.0) in R (version 4.3.3). The 644 DEGs were used to construct a co-expression network using the WGCNA package in R with the blockwiseModules function and parameters “power= 12, maxBlockSize = 5000, networkType = “signed”, TOMType = “signed”, minModuleSize = 10, mergeCutHeight = 0.15” to construct a signed network and generate co-expressed gene modules (Langfelder and Horvath, 2008).

For integration of metabolomic and transcriptomic data, Pearson correlation coefficients were calculated for each gene-metabolite pair using the ‘cor’ package in R (version 4.3.3). The gene-metabolite network was constructed with R² > 0.9, where the nodes correspond to genes/metabolites, and the edges between the nodes represent the correlation coefficients calculated ( Supplementary Data Sheet 4 ).