Bakuchiol and Ethyl (Linoleate/Oleate) Synergistically Modulate Endocannabinoid Tone in Keratinocytes and Repress Inflammatory Pathway mRNAs

BAK + ELN is a potent FABP5 and FAAH inhibitor

We used in vitro assays to quantity the inhibitory activity of compounds against FABP3, FABP5, and FAAH (Figure 2). These studies were performed using a mixture of fatty acid esters, ELN, consisting of approximately 70% EL and 15% ethyl oleate. All compounds showed FABP3 inhibitory capacity, although average half-maximal inhibitory concentrations (IC₅₀) were significantly lower for CBD and BAK than for ELN (P < 0.05) (Figure 2a‒d). Consistent with docking simulations (Figure 1), BAK and BAK + ELN had significantly lower FABP5 inhibition IC₅₀ values than CBD and ELN (P < 0.05) (Figure 2e‒h). FABP5 inhibition IC₅₀ values were lowest for BAK + ELN (Figure 2e‒g). BAK, ELN, and BAK + ELN were all more effective FAAH inhibitors than CBD (Figure 2l). BAK + ELN was the most potent FAAH inhibitor (Figure 2i‒k), although IC₅₀ values did not differ significantly among BAK, ELN, and BAK + ELN (Figure 2l). The combination of BAK + ELN with CBD did not increase FAAH inhibition potency (Figure 3a‒d) nor did the combination of BAK + ELN with hemp seed oil (Figure 3e‒h).

BAK + ELN has a stronger anti-inflammatory transcriptional effect than CBD, BAK, and ELN

Microarrays were used to evaluate the direct effects of each compound on gene expression in TNF-stimulated KCs. Treatment of KCs with TNF altered the expression of 2,293 genes (1,137 increased, 1,156 decreased; false discovery rate < 0.10 with fold change [FC] > 1.25 or FC < 0.80) (Figure 4a and b). In TNF-stimulated cells, the relative effects of CBD, BAK, ELN, and BAK + ELN were smaller but significant. CBD altered the expression of 319 genes (205 increased, 114 decreased) (Figure 4e and f), whereas BAK + ELN altered the expression of 490 genes (202 increased, 288 decreased) (Figure 4q and r). BAK and ELN individually had smaller effects on gene expression (Figure 4i, j, m, and n). These effect size differences were apparent from changes in principal component (PC) scores (Figure 5a‒d). CBD and BAK + ELN had strong effects on PC scores, with a large shift in the PC centroid (Figure 5a and d), whereas BAK and ELN had smaller effects on the PC centroid (Figure 5b and c).

The effects of CBD correlated weakly with those of TNF (r_s = ‒0.02) (Figure 4, Figure 5g). In contrast, the effects of BAK and ELN were negatively correlated with those of TNF (r_s = ‒0.32 and ‒0.20, respectively) (Figure 5e and g). However, this negative correlation was even stronger for the BAK + ELN treatment (r_s = ‒0.45) (Figure 5e and g). Inspection of self-organizing maps confirmed a disparity between the effects of CBD and BAK + ELN (Figure 5f and h).

CBD upregulates metallothionein gene expression and downregulates genes mediating alcohol, steroid, and cholesterol synthesis

CBD altered the expression of 319 genes but did not reverse changes seen with TNF treatment (Figure 6e). Genes increased by CBD included CCL8, ATF3, SLC3A2, and PLIN2 (Figure 6b, c, and e). CBD increased metallothionein family genes associated with stress response to metal ions (Figure 6g) as well as genes associated with cytokine receptor interaction (Figure 6i). Genes decreased by CBD included FABP3, HMGCS1, DHCR7, and TMEM97 (Figure 6b, d, and f). Such genes were strongly associated with the synthesis of lipids, steroids, and alcohols (Figure 6h and j).

BAK + ELN upregulates proliferative and DNA repair pathways while downregulating type I IFN responses

BAK, ELN, and BAK + ELN each altered gene expression in a direction contrary to TNF (Figure 5e). The genes increased by BAK (Figure 7a) were associated with cell cycle checkpoints and DNA replication (Figure 7g). The genes decreased by BAK (Figure 7a and d) overlapped with those decreased by CBD (Figure 7b) and were associated with alcohol, cholesterol, and steroid synthesis (Figure 7h).

BAK + ELN had stronger effects on gene expression than BAK (Figure 7a, d, g, and h) or ELN (Figure 7b, e, i, and j). The genes increased by BAK + ELN (Figure 7c) were associated with DNA replication and cell division (Figure 7k). Conversely, the genes decreased by BAK + ELN included MX1, FABP3, and RSAD2 (Figure 7c), and such genes were associated with type I IFN signaling and viral response (Figure 7l). This downregulation of type I IFN genes was not seen with CBD treatment (Figure 8); however, downregulation of MX1 and OAS1 by both CBD and BAK + ELN was observed on the basis of RT-PCR analysis of an independent sample set (Figure 9).

BAK + ELN represses the expression of genes induced during KC differentiation

Activation of the AEA‒CB1 axis blunts KC differentiation (Maccarrone et al., 2003; Paradisi et al., 2008; Roelandt et al., 2012). Markers of early KC differentiation (DSC1, keratin 1 gene K1, keratin gene K10) were slightly increased by CBD, BAK, and ELN (Figure 10a, e, and i) but more often decreased by BAK + ELN (Figure 10m). We evaluated genes having increased expression in KCs undergoing differentiation in a devitalized human dermis model (GSE52651) (Lopez-Pajares et al., 2015). Genes upregulated after 3‒6 days of differentiation were biased toward decreased expression in TNF-stimulated KCs treated with BAK + ELN (P < 0.05) (Figure 10b, f, j, and n). Consistent with this, genes upregulated on each day of differentiation were biased toward decreased expression in BAK + ELN‒treated KCs (P < 0.05) (Figure 10c, g, k, and o). Finally, genes exhibiting a linear trend toward increasing expression over the 7-day differentiation time course overlapped with genes decreased by BAK + ELN (Figure 10d, h, l, and p).

CBD activates MTF-1 pathway genes, whereas BAK + TNF represses an IFN‒signal transducer and activator of transcription/IFN regulatory factor axis

We evaluated the enrichment of DNA motifs (Swindell et al., 2015) in upstream regions of genes altered by CBD and BAK + ELN. Sequences upstream of CBD-increased genes were enriched with IFN regulatory factor and MTF-1 transcription factor motifs (Figure 11a‒c), whereas sequences upstream of CBD-decreased genes were enriched with motifs interacting with E2F, SP2, and activating protein 1 transcription factors (Figure 11d‒f). We identified eight CBD-increased genes, with MTF-1 recognition sequences clustered near the transcription start site (Figure 11g).

Sequences upstream of BAK + ELN‒increased genes were enriched with motifs interacting with ARID3A and SFT2D1 (Figure 12a‒c), whereas sequences upstream of BAK + ELN‒decreased genes were enriched with signal transducer and activator of transcription and IFN regulatory factor family motifs (Figure 12d‒f). Motifs enriched in BAK + ELN‒decreased promoters shared a core sequence (5-AGTTTCnnTTTC/GAAAnnGAAACT-3) concentrated near the transcription start site of downregulated IFN targets (Figure 12g).

BAK + ELN represses the expression of inflammatory genes elevated in diseased skin

An unhealthy skin signature shared by multiple skin diseases has been reported (Mills et al., 2018). Whereas CBD led to a small but significant 5% average increase in the expression of unhealthy skin signature genes, BAK + ELN had no significant effect (Figure 13a and b). We next evaluated genes with expression nonspecifically altered in psoriasis lesions (i.e., having expression altered in psoriasis as well as multiple other skin conditions) (Swindell et al., 2016). Among such genes, those increased in psoriasis lesions were on average increased by CBD (P < 0.01) (Figure 13c), whereas the genes repressed in psoriasis lesions were repressed by CBD (P < 0.01) (Figure 13e). In contrast, BAK + ELN downregulated the genes nonspecifically elevated in psoriasis lesions (P < 0.01) (Figure 13c and d). BAK + ELN also downregulated coexpressed type I IFN‒regulated genes activated in psoriasis lesions, wounds, and multiple types of skin cancer (Swindell et al., 2013b) (Figure 13f, g, h, and m) and upregulated modules containing cell cycle regulators (Figure 13i‒o).

BAK + ELN represses the expression of proinflammatory eCB signaling genes upregulated by TNF

Genes associated with retrograde eCB signaling (Kyoto Encyclopedia of Genes and Genomes pathway hsa04723) were disproportionately downregulated by TNF compared with those of nontreated control (CTL) KCs (P = 0.0063) (Figure 14a) but were upregulated by BAK + ELN + TNF compared with those of TNF-treated KCs (P < 0.001) (Figure 14e). Otherwise, there was no significant trend to indicate that CBD, BAK, or ELN upregulated or downregulated eCB genes in TNF-treated KCs (P ≥ 0.093) (Figure 14b‒d). Genes upregulated by TNF but downregulated by BAK + ELN included PTGS2/COX2, which belongs to a pathway by which 2-AG is metabolized to proinflammatory prostaglandin glycerol esters (Figure 15). Several genes downregulated by TNF but upregulated by BAK + ELN (Figure 14f, h, and i) were associated with pathways downstream of the CB1 receptor (Figure 15).

Molecular docking simulations

The protein data bank file of FABP5 was retrieved from the protein data bank using the PDB identifier 1B56 (http://www.rcsb.org). Subsequently, protein preparation was achieved using the default workflow of Protein Preparation Wizard implemented in Maestro (Protein Preparation Wizard 2019; Schrodinger, New York, NY). Briefly, all missing hydrogen atoms were added after deleting original hydrogen atoms, bond orders were assigned, and water molecules >5Å away from het-groups were removed, whereas missing loops and side chains were added with Prime (Prime 2019; Schrodinger). Protonation and tautomeric states of amino acids were adjusted to pH 7 ± 2.0 with EpiK. Thereafter, hydrogen bond network was optimized with PROPKA at pH 7.0, followed by a final restrained minimization with a convergence of heavy atoms to a root-mean-square deviation of 0.3 Å. The receptor active site was defined by covering a volume of the binding site through receptor grid generation in Glide (Receptor Grid Generation 2019; Schrodinger) and OPLS3e force field. The area surrounding the cocrystalized ligand was specified as the active site, thus setting the grid box to be the centroid of the cocrystalized ligand.

For ligand preparation, LigPrep module (LigPrep 2019; Schrodinger) was used to generate low-energy three-dimensional structures of BAK, EL, and CBD using the OPLS3e force field. During this process, their ionization states were assigned at pH 7 ± 2.0 using EpiK (EpiK 2019, Schrodinger). For the protein grid generation step, hydrogen bond constraints were applied to critical protein residues on the basis of biochemical data. Molecular docking simulations of each compound at the FABP5-active site were subsequently performed with Glide XP (Ligand Docking 2019, Schrodinger). The best poses of the protein‒ligand complex were selected on the basis of their glide energy and glide model, whereas the free binding energy of selected docked complexes was calculated with Prime Molecular Mechanics-Generalized Borne Surface Area.

FABP3, FABP5, FAAH, and PTGS2 (COX-2) inhibitor activity

FABP3 inhibitory activity was evaluated using a commercial FABP3 ELISA Kit (catalog number ELH-FABP3, Raybiotech, Peachtree Corners, GA). We dissolved 0.1 g of sample in 1 ml DMSO and diluted it with diluent C from the kit to make a stock solution (1 mg/ml). Serial dilution (1:2) of the sample was then made using diluent C to determine the IC₅₀. We added 50 μl of sample and 50 μl of FABP3 (200 ng/ml) (catalog number 230-00037-10, Raybiotech). FABP5 inhibitory activity was evaluated using a commercial FABP5 ELISA Kit (catalog number ELH-FABP5, Raybiotech). For the FABP5 assay, 0.1 g of sample was dissolved in 1 ml DMSO and diluted with diluent C from the kit to make a stock solution (800 μg/ml). Serial dilution (1:2) of the sample using diluent C was then made to determine the IC₅₀. We added 50 μl of sample and 50 μl of FABP5 (200 ng/ml) (catalog number 268-10276-1, Raybiotech) to each well. FAAH inhibitory activity was evaluated using a commercial FAAH inhibitor screening assay kit (catalog number 10005196, Cayman Chemical, Ann Arbor, MI). We dissolved 0.01 g of sample in 1 ml DMSO to make a stock solution. Serial dilution (1:5) of the sample using DMSO was then performed to determine the IC₅₀. PTGS2 (COX-2) inhibition was evaluated using a COX (bovine/human) inhibitor screening assay kit (catalog number 560131, Cayman Chemical). Stock solution was made by dissolving 0.01 g of the sample in 1 ml DMSO. Analyses were then performed using 1-to-3 serial sample dilutions. The IC₅₀ for each serial dilution assay was calculated using robust regression as described previously (Swindell et al., 2020).

KC experiments

Neonatal human epidermal KCs (passage 8) were plated in KC growth medium at late subconfluence density (∼220,000 cells/well). We generated 12 samples with treatments corresponding to nontreated CTL KCs (CTL, n = 2), TNF-treated KCs (TNF, n = 2), KCs treated with TNF and CBD (TNF + CBD, n = 2), KCs treated with TNF and BAK (TNF + BAK, n = 2), KCs treated with TNF and ELN (TNF + ELN, n = 2), and KCs treated with TNF and 1:1 (w/w) combination of BAK and ELN (BAK + ELN, n = 2). Cells were incubated with test materials in KC maintenance medium for 24 hours. After the incubation period, RNA was extracted using the Qiagen (Hilden, Germany) RNAeasy Mini Plus kit with QIAcube Connect automated nucleic acid purification instrument. Total RNA was extracted with yields of 10.1‒17 ng/μl per sample (Figure 18a) and 260/280 absorbance ratios between 1.8 and 2.4 for all samples (Figure 18b).

Microarray analysis

Microarray hybridizations were performed using standard protocols by Thermo Fisher Scientific (Waltham, MA). RNA was shipped on dry ice and further processed upon receipt with complementary RNA yields between 50.0 and 65.6 μg (Figure 18c) and cDNA yields between 14.7 and 20.9 μg (Figure 18d). Gene expression was assayed using the Affymetrix Clariom S array, which is a full-genome platform with 300,304 probe features organized into 27,189 probe sets. Visual inspection of microarray pseudoimages did not reveal any spatial artifacts (Figure 18e‒p). Eukaryotic hybridization CTLs (spike CTLs) were detected on each array with the expected ordering of signal intensities in all cases (Cre > BioD > BioC > BioB) (Figure 18q). Likewise, poly-A RNA in vitro‒synthesized labeling CTLs (Bacillus subtilis genes) were detected in each sample and approximated the expected pattern (dap > thr > phe > lys) (Figure 18r). Normalized unscaled standard error (Bolstad et al., 2005) medians were close to one for each array with a similar interquartile range (Figure 18s). Similarly, relative log expression medians (Bolstad et al., 2005) were approximately zero with a similar interquartile range across samples (Figure 18t). Receiver operating curve area under the curve statistics were calculated for each array as a summary measure of signal intensity differences between exonic (positive CTL) and intronic (negative CTL) probes. Area under the curve statistics were close to one for all samples (0.96‒0.97), indicating a good separation between probes targeting intronic and exonic regions (Figure 18u).

Raw CEL files were normalized using the robust multichip average algorithm to generate log₂-scaled signal intensities (Irizarry et al., 2003). Of the 27,189 probe sets, 19,937 were annotated with a human gene symbol (18,088 unique gene symbols). For those genes represented by multiple probe sets, a single representative was chosen to limit redundancy in the analyses. This was done by choosing the probe set with the highest average expression across the 12 microarray samples. This yielded a filtered set of normalized expression values for 18,088 unique protein-coding human genes. A gene was considered to be expressed in a given sample if its expression was above the 20th percentile. Hierarchical cluster analysis revealed good separation between CTL and TNF-treated samples, with BAK + ELN + TNF samples separated from other TNF-treated samples (Figure 18v). When samples were plotted with respect to the first two PC axes, CTL samples were separated from TNF-treated samples, and most TNF-treated samples were grouped together (Figure 18w). There were no statistically significant outliers with respect to the first PC axis (Grubb’s P = 0.09) (Figure 18x).

Differential expression analyses

Differential expression analyses were performed using linear models with empirical Bayes moderation of standard errors (R package: limma, functions: lmFit and eBayes) (Smyth, 2004). A total of five differential expression comparisons were performed (TNF vs. CTL, 13,910 genes; CBD + TNF vs. TNF, 13,897 genes; BAK+ T NF vs. TNF, 13,881 genes; ELN + TNF vs. TNF, 13,908 genes; BAK + ELN + TNF vs. TNF, 13,888 genes). Each comparison only included genes expressed in at least two of the four samples from both treatments with an SD of median-scaled expression above the 5th percentile. Raw P-values from differential expression analyses were adjusted using the Benjamini‒Hochberg method (Benjamini and Hochberg, 1995). Differentially expressed genes were identified on the basis of a false discovery rate threshold of 0.10 with FC > 1.25 or FC < 0.80. Volcano plots showed similar numbers of increased and decreased genes for each comparison (Figure 4a, e, i, m, and q). MA plots showed that FC estimates did not differ systematically on the basis of absolute expression levels (Figure 4b, f, j, n, and r). For most comparisons (except TNF + ELN vs. TNF), quantile‒quantile plots revealed an overabundance of small or large moderated t-statistics (Figure 4c, g, k, o, and s) and left-skew of the raw P-value distribution (Figure 4d, h, l, p, and t), consistent with significant differential expression.

RT-PCR (replication study)

The experiment was repeated using neonatal human epidermal KCs (passage 2) pooled from 10 donors (catalog number FC0064, Lifeline Cell Technology, Frederick, MD). Test materials were stored at room temperature, and stock solutions were prepared on the day of the experiment (20 mg/ml in DMSO). Cells were plated in DermaLife Basal Medium (catalog number LM-0004, Lifeline Cell Technology) in a 24-well plate and allowed to grow to subconfluence. Growth media was replaced by maintenance/differentiation medium containing TNF-α (50 ng/ml), and test compounds were applied for 24 hours (BAK, ELN, CBD, or BAK + ELN). After 24 hours, KCs were observed under an inverted microscope (Nikon Eclipse TS100) and RNA was extracted using a commercial kit (Cat no. 740955.240C, Macherey-Nagel, Bethlehem, PA). Purified total RNA was assessed at 260 and 280 nm using the NanoDrop Lite (ThermoFisher Scientific, Waltham, MA). cDNA was prepared using the AzuraQuant cDNA kit (Azura Genomics, Raynham, MA). The expression of selected genes was evaluated using real-time quantitative PCR with AzuraQuant Green Fast qPCR master mix (Azura Genomics) and the BioRad iCycler iQ Detection System. PCR primers were obtained from Thermo Fisher Scientific (TaqMan assays: HMGCS1, Hs00940425_g1; MX1, Hs00895604_m1; OAS1, Hs00973635_m1; GAPDH, Hs02786624_g1). Relative gene expression was evaluated using the ΔΔCt method (Livak and Schmittgen, 2001) with GAPDH used as a reference (housekeeping) gene.

Motif enrichment analysis

Genes significantly regulated by CBD or BAK + ELN were assessed to identify DNA motifs enriched in sequences upstream of transcription start sites. The analysis was performed using a precompiled dictionary of 2,935 experimentally determined motifs associated with human and/or mouse transcription factors or unconventional DNA binding proteins (Swindell et al., 2015). Motif enrichment was evaluated using semiparametric generalized additive logistic models (Swindell et al., 2013a). For each motif and set of differentially expressed genes, this approach generates a Z statistic, which is proportional to the degree of motif over-representation in sequences upstream of differentially expressed genes compared with that in sequences upstream of KC-expressed nondifferentially expressed genes (Swindell et al., 2013a). To CTL the false discovery rate among tests performed for the 2,935 motifs, raw P-values associated with Z statistics were adjusted using the Benjamini‒Hochberg method (Benjamini and Hochberg, 1995).

Cortisol and cytokine levels in RHE-conditioned medium

We evaluated the effects of CBD, BAK, and BAK + ELN on cortisol and IL-8 production in an RHE-conditioned medium. RHE tissues were obtained from Zen-Bio (catalog number RHE-24, order number 51500, Durham, NC) and stored at 4 ^oC. The following day, tissues were transferred to six-well plates and equilibrated for 4 hours in RHE medium from Zen-Bio (ZenSkin RHE Assay Medium, lot number 031620). The medium was changed before the application of test materials. Test materials were solubilized in caprylic acid/caprylic triglycerides and spread evenly on top of tissues (2 mg/cm²) using a positive displacement pipette. Solvent alone was used as the negative CTL treatment (CTL). After a 48-hour incubation period, cortisol output in a conditioned medium was measured using the Cortisol Parameter Assay Kit (catalog number KGE008B, R&D Systems, Minneapolis, MN). Tissue viability was assessed using the MTT technique following the manufacturer’s instructions and reagents. IL-8 output into the conditioned medium was measured using a sandwich ELISA assay with the antibodies from Invitrogen (catalog number CHC1303, Thermo Fisher Scientific) and BioLegend (catalog number 501101 and 501201), respectively. Colorimetric measurements were performed using a high-performance spectrophotometer (SpectraMax 190 Microplate Reader, Molecular Devices, San Jose, CA) with SoftMax3.1.2PRO software.

A linear standard curve was used to estimate cortisol level (ng/ml) from raw quantification signals. The cortisol level was multiplied by 10 to account for sample dilution, and the signal arising from blank wells was subtracted. To adjust for cell viability, the estimated cortisol level was divided by signals corresponding to the formazan conversion product in MTT assays. IL-8 levels were estimated directly from signal absorbance (450 nm). All estimates were expressed as a percentage of the CTL (solvent-only) treatment. Statistical significance was assessed using Fisher’s least significant difference test with a type I error rate of 0.05.

We thank the National Research Foundation of South Africa for a Competitive Grant for unrated Researchers (grant number 121276) and the Centre for High-Performance Computing (Cape Town, South Africa)) for permitting the use of supercomputing facilities. The authors would like to thank Tony Chang and Marsha Sintara (both from International Chemistry Testing, Hopkinton, MA) for technical assistance as well as Nosipho Cele and Paul Awolade (both from School of Chemistry and Physics, University of KwaZulu-Natal, Durban, South Africa) for assistance with molecular docking studies.