Search
Close this search box.

Integrated omics of Saccharomyces cerevisiae CENPK2-1C reveals pleiotropic drug resistance and lipidomic adaptations to cannabidiol – npj Systems Biology and Applications

CBD treatment incurs a growth inflection after controls have established their stationary phases in yeast cell cultures

We first investigated the effect of cannabidiol treatment on yeast using continual growth monitoring to test if there was any evident growth inhibition or other visible phenotypic effects. To accomplish this, we monitored cell growth of S. cerevisiae CENPK2-1C cells inoculated at 0.1 ODU/mL in Erlenmeyer flasks at 30 °C and orbital rotation at 200 rpm continually throughout the logarithmic and latent phases. Treated cells were fed cannabinoids at 0.5 mM at the same time as the cellular inoculation of the flasks. Cannabinoids were introduced in an organic solvent matrix, and appropriate positive controls were implemented (depending on which experiment) to control for the residual effects of the respective organic solvent, along with negative controls with no treatment but identical culture conditions. As a general trend, we observed that the backscatter signal decreases for ~6 h (Fig. 2). However, it is important to note that these growth curves do not represent cell viability density, but signal as measured by a backscatter detector. Given that the highly lipophilic cannabinoids immediately precipitated into globules upon delivery to the aqueous cell media, we hypothesized that this decrease in backscatter represented the assimilation of the cannabidiol into and onto the cells, which would decrease the total number of free particles in the media.

Fig. 2: Evaluation of Growth in CBD-administered samples across four experimental setups.
figure 2

In each subplot, the black trace represents the CBD-exposed sample, the orange trace represents the positive control, and the blue trace denotes the untreated sample. Biomass was quantified through backscatter measurements facilitated by an Aquila continual biomass monitoring sensor placed beneath the shake flasks. The dashed gray line in each subplot indicates a critical point of growth inflection for the CBD-treated cultures. a Methanol solvent matrix and methanol positive control. Approximately at the 45-h mark, a phase transition can be observed, triggering a second growth phase for the CBD-administered cell cultures. b DMSO solvent matrix and DMSO positive control. A metabolic phase transition is again observed. c sucrose-fed cultures (methanol solvent matrix). Sucrose was selected because it is an α-linked heterodisaccharide. d minimal media (glucose, methanol solvent matrix). In all experiments, there is a persistent growth inflection found uniquely in the CBD-fed cell cultures. All experiments have been replicated a minimum of three times.

A growth inflection point persists in CBD treated yeast cultures independently of the solvent used to introduce the CBD and the nutrient environment

We evaluated the yeast growth curves using a methanol solvent matrix to introduce the cannabidiol into the liquid culture (Fig. 2a). We observed a growth inflection point at around 45 h, which did not occur in either the untreated group or the positive control (0.5% methanol). We thought that this might occur because of a post-diauxic shift to another source of carbon for growth, so we moved to manipulate the nutrient environment of the yeast through a series of experiments to determine if this phenotype persisted across multiple carbon sources.

It is well known that methanol can be metabolized by S. cerevisiae, so we decided to further evaluate this effect using DMSO as the solvent matrix (Fig. 2b). In both groups, as the cell growth curves proceeded, it was clear that both the treatment group and positive control grew to a higher density than the negative control, although it was unclear if this was due to carbon utilization of the organic solvent or if the addition of 0.5% organic solvent incurs a physical effect that promotes the growth of the yeast cultures. Given that there is no described mechanism to utilize DMSO as a carbon source in S. cerevisiae, we prefer the latter explanation, but further investigation here was outside of the scope of this study. Overall, the total density of the yeast cultures was uninteresting for the study given that the positive controls grew to similar densities. We focused instead on the growth inflection point unique to CBD-fed cultures. After establishing that this growth inflection occurs in both methanol- and DMSO- fed yeast cultures, we concluded that cannabidiol-fed cultures were able to utilize additional carbon and continue growing after the onset of stationary phase in the control groups. Therefore, we investigated this effect further by manipulating the growth media.

We replicated the study in sucrose-containing media (Fig. 2c) to determine if the effect was still visible after switching to a heterodimeric sugar and we observed that the latent stationary phase metabolic shift persisted. We then replicated the study using minimal media yeast nitrogen base (YNB) and once again, observed the same effect (Fig. 2d). These results led us to conclude that the growth inflection point caused by CBD treatment occurs independently of the nutrient environment of the yeast cells.

