Essential oils diversely modulate genome-wide gene expression in human dermal fibroblasts

Abstract The increasing popularity of essential oils for skincare has led to investigation of their biological effects in human skin cells. In this study, we investigated the biological activities of three commercially available essential oils, i.e. rosemary oil, wild orange oil, and a blend (commercial name: Deep Blue) composed of oils from wintergreen, camphor, peppermint, blue tansy, German chamomile, Helichrysum, and Osmanthus, in a pre-inflamed human dermal fibroblast culture model, simulating chronic inflammation. The impact of essential oils on proteins associated with inflammation and tissue remodeling and on the genome-wide expression of 21,224 genes was investigated. The three essential oils diversely modulated global gene expression. Ingenuity Pathway Analysis showed that the oils affected numerous critical genes and signaling pathways. Specifically, rosemary oil influenced processes involved in cancer signaling and metabolism; orange oil affected processes related to cancer signaling, immunomodulation, and metabolism; the blend influenced inflammation, immunomodulation, and wound healing. These findings are largely consistent with the existing literature, supporting the beneficial biological activities of these essential oils. Our study provides the first evidence indicating how these essential oils affect genome-wide gene expression in human skin cells and establishes a basis for further research into their biological mechanisms of action.


ABOUT THE AUTHORS
Our group studies the health benefits of essential oils, particularly the efficacy and safety of essential oils and their active components. Our in vitro and clinical studies of essential oils utilize various experimental approaches, including analytical, biological, biochemical, and biomedical methodologies. The research work discussed herein represents one part of a large research project designed to extensively examine the impact of essential oils on human cells. Our work will further our understanding of the health benefits of essential oils and attract a wide research audience. We believe that a full understanding of these health benefits will ultimately lead to the evaluation and use of essential oils as an adjunctive therapy for a variety of diseases. Dr Han holds a PhD in Biological Sciences and is a Fellow of the American College of Nutrition. Dr Parker holds a PhD in Nutritional Sciences.

PUBLIC INTEREST STATEMENT
Essential oils have been widely used globally owing to their health benefits. Our study examined the effects of three essential oils (rosemary, wild orange, and an essential oil blend) in a human skin disease model of chronic inflammation. The effects of these oils on global gene expression (21,224 genes) were also analyzed. These oils diversely affected numerous genes and pathways, many of which are critically involved in the processes of inflammation, immune responses, tissue remodeling, cancer signaling, and metabolism. These findings suggested that essential oils can be biologically active and impact human gene expression. Additional research focusing on how these oils impact gene expression and biological processes is recommended. A comprehensive exploration of the health benefits of essential oils may lead to viable options for treating many diseases. Thus, this study provides an important basis for further research on essential oils and their health benefits for humans.

Introduction
Scientific research on the biological effects of essential oils in human cells is scarce. Given the common use of essential oils for skincare, scientific studies of their effects on human skin cells are needed. In this study, we investigated the biological activities of three commercially available essential oils, rosemary essential oil (REO), wild orange essential oil (WEO), and an essential oil blend (EOB, commercial name Deep Blue) in human dermal fibroblasts. REO is commonly used as a food flavoring agent. In folk medicine, it is valued for its antispasmodic, analgesic, antirheumatic, carminative, cholagogue, diuretic, expectorant, and anti-epileptic properties (de Melo et al., 2011;Takaki et al., 2008). The main components of REO are 1,8-cineole, α-pinene, and camphor. WEO is composed almost entirely of limonene and possesses anti-inflammatory and anticancer properties (Kummer et al., 2013). Finally, EOB is primarily composed of essential oils from wintergreen, camphor, peppermint, blue tansy, German chamomile, Helichrysum, and Osmanthus, based on the product label.
Here, we studied the effects of these oils on important protein biomarkers associated with inflammation, immune response, and tissue remodeling. We then analyzed the impact of these oils on genome-wide gene expression. This study provides the first evidence of the effects of these essential oils on global gene expression in human skin cells and establishes a basis for further research on the biological mechanisms underlying the action of these essential oils.

Materials and methods
All experiments were conducted using a BioMAP system HDF3CGF, which was designed to model the pathology of chronic inflammation in a robust and reproducible manner. The system comprises three components: a cell type, stimuli to create the disease environment, and a set of biomarker (protein) readouts to examine how the treatments affected the disease environment (Berg et al., 2010).

Cell culture
Primary human neonatal fibroblasts were prepared as previously described (Bergamini et al., 2012) and were plated under low serum conditions for 24 h before stimulation with a mixture of interleukin (IL)-1β, tumor necrosis factor (TNF)-α, interferon (IFN)-γ, basic fibroblast growth factor (bFGF), epidermal growth factor (EGF), and platelet-derived growth factor (PDGF). The cell culture and stimulation conditions for the HDF3CGF assays have been described in detail elsewhere and were performed in a 96-well plate (Bergamini et al., 2012; R Development Core Team, 2011).

Protein-based readouts
A direct enzyme-linked immunosorbent assay (ELISA) was used to measure the biomarker levels of cell-associated and cell membrane targets. Soluble factors in the supernatants were quantified using either homogeneous time-resolved fluorescence detection, bead-based multiplex immunoassay, or capture ELISA. The effects of the test agents on cell proliferation and viability (cytotoxicity) were measured using the sulforhodamine B (SRB) assay. For proliferation assays, the cells were cultured and measured after 72 h, which is optimal for the HDF3CGF system, and the detailed procedure has been described in a previous study (Bergamini et al., 2012). Measurements were performed in triplicate wells, and a glossary of the biomarkers used in this study is provided in Supplementary Table S1.
Quantitative biomarker data are presented as the mean log 10 relative expression level (compared to the respective mean vehicle control value) ± standard deviation of triplicate measurements. Differences in biomarker levels between essential oil-and vehicle-treated cultures were tested for significance with the unpaired Student's t test. A p value <0.05, outside of the significance envelope, with an effect size of at least 10% (more than 0.05 log 10 ratio units), was regarded as statistically significant.

RNA isolation
Total RNA was isolated from cell lysates using the Zymo Quick-RNA MiniPrep kit (Zymo Research Corp., Irvine, CA, USA) according to the manufacturer's instructions. RNA concentration was determined using a NanoDrop ND-2000 system (Thermo Fisher Scientific, Waltham, MA, USA). RNA quality was assessed using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) and an Agilent RNA 6000 Nano kit. All samples had an A260/A280 ratio between 1.9 and 2.1 and an RNA integrity number score >8.0.

Microarray analysis of genome-wide gene expression
The effect of essential oils on the expression of 21,224 genes was evaluated in the HDF3CGF system after a 24 h treatment. Samples for microarray analysis were processed by Asuragen, Inc. (Austin, TX, USA) according to the company's standard operating procedures. Biotin-labeled cRNA was prepared from 200 ng of total RNA using an Illumina TotalPrep RNA Amplification kit (Thermo Fisher Scientific) and one round of amplification. The cRNA yields were quantified using ultraviolet spectrophotometry, and the distribution of the transcript sizes was assessed using the Agilent Bioanalyzer 2100. Labeled cRNA (750 ng) was used to probe Illumina human HT-12 v4 expression bead chips (Illumina, Inc., San Diego, CA, USA). Hybridization, washing, staining with streptavidin-conjugated cyanine-3, and scanning of the Illumina arrays were carried out according to the manufacturer's instructions. The Illumina BeadScan software was used to produce the data files for each array; the raw data were extracted using Illumina BeadStudio software.
The raw data were uploaded into R (R Development Core Team, 2011) and analyzed for qualitycontrol metrics using the beadarray package (Dunning, Smith, Ritchie, & Tavare, 2007). The data were normalized using quantile normalization (Bolstad, Irizarry, Astrand, & Speed, 2003), and then re-annotated and filtered to remove probes that were non-specific or mapped to intronic or intragenic regions (Barbosa-Morais et al., 2010). The remaining probe sets comprised the data-set for the remainder of the analysis. The fold-change expression for each set was calculated as the log 2 ratio of essential oils to the vehicle control. These fold-change values were uploaded onto ingenuity pathway analysis program (IPA, Qiagen, Redwood City, CA, USA, www.qiagen.com/ingenuity) to generate the networks and pathway analyses.

Reagents
Essential oils (dōTERRA Intl., Pleasant Grove, UT, USA) were diluted in dimethyl sulfoxide (DMSO) to 8 × the specified concentrations (final DMSO concentration in culture media was no more than 0.1% [v/v]). Then, 25 μL of each 8 × solution was added to the cell culture to obtain a final volume of 200 μL, and DMSO (0.1%) served as the vehicle control.

Bioactivity profile of essential oils in pre-inflamed human dermal fibroblasts
We analyzed the activities of essential oils in a dermal fibroblast model, HDF3CGF, which features the microenvironment of inflamed human skin cells presenting with elevated inflammation and immune responses. Cell viability was assessed after treatment with four concentrations of essential oils (0.01, 0.0033, 0.0011, and 0.0004%, v/v). None of the tested concentrations were cytotoxic to these cells, and thus, the high concentration (0.01%) was applied for additional analyses. We next evaluated the effects of essential oils on 17 biomarkers. We found that the essential oils significantly affected the expression of certain biomarkers if the biomarker values were significantly different (p < 0.01) from those in cells treated with vehicle (control), with an effect size of at least 20% (more than 0.1 log ratio units). None of the three studied essential oils significantly affected any of the 17 studied biomarkers at a concentration of up to 0.01% (data not shown). The protein expression data were not consistent with previously reported activities of these essential oils or their major active components. This lack of effect may be partially explained by the concentration (0.01%) chosen to avoid cytotoxicity on fibroblasts, which was lower than that used in previous studies.

Effects of essential oils on genome-wide gene expression
We then analyzed the effects of essential oils at a concentration of 0.01% on the RNA expression of 21,224 genes in the HDF3CGF model. The results showed different effects of essential oils on the expression of human genes. Specifically, among the 130 most highly regulated genes (with a fold change ratio of expression over vehicle control of ≥ 1.5) by REO, most genes (79 of 130 genes) were significantly upregulated, and the remaining genes were downregulated (Table S2). Among the top 104 most highly regulated genes by WEO, over half (58 of 104 genes) were significantly downregulated, and the remaining genes were upregulated (Table S3). Among the 129 most highly regulated genes by EOB, most genes (99 of 129 genes) were significantly upregulated, and the remaining genes were significantly downregulated (Table S4).
IPA studies showed that the bioactivity of these essential oils significantly matched with many canonical signaling pathways from the literature-validated database (Figures 1-3). REO mainly affected Notes: The p-value was calculated using right-tailed Fisher's exact tests. The p-value measures how likely the observed association between a specific pathway and the data-set would be if it was only due to random chance. The smaller the p-value (larger value −ln [p-value], indicated by black bars) of the pathway is, the more significantly it matches the bioactivity of REO. The ratio, indicated using the gray bar, was calculated by dividing the number of genes from the REO data-set that participate in the canonical pathway by the total number of genes in that pathway. NAD, nicotinamide adenine dinucleotide; MAPK, mitogenactivated protein kinase.
pathways and genes associated with processes of cancer biology and metabolism, among others ( Figure 1). For example, the four most well-matched pathways for REO were the mitotic roles of the polo-like kinase pathway, glycolysis I pathway, serotonin degradation pathway, and cell cycle G 2 /M DNA damage checkpoint regulation pathway. These findings suggested that REO may play important roles in modulating these signaling pathways and thus may regulate processes involved in cancer biology.
Interestingly, inhalation of REO decreases the level of the stress hormone cortisol in humans (Atsumi & Tonosaki, 2007), suggesting an anxiolytic effect. REO possesses analgesic effects, primarily due to its high content of eucalyptol, camphor, and alpha-pinene (Raskovic et al., 2015).
WEO modulated pathways and genes associated with cancer signaling, immune modulation, and metabolism ( Figure 2). The four most well-matched pathways for WEO were the pyridoxal 5′-phosphate (PLP) salvage pathway, p53 signaling pathway, tumor necrosis factor receptor 1 (TNFR1) signaling pathway, and salvage pathways of pyrimidine ribonucleotides. These results indicate that WEO may affect cancer biology, immune responses, and metabolism.
WEO is composed almost entirely of limonene (>90%). In vitro and animal studies suggested that d-limonene possesses anti-inflammatory (d'Alessio et al., 2013;Kummer et al., 2013) and anticancer properties (Crowell, 1999;Chaudhary, Siddiqui, Athar, & Alam, 2012;Jia et al., 2013). In vitro, Notes: The p-value was calculated using right-tailed Fisher's exact tests. The p-value measures how likely the observed association between a specific pathway and the data-set would be if it was only due to random chance. The smaller the p-value (larger value -ln [p-value], indicated by black bars) of the pathway is, the more significantly it matches the bioactivity of WEO. The ratio, indicated using the gray bar, was calculated by dividing the number of genes from the WEO data-set that participate in the canonical pathway by the total number of genes in that pathway. HIV1, human immunodeficiency virus type 1.
WEO inhibits cell proliferation and induces apoptosis of colon cancer cells in a concentration-dependent manner (Chidambara Murthy, Jayaprakasha, & Patil, 2012). Another study indicated that limonene may enhance immune responses and act as an immunomodulator (Raphael & Kuttan, 2003). Interestingly, a recent study found that limonene induces the brown fat-like phenotype in 3T3-L1 white adipocytes, suggesting its modulatory role on lipid metabolism, and may have some potential therapeutic effects on metabolic conditions such as obesity (Lone & Yun, 2016). EOB mainly affects pathways and genes associated with inflammation, immunomodulation, and wound healing (Figure 3). The four most well-matched pathways for EOB were LXR/RXR activation, graftversus-host disease signaling, T-cell and B-cell signaling in rheumatoid arthritis, and RAR activation. The major active components of EOB, including methyl salicylate, menthol, eucalyptol, and camphor, possess anti-inflammatory and analgesic properties and are widely used in over-the-counter drugs.
There are several limitations to the current study. First, the in vitro study results cannot be directly applied to more complex human systems. Furthermore, the impact of the studied essential oils on global gene expression was only evaluated after short-term exposure. How these essential oils affect genome-wide gene expression after long-term application remains unclear.
In summary, the gene expression data obtained in this study are largely consistent with the existing literature. The study provides important evidence regarding the potential mechanisms through which these essential oils impact global gene expression in human skin cells for the first time. Moreover, our data suggested that these essential oils possess promising therapeutic benefits, warranting further investigation. Thus, additional careful research into the biological mechanisms of these essential oils is also recommended. Notes: The p-value was calculated using right-tailed Fisher's exact tests. The p-value measures how likely the observed association between a specific pathway and the data-set would be if it was only due to random chance. The smaller the p-value (larger value -ln [p-value], indicated by black bars) of the pathway is, the more significantly it matches the bioactivity of EOB. The ratio, indicated using the gray bar, was calculated by dividing the number of genes from the EOB data-set that participate in the canonical pathway by the total number of genes in that pathway. LXR, liver X receptor; RXR, retinoid X receptor; RAR, retinoic acid receptor; cAMP, cyclic adenosine monophosphate; Cdc 42, cell division control protein 42; FXR, farnesoid X receptor; CTLA4, cytotoxic T-lymphocyte-associated protein 4.