Acidic cannabinoids and cannabidiol analogs with shortened alkyl chains create growth perturbations when treated to yeast cell cultures

Last, we were interested to determine if other cannabinoids behaved similarly to the effects that we observed from cannabidiol treatment. We examined the reproducibility of this effect in yeast cultures which were treated with acidic cannabinoids, as well as the alkyl chain substituted analogs of cannabidiol (CBD-C5), cannabidivarin (CBD-C3), cannabidiorcol (CBD-C1), and cannabidiresoricol(CBD-C0) (Supplementary Information Fig. 1). Interestingly, each acidic cannabinoid incurred a unique effect on the yeast cultures that was distinct from the effect caused by cannabidiol. CBGA-fed cultures behaved like the positive and negative controls except that the beginning of the exponential growth phase was delayed by about 6 h. Each subsequent metabolic shift was similarly delayed, however we did not observe any additional abnormal growth inflections. THCA also behaved almost normally except for a subtle change occurring around 19 h. The pattern of backscatter signal in the CBDA treated cultures was consistent with the controls, however, interestingly, CBDA treatment led to the only example of overall lower backscatter density than the positive and negative controls. Finally, the short-alkyl chain homologs of cannabidiol suspended growth entirely until 40-52 h, at which point growth proceeded normally. Overall, these results demonstrated that each cannabinoid affects yeast culture growth uniquely.

After we determined that cannabinoids exert phenotypic influence on the growth patterns of our yeast, we hypothesized that changes in the central metabolism of the cells were the cause of this observed post-diauxic shift, so we proceeded with a broad, non-targeted integrated omics approach to define the cellular changes that occurred during this growth shift. Therefore, we grew more yeast cultures under continual biomass monitoring and harvested the cultures at the growth inflection point as it occurred. We then harvested RNA from part of the culture suspensions for transcriptomic sequencing and diverted the rest for metabolomic analysis using a UPLC/HRMS.

Transcriptomics and genomics indicate activation of a pleiotropic drug resistance response upon cannabidiol treatment at the growth inflection point

S. cerevisiae CENPK2-1C genome sequence enables strain-specific transcriptomic assembly

We first sequenced and assembled the genome of our strain to ensure that our transcriptomic data could be accurately quantified. S. cerevisiae CENPK2-1C was chosen for this study because it is a common strain used for metabolic engineering, owing to its deletions of URA3, HIS3, LEU2, and TRP1, which can be used as selection markers. However, a well sequenced and annotated genome for this strain is not available. Moreover, as big-data becomes more widely accessible to laboratories, we advocate an approach of using strain-specific curations to validate assumptions that were previously made about yeast genome stability. Given the fundamental importance of having an accurate representation of gene copy number in transcriptomics analyses, we sequenced our internal laboratory strain genome of CENPK2-1C (Euroscarf, MATa; ura3-52; trp1-289; leu2-3,112; his3-Δ1; MAL2-8c; SUC2) using Nanopore long reads for assembly and Illumina short reads for polishing. After filtering reads less than 1000 bp, our Nanopore run sequenced 477,258 total reads. The mean read length was 10,566.8. Additional information, including a histogram of the read lengths and the Nanoplot summary statistics, is available in the supplementary material. After Flye assembly29, our genome consisted of 28 contigs with an N50 of 800.5 Kbp and a genome size of 11.98 Mbp. Upon realigning the reads to the assembly using Minimap230, we observed an average coverage depth of 958. As these quality checks were satisfactory, we moved to polish the assembly with ONT’s medaka© (using long reads) and then with pilon31 (using trimmed short Illumina reads). Upon assembly and polishing, a BUSCO32 analysis indicated 99.4% completeness of the assembly (Supplementary Information Fig. 6). Last, we used Funannotate (software, open access) and the annotation transfer tool in Geneious to annotate a total of 4459 genes. To evaluate and benchmark the differences between the genomes, we created a 2D dot plot using Nucmer and DNAdiff, the results from which we have included in the supplementary material as Supplementary Fig. 62 and Supplementary Table 6. We identified a total of 23468 SNP’s, 6524 indel’s, 12 relocations, 23 translocations, and 305 insertions. While these findings indicate an overall high degree of similarity evidenced by the genome alignments, the differences are nevertheless relevant to our transcriptomics study.

Differential gene expression analysis shows eight genes are overexpressed in response to cannabidiol

We pursued differential gene expression analysis of our CBD-fed and MeOH-fed positive controls because we wanted to determine which genes were being expressed at the time point of the growth inflection. After sequencing, we evaluated our transcriptomic data sets using Nanoplot to ensure that read counts would be sufficient for quantitative analysis. The Nanoplot data for the transcriptomic cDNA libraries is provided in the supplementary information. A total of 1,448,164 reads, derived from both CBD-fed and positive control (MeOH) group libraries, were stranded using the UNAGI33 pipeline and assembled to the S. cerevisiae CENPK2-1C genome using Minimap2. In essence, UNAGI is a strand-aware transcriptome assembly pipeline designed for handling cDNA reads generated by Nanopore sequencing. Subsequently, gene expression quantification was performed using Geneious. The alignments resulted in read count tables, encompassing the analysis of 4435 genes. A significant subset of these genes was associated with biosynthetic processes, prompting the inclusion of an additional set of reads from a mid-logarithmic growth phase cDNA library data set. This additional set comprised 2,966,571 reads also assembled to the S. cerevisiae CENPK2-1C genome, providing a more precise perspective for our DESeq2 analysis.

Our comparative analysis of read sets was performed using DESeq234 (n = 6), where the CBD-fed group was the experimental focus and was compared against the MeOH-treated group sampled in stationary phase at the same time as the experimental group (positive control) and an untreated, midlogarithmic-phase group (negative control). DESeq2 was selected for its ability to effectively handle small sample sizes using a shrinkage estimator for dispersion and fold changes. We hypothesized that the mean fold change should be 0 because most genes should be expressed at the same levels. We observed a mean fold change of 0.006, which confirmed that DESeq2 was an appropriate method to compare our data sets given our experimental design. A cut-off for significance was established at a log2foldChange of >1.5 and an adjusted p-value of <0.05. Significant overexpression was observed for the following genes: CIS1, RDN37-1, PDR5, RDN25-1, RDN18-1, TAR1, IMT1, YGR035C. These results are visually represented in a volcano plot (Fig. 3).

Fig. 3: Volcano plot visualization of differential gene expression analysis.
figure 3

Points which are labeled and shown in red represent genes that were identified to have greater than 1.5 log2 fold change (LFC) and an adjusted p-value less than 0.05. Most transcripts were expressed at or near the same level. Eight transcripts exceeded the significance threshold; CIS1, RDN37-1, PDR5, RDN25-1, RDN18-1, TAR1, IMT1, YGR035C.

Pleiotropic drug resistance is expressed in CBD-treated yeast

This set includes three notable genes which are a part of the pleiotropic drug resistance35 response in yeast; PDR536, CIS137, and its paralog YGR035C. Little is known about CIS1, other than it is under regulatory control of the PDR-involved transcription factor Yrr1p and it is localized to mitochondria. It shares 21.6% sequence identity and little structural similarity with its paralog YGR035C38 (Supplementary Information Fig. 8). YGR035C is also under the regulatory control of PDR-linked transcription factors Yrm1p and Yrr1p. Interestingly, WoLFPSORT predicts subcellular localization to the nucleus (Supplementary Information Table 4). While the two proteins have little in common, there are 26 conserved residues across the peptides and a 15 amino-acid region which contains 12 conserved residues, PXKITRYDLXKAAXE. The alignment and the predicted structures with the residues highlighted in red are available in the supplementary information.

In contrast to CIS1 and YGR035C, much is known about the promiscuous ABC steroid and cation transporter PDR539. PDR5 is localized to the plasma membrane and is under the positive regulation of the transcription factors Pdr1p and Pdr3p40, however many transcription factors have been shown to induce a PDR5 response both positively and negatively. PDR5 is known to facilitate removal of drugs, steroids, antifungals and phospholipids from yeast cells and has been implicated in many studies as a critical part of the pleiotropic drug response41,42. PDR5 is known to facilitate cation transport43 and lipid translocation44,45. PDR5 is an important gene which has been used as a model to study the accumulation of drugs into yeast cells and the induction of pleiotropic drug resistance46,47. PDR5 has several homologs, including PDR15 and PDR10. In our dataset, no other genes involved in PDR or vacuolar sequestration were overexpressed. These results suggest that under normal growth conditions, PDR5 is targeted and specific towards cannabidiol.

Mitochondrial regulation and cellular biosynthetic processes are expressed at the growth inflection point

The remaining significant genes in the dataset were TAR1, RDN25-1, RDN37-1, RDN18-1, and IMT1. TAR1 and RDN25-1 are opposite adjacent neighbors in the genome, and their expression is tightly coregulated48. TAR1 is a mitochondrial activity regulator49, and the remaining genes are rRNA’s and the tRNA IMT1. These genes, which are related to biosynthetic processes, remain significant after we controlled for the influence of cellular proliferation using an RNA library collected during normal mid-logarithmic growth phase to compare to our treatment group in our DESeq2 analysis. One possible explanation is that they may be connected to a change in the source of carbon for the cell cultures, consistent with a post-diauxic shift.

Metabolomics corroborates a post-diauxic shift model and identifies unique and overabundant lipids at the growth inflection point

Cannabidiol treatment instigates an intracellular metabolic response in S. cerevisiae CENPK2-1C at the observed phenotypic change in a growth curve

We were interested in metabolomics because we hypothesized that changes in central and secondary metabolism would underpin the cause of the observed post-diauxic shift. Our questions regarding metabolism were if samples within treatment groups were similar to each other, especially if CBD-treated samples shared high similarity; and which compounds were overabundant within CBD-fed samples. Principal Component Analysis (PCA) revealed distinctive shifts in the intracellular metabolite composition of Saccharomyces cerevisiae upon cannabidiol (CBD) treatment. The cell pellets demonstrated a higher degree of clustering influenced by CBD treatment compared to supernatants, suggesting a pronounced intracellular response to CBD feeding. We base this conclusion on the clustering patterns we observed in our PCA, where supernatants clustered tightly regardless of their condition while cell pellets demonstrated separation across conditions (Fig. 4). This PCA showed that 30 principal components explained a variance (R2) of 0.957 indicating a robust data model.

Fig. 4: Visualization of LC-MS(qTOF) untargeted metabolomics data with PCA.
figure 4

PC1 is plotted against PC2. 30 principal components explained a variance (R2) of 0.957.

After PCA, we conducted a two-way Analysis of Variance (ANOVA), which facilitated the discovery and annotation of numerous significant metabolites. Out of the 2798 spectra identified in the entire data set, 417 were found in the solvent blanks, leaving 2381 spectra that we found in our supernatant and cell extracts. We wanted to know which of these were abundant in CBD-fed cell pellets, so we removed the supernatant samples from the dataset and ran our ANOVA. 480 spectra showed significant alpha values in a two-way ANOVA, with condition (CBD-fed, MeOH-fed, unfed, or blank) as the experimental group and substrate (pellet or blank) as the fixed factor. Each condition was tested in biological triplicate (n = 9) with each biological replicate sampled three times for technical replication. We used ANOVA to evaluate our data as it enabled simultaneous comparison across multiple conditions and accounted for variable subsets within each condition and allowed us to identify individual compounds which were important to the intercellular response from cannabidiol. Among the significant compounds, we were especially interested in 48 spectra because of their overabundance in the CBD-fed cell pellets. We determined overabundance based on a comparison of their intensity counts to the other samples. These intensity plots are included in the supplementary information. From these 48 spectra, twenty non-cannabidiol compounds were identified via HRMS parent ions and in-source fragments, all of which were within 0.003 of the theoretical m/z value. Detailed information regarding the retention times, m/z, parent masses (M), and molecular formulas are presented in Table 1. A compound of particular note was 1-docosanoyl-glycero-3-phosphoserine, which was present exclusively in CBD-fed cells (Fig. 6). The identification of 1-docosanoyl-glycero-3-phosphoserine was substantiated by in-source fragments. These metabolites were then used for a KEGG pathway analysis.

Table 1 Compounds marked with a are considered central metabolites and compounds marked with b act as central or secondary metabolites, depending on the pathway in which they participate

Central and secondary metabolites reveal enrichment of pathways related to cell growth

The significant metabolic pathways delineated by our KEGG pathway analysis (Fig. 5) yielded vital insights into the cellular machinations instigated by CBD exposure. We observed upregulation in pathways such as aminoacyl-tRNA biosynthesis; phenylalanine, tyrosine, and tryptophan biosynthesis; valine, leucine, and isoleucine degradation and biosynthesis, as well as ubiquinone and other terpenoid-quinone biosynthesis. The convergence of these pathways, particularly those involved in amino acid metabolism and the biosynthesis of essential coenzymes, suggests a response of primary metabolism to CBD treatment, with potential implications for cellular growth and stress adaptation mechanisms. Overall, an enrichment of tRNA biosynthesis corresponds tightly to transcriptomics dataset, which found significant overexpression of the tRNA initiation factor IMT1. In parallel to the changes in amino acid metabolism and glycerophospholipid metabolism, we can infer that central carbon metabolism does in fact showcase that a cellular growth response takes place at the sampled time point, corroborating the post-diauxic shift hypothesis.

Fig. 5: Comprehensive Analysis of LC-MS(qTOF) metabolomics data reveals key pathways.
figure 5

Asterisks denote pathways of significant alteration. Predominant among these is aminoacyl-tRNA biosynthesis, highlighted by an elevated presence of metabolites including phenylalanine, leucine, tryptophan, methionine, tyrosine, alanine, and histidine. This observation aligns with transcriptomic data (Fig. 3), particularly the upregulation of IMT1. The biosynthesis of branched-chain amino acids (valine, leucine, isoleucine) is also indicated by the notable levels of leucine and isoleucine. Phenylalanine, tyrosine, and tryptophan biosynthesis are further corroborated by their significant presence. While biotin metabolism and glycerophospholipid metabolism did not show significant changes, the detection of 8-Amino-7-oxononanoate and specific lipids (unrepresented in the KEGG database) such as 9-Hydroxy-12-octadecenoic acid and 1-docosanoyl-glycero-3-phosphoserine suggests underlying activity. Detailed data are presented in Supplementary Table 5.

Cannabidiol treated S. cerevisiae CENPK2-1C uniquely produce a 22 carbon monoacylglycerophosphoserine

We report the discovery of high production of 1-docosanoyl-glycero-3-phosphoserine which was exclusive to our CBD-fed cells (Fig. 6). We confirmed the identity of this molecule with four fragment matches we found in Metfrag. Moreover, we report the existence of an additional fragment at 528.28483 m/z (negative mode), which is made obvious by a shared abundance pattern across all samples for this spectrum across all significant spectra we identified as this retention time (see Supplementary Information Fig. 2). Further, when we ran an independent analysis using reverse-phase chromatography to check for cannabidiol metabolites, we replicated the discovery of this compound (however, metabolites of cannabidiol were not identified- see Supplementary Information Fig. 3). Within the cell-pellet analysis, glycerophosphoglycerol and sn-glycero-3-phosphoethanolamine were exciting to identify because it furthers the body of evidence which points towards dynamic lipidomic remodeling in the CBD-treated cells and might be related to the production of 1-docosanoyl-glycero-3-phosphoserine. While the pathway analysis did not reach a significant p-value for the category of glycerophospholipid metabolism, we note that this analysis excluded 9-Hydroxy-12-octadecenoic acid and 1-docosanoyl-glycero-3-phosphoserine from the input, as they were not available in the KEGG database. Accordingly, we conclude that this result is currently inconclusive and warrants further investigation.

Fig. 6: Intensity (y-axis, a.u.) plot of the m/z peak at 580.36028 across all samples (x-axis) in the analysis.
figure 6

(Top) This compound was found only in cell pellets from CBD-treated samples, indicating that cannabidiol initiates the production of this compound in vivo. (Bottom) A visual description of 1-docosanoyl-glycero-3-phosphoserine, the identity of the compound associated with the M. 581.3676 (see Table 1). 1-docosanoyl-glycero-3-phosphoserine is a monoacylglycerophophatidylserine with 22 carbons in the alkyl chain. We predict it is associated with the plasma membrane because it is a glycerophospholipid, a group of lipids commonly associated with plasma membrane construction.

Knock-out and complementation of PDR5 results in deletion-and-rescue of the post-diauxic shift Phenotype

To further underpin the contribution of the overexpressed transcripts to the observed phenotypic changes in our study, we performed a knock-out and complementation assay on PDR5 and a knockout assay of CIS1. CIS1 was deleted and selected for using a uracil marker. PDR5 was knocked out using a kanamycin marker, or for complementation was knocked out and selected for using a PDR5-kanamycin insertion cassette.

In continual biomass monitoring experiments, Δcis1::URA3 behaved as the wildtype, demonstrating that the post-diauxic shift phenotype occurs independently of the presence of a functional copy of CIS1 in the genome. However, Δpdr5::KAN lost the associated post-diauxic shift and Δpdr5::PDR5-KAN regained the function. Data are shown in Supplementary Information Fig. 63.