Transcriptomic analysis of tobacco-flavored E-cigarette and menthol-flavored E-cigarette exposure in the human middle ear

Electronic cigarettes (e-cigarettes) are the most widely used electronic nicotine delivery systems and are designed to imitate smoking and aid in smoking cessation. Although the number of e-cigarette users is increasing rapidly, especially among young adults and adolescents, the potential health impacts and biologic effects of e-cigarettes still need to be elucidated. Our previous study demonstrated the cytotoxic effects of electronic liquids (e-liquids) in a human middle ear epithelial cell (HMEEC-1) line, which were affected by the manufacturer and flavoring agents regardless of the presence of nicotine. In this study, we aimed to evaluate the gene expression profile and identify potential molecular modulator genes and pathways in HMEEC-1 exposed to two different e-liquids (tobacco- and menthol-flavored). HMEEC-1 was exposed to e-liquids, and RNA sequencing, functional analysis, and pathway analysis were conducted to identify the resultant transcriptomic changes. A total of 843 genes were differentially expressed following exposure to the tobacco-flavored e-liquid, among which 262 genes were upregulated and 581 were downregulated. Upon exposure to the menthol-flavored e-liquid, a total of 589 genes were differentially expressed, among which 228 genes were upregulated and 361 were downregulated. Among the signaling pathways associated with the differentially expressed genes mediated by tobacco-flavored e-liquid exposure, several key molecular genes were identified, including IL6 (interleukin 6), PTGS2 (prostaglandin-endoperoxide synthase 2), CXCL8 (C-X-C motif chemokine ligand 8), JUN (Jun proto-oncogene), FOS (Fos proto-oncogene), and TP53 (tumor protein 53). Under menthol-flavored e-liquid treatment, MMP9 (matrix metallopeptidase 9), PTGS2 (prostaglandin-endoperoxide synthase 2), MYC (MYC proto-oncogene, bHLH transcription factor), HMOX1 (heme oxygenase 1), NOS3 (nitric oxide synthase 3), and CAV1 (caveolin 1) were predicted as key genes. In addition, we identified related cellular processes, including inflammatory responses, oxidative stress and carcinogenesis, under exposure to tobacco- and menthol-flavored e-liquids. We identified differentially expressed genes and related cellular processes and gene signaling pathways after e-cigarette exposure in human middle ear cells. These findings may provide useful evidence for understanding the effect of e-cigarette exposure.


Results
The main components of e-liquids are propylene glycol (PG), vegetable glycerin (VG), and flavoring agents, with/without the addition of nicotine. In previous studies, when we analyzed the toxicity of the solvent alone (PG:VG = 5:5) versus tobacco-flavored e-liquids (without nicotine) and menthol-flavored e-liquid (without nicotine), the average IC50 values in the flavored e-liquids were decreased 5 , indicating that flavor-added e-liquids are toxic and that between these two flavors, menthol-flavored e-liquids are more cytotoxic. In the present study, the effect of the components of e-liquids on cell viability was measured again via CCK-8 analysis, and the IC50 values were obtained. The IC50 value of the solvent of the e-liquids (PG:VG = 5:5) was 4.50 ± 0.14%, and that of nicotine was 0.07 ± 0.01 mg/mL. Additionally, the IC50 values of tobacco-flavored e-liquid and menthol-flavored e-liquids were 3.02 ± 0.16% and 1.62 ± 0.01%, respectively (Fig. 1). These IC50 values were subjected to further analysis. The IC50 value of the menthol-flavored e-liquid was lower than that of the tobacco-flavored e-liquid. However, when we identified the expression levels of inflammatory cytokines in HMEEC-1 after e-liquid exposure, the expression levels of COX-2 and TNF-α were found to be higher in the tobacco-flavored e-liquid than in the menthol-flavored e-liquid, indicating that the mechanism underlying the observed cytotoxicity differs between the two flavors (Fig. 2).

E-liquid-related gene expression profile in HMECC-1.
To determine whether gene expression is altered in response to exposure to different flavored e-liquids, the solvent or nicotine, we investigated gene expression levels through mRNA sequencing analysis and generated a heat map using the average values of two replicates for all five samples (Fig. 3A). The heatmap of the expression values of the selected DEGs in log10 (FPKM) units was compared across genes and samples (fold change > 2.0 and Q < 0.05). The differentially expressed genes (DEGs) between the two selected biological conditions were analyzed (Control-VS-Menthol, Control-VS-Nicotine, Control-VS-PV, Control-VS-Tobacco, Menthol-VS-Nicotine, Menthol-VS-PV, Menthol-VS-Tobacco, Nicotine-VS-PV, Nicotine-VS-Tobacco, PV-VS-Tobacco). Since our study was aimed to identify the gene expression according to two different flavored e-liquids, volcano plots were drawn for comparisons across control and two different flavored e-liquids (Fig. 3B). Both tobacco-flavored and menthol-flavored e-liquids induced significant changes in gene expression. Under tobacco-flavored e-liquid exposure, a total of 843 genes were differentially expressed with fold changes of | Log2 |> 1 in response to e-liquid exposure (Q < 0.05) (Fig. 3C, Supplementary dataset S1). Among these genes, 262 were upregulated, and 581 were downregulated (Fig. 3C). Menthol-flavored e-liquid exposure resulted in 228 upregulated genes and 361 downregulated genes with fold changes of | Log2 |> 1 (Q < 0.05) (Fig. 3C, Supplementary dataset S2).
KEGG pathway analysis revealed that tobacco-flavored e-liquid exposure caused a prominent change in the expression of genes involved in translation, infection, cell cycle/apoptosis-related signal transduction, cancer and the regulation of metabolic pathways (Table 1). Menthol-flavored e-liquid exposure resulted in a prominent Scientific Reports | (2020) 10:20799 | https://doi.org/10.1038/s41598-020-77816-2 www.nature.com/scientificreports/ change in the expression of genes involved in metabolic pathways, cancer, inflammation, and focal adhesion ( Table 2). For the tobacco-flavored e-liquid group, the GO annotations of the predicted targets enriched among the 838 genes that were mappable to DAVID were selected according to a |fold-change|≥ 2 and Q-value ≤ 0.05 compared to the control. For the menthol-flavored e-liquid group, the GO annotations of the predicted targets enriched among the 586 genes were also selected according to the above criteria. The functional annotations were categorized into biological processes, cellular components and molecular functions, and only the top 10 GO terms The IC50 values for each type of solution were 3.02 ± 0.16%, 1.62 ± 0.01%, 4.50 ± 0.14% and 0.07 ± 0.01 mg/mL, respectively, and these concentrations were used for further analysis (p < 0.05).

Figure 2.
E-liquids stimulated the expression of inflammatory cytokines in HMEEC-1. Cells were treated with tobacco-flavored e-liquid and menthol-flavored e-liquid at the IC50 concentration for 24 h. Quantitative real-time PCR was performed to evaluate inflammatory cytokine gene expression levels. The expression levels of COX-2 and TNF-α increased when cells were exposed to the e-liquids, and the mRNA expression of cytokines was significantly elevated in the tobacco-flavored e-liquid group (p < 0.05).    (Tables 3, 4). Among the genes affected by the tobacco-flavored e-liquid, the top 20 up-and downregulated genes are listed in Table 5. The upregulated genes included CYP4F3, IL1RL1, PTGS2, and CXCL8, which are related to inflammation; HKDC1, MYEF2 and CYP1A1, which are related to hepatocarcinoma or lung carcinoma; and SLC7A11 and NEFM, which are related to neuronal damage. The downregulated genes included inositol polyphosphate-5-phosphatase D, carbonic anhydrase, proline dehydrogenase, phospholipase, pyrroline-5-carboxylate reductase-like, retinol dehydrogenase 16, and glucokinase. The top 20 genes that were upregulated and downregulated in response to menthol-flavored e-liquids are listed in Table 6. The upregulated genes included IL24, CLU, which are related to apoptosis, and CYP4F3, CCL26, and IL1RL1, which are related to inflammation. The top 20 downregulated genes associated with the menthol-flavored e-liquid included SOX18, which is related to embryonic development; CDH8 and SDK2, which are related to cellular adhesion; and IFITM1, which is related to the immune response.
Direct signaling pathways among DEGs following e-liquid exposure. The molecular signaling networks among the genes that were differentially expressed in response to the tobacco-flavored e-liquids were analyzed to predict the relevant molecular pathways. The direct interaction pathways among the upregulated and downregulated genes are demonstrated in Fig. 4 and the supplementary information (SI). Various genes were related to each other through several key regulator genes. Among the upregulated genes, IL6, PTGS2, FOS, CXCL8, and JUN showed high connectivity in the pathways; among the downregulated genes, TP53 showed high connectivity in the pathways. The molecular signaling networks among the genes that were differentially expressed in response to the menthol-flavored e-liquid were also analyzed to predict the relevant molecular pathways ( Fig. 5 and SI). MMP9 (matrix metallopeptidase 9), PTGS2 (prostaglandin-endoperoxide synthase 2), HMOX1 (heme oxygenase 1), and NOS3 (nitric oxide synthase 3) were upregulated, while MYC (MYC protooncogene, bHLH transcription factor) and CAV1 (caveolin 1) were downregulated. These genes were predicted to be key genes of the direct signaling pathways among the DEGs following menthol-flavored e-liquid exposure.
Cell process-and disease-related biological pathways among genes identified in response to e-liquid exposure. To demonstrate the effect of the tobacco-flavored e-liquid on human middle ear epithelial cells, we analyzed cellular process-and disease-related pathways. We found that the tobacco-flavored e-liquid-related cellular processes included monocyte chemotaxis, leukocyte recruitment, the inflammatory response, the immune response, apoptosis, SMC proliferation, tumor growth, endothelial cell migration, oxidative stress and ROS generation. Several cell processes and diseases related to inflammation (red rectangle), carcinogenesis (blue rectangle), oxidative stress (orange rectangle), lung disease (brown rectangle) and arterial disease (purple rectangle) showed relationships with the biological signaling networks among the DEGs of the tobacco-flavored e-liquid-treated group (Fig. 4). In the menthol-flavored e-liquid-treated group, numerous cell processes and diseases related to cancer (blue rectangle), inflammation (red rectangle), oxidative stress (green Table 2. Analysis of enrichment of KEGG pathways and involved genes against menthol-flavored e-liquids.

qRT-PCR expression levels of the potential biomarkers of e-liquid exposure.
Based on the RNAsequence and pathway analyses, we investigated several genes (IL6, PTGS2, FOS, CXCL8, JUN and TP53) as candidate tobacco-flavored e-liquid-responsive biomarkers. The selection of key genes in the analyzed pathways was performed based on the local connectivity; Pathway Studio provides the number of relationships with neighbors as local connectivity, which is considered as a parameter for scoring the significance of the entity in the pathway. The top six genes that showed high local connectivity values were selected as the key genes in the pathway. To validate the RNA sequencing results, we examined the expressed transcript levels by qRT-PCR. In accordance with the RNA sequencing results, the expression levels of IL6, PTGS2, FOS, CXCL8, and JUN were found to be increased, whereas that of TP53 was found to be decreased (Fig. 6).

Discussion
E-cigarettes are evolving and diversifying, and their use among youth and young adults is a major public health concern. While e-cigarette products and their patterns of use are changing quickly, the associated health effects are not completely understood. Since flavoring agents are reported to contribute to the observed cytotoxicity, among other components of e-liquids [18][19][20][21] , it would improve the understanding of the effect of e-cigarettes on human organs and bodies to study their transcriptomes at the gene level according to exposure to different flavors. In the current study, we identified the genes and pathways showing alterations in expression upon exposure to e-cigarettes. Furthermore, we investigated the direct interactions between genes and related cellular processes Table 3. Go annotation of predicted targets in tobacco-flavored e-liquid. The top 10 most enriched GO terms are listed in terms for biological process, cellular component, and molecular function based on p-values. www.nature.com/scientificreports/ or diseases after exposure to e-cigarettes. Our data demonstrated that e-liquids provoke changes in gene expression in the human middle ear and that there are both similarities and differences in gene expression associated with different flavors.
In the KEGG pathway analysis, both tobacco-flavored and menthol-flavored e-liquid exposure was associated with genes that are related to bladder cancer: CXCL8, HBEGF, CDKN2A, FGFR3, THBS1, CDK4, MYC, TYMP, MMP9, and MMP1. Additionally, genes related to small cell lung cancer, such TRAF1, PTGS2, COL4A2, LAMA5, LAMC3, LAMB3, RXRA, TP53, CDK4, MYC, and PIK3R2, were identified under exposure to both flavors. In the menthol-flavored e-liquid group, genes that are linked with ECM-receptor interactions in small cell lung cancer (COL4A2, LAMB3, LAMA5, and LAMC3) and genes associated with proteoglycans in cancer (CAV1, WNT10B, TFAP4, MMP9, HSPG2, FLNC, ITPR1, CTSL, SMO, ANK1, PLCG2, HBEGF, THBS1, MYC, and PIK3R2) were identified. In the network analysis, diseases such as non-small-cell lung cancer, lung cancer, malignant transformation and prostate neoplasm were associated with exposure to tobacco-flavored e-cigarettes. Moreover, when cells were exposed to menthol-flavored e-liquid, cellular processes related to carcinogenesis was enhanced compared to the results of exposure to the tobacco-flavored e-liquid. Currently, there is no strong evidence or studies suggesting a role of electronic cigarettes in the pathogenesis of cancer 3,22,23 . One in vitro study indicated that electronic cigarette vapors induced double-strand breaks in DNA in head and neck squamous cell carcinoma cell lines that, if unrepaired, would result in chromosomal rearrangement and carcinogenesis 24 . Considering the relatively brief history of these products compared to the prolonged progression of carcinogenesis, it would be premature to discuss the long-term effect of electronic cigarette use. Nevertheless, a transcriptomic analysis can help to identify the relevant effects by revealing changes in gene expression levels. www.nature.com/scientificreports/ In the current study, when HMEEC-1 were exposed to tobacco-flavored e-liquid, the top twenty upregulated genes included CYP4F3, IL1RL1, PTGS2, and CXCL8, which are related to inflammation. When cells were treated with menthol-flavored e-liquid, genes associated with inflammation (CYP4F3, CCL26, and IL1RL1) were upregulated, but genes related to apoptosis (IL24, CLU) also showed increases. When we analyzed the molecular signaling network among genes that were differentially expressed in response to both tobacco-flavored and menthol-flavored e-liquids, both groups were found to be enriched in cellular processes such as inflammation, carcinogenesis, and oxidative stress. However, we found that the specific genes involved were different; IL6, PTGS2, FOS, CXCL8, JUN and TP53 were the main key regulator genes in the tobacco-flavored e-liquid group, while MMP9, PTGS2, HMOX1, NOS3, MYC, and CAV1 were the main regulator genes in the menthol-flavored e-liquid group. These findings suggest that e-liquids affect various activities and that different flavors impact numerous pathways.
In the pathway analysis, interstitial fibrosis was a newly predicted pathway among the cellular processes associated with menthol-flavored e-liquid exposure. Additionally, in the KEGG analysis, unlike tobacco-flavored e-liquid exposure, menthol-flavored e-liquid exposure was related to ECM-receptor interactions and focal www.nature.com/scientificreports/ adhesions. Cell proliferation is subject to many levels of control, but it is becoming apparent that mechanical signaling through the cytoskeleton link between focal adhesions and regulators of cellular contractility contributes to the regulation of cell proliferation 25,26 . These findings might explain the difference in cell viability between tobacco-versus menthol-flavored e-liquid exposure. E-cigarettes might be less cytotoxic than conventional cigarettes [27][28][29][30] ; however, we cannot overlook the unexpected health effects resulting from e-cigarette exposure. Other researchers have reported incompatible findings regarding the significant effects of e-cigarettes found in transcriptomic analysis. In studies with HBE cells (BEAS-2B), some authors have reported that e-cigarette vapor (e-vapor) exposure results in lower toxicity than conventional cigarette smoking and that no significant differences in gene expression could be found 27 . Another study reported that e-vapors are not benign and that they elicit discrete transcriptomic signatures in HBE cells, such as alterations in phospholipid and fatty acid triacylglycerol metabolism pathways, cytochrome P450 function, retinoid metabolism, and nicotine catabolism 17 . In studies conducted in mouse heart tissue and a human gingival epithelial organotypic culture system with Tobacco Heating System 2.2 (THS2.2), which heats tobacco instead of burning it and generates an aerosol similar to that of an e-cigarette, the authors reported no changes www.nature.com/scientificreports/ or fewer changes in gene expression compared to conventional cigarette exposure 17,30 . However, these studies adopted different doses of e-cigarette exposure, analysis time points and delivery methods, and they mostly compared their data with the results of conventional cigarette exposure.
In the present study, both tobacco-flavored and menthol-flavored e-liquids affected KEGG pathways involved in glucose, amino acid and lipid metabolism, which could be related to the regulation of many physiological processes. Additionally, the solvent (PG, VG) itself showed cytotoxicity in the CCK analysis cell viability test and induced differential gene expression compared to the control according to the heatmap and the volcano plot. Another report showed that the e-cigarette exposure altered the expression of rhythmic genes in a way that could be translated into systemic biological alterations and that the major solvents used in e-cigarettes (PG, VG) had unsuspected effects on gene expression related to the molecular clock 31 .
Our experiment was performed in the HMEEC line to identify the effect of e-cigarettes on gene expression, related cellular processes and gene signaling pathways. One limitation of this study is that we did not identify the effect of the vaporized e-liquid, so the data may not reflect the changes that occur in human tissues. However, the results of the present study will help future studies examine the mechanisms underlying the effects of e-cigarettes on the human middle ear and upper airway.
To study the effects of the e-liquids, the cells were grown to 60% confluence in 96-well culture plates (SPL Life Science, Korea) at 37 °C in a carbon dioxide-enriched (95% air, 5% CO 2 ) humidified atmosphere. The HMEEC-1 cells were plated in 96-well plates (1 × 10 4 cells/well) for cell viability assays or 6-well plates (5 × 10 5 cells/well) for RNA-seq analysis. The next day, when the cells reached 80-90% confluence, the cells were starved with serumfree medium for 2 h, then exposed to the e-liquids and subsequently incubated for 24 h. The control groups were not exposed to the e-liquids or PG/VG treatment.
Total RNA was isolated from the tissue samples using the phenol-based (TRIzol) method. One microgram of total RNA was processed to prepare the mRNA sequencing library using a TruSeq stranded mRNA sample  Transcriptome data analysis. Potential sequencing adapters within the raw reads were trimmed by Skewer ver 0.2.2 33 . The cleaned high-quality reads after the trimming of sequencing adapters were mapped to the human reference genome (EnsEMBL, GRC37) using STAR (ver 2.5) 34 . Since the sequencing libraries were prepared in a strand-specific manner by using Illumina's strand-specific library preparation kit, the strand-specific library option -library-type = fr-firststrand was applied in the mapping process.
To quantify the mapped reads against the reference genome to obtain gene expression values, Cuffquant (Cufflinks ver 2.2.1) was used with the strand-specific library option -library-type = fr-firststrand and other default options. The gene annotation of the reference genome mm10 from the UCSC Genome Browser (https ://genom e.ucsc.edu) in GTF format was used for a gene model assignment, and the expression values were calculated in fragments per kilobase of transcript per million fragments mapped (FPKM) units. The differentially expressed genes (DEGs) between the two selected biological conditions were analyzed with Cuffdiff software in the Cufflinks package with the strand-specific library option-library-type = fr-firststrand and remaining default options (Control-VS-Menthol, Control-VS-Nicotine, Control-VS-PV, Control-VS-Tobacco, Menthol-VS-Nicotine, Menthol-VS-PV, Menthol-VS-Tobacco, Nicotine-VS-PV, Nicotine-VS-Tobacco, PV-VS-Tobacco). To compare the expression profiles among the samples, the normalized expression values of selected DEGs were clustered in an unsupervised manner with in-house R scripts. DEGs were identified according to a significance threshold of a q-value (false discovery rate) < 0.05. To identify the biological functional roles of the differences in gene expression between the compared biological conditions, a gene set overlap test between the analyzed DEGs and the functionally categorized genes (according to GO biological processes, KEGG pathways, and transcription factor binding target gene sets) was performed using the DAVID tool 35 . The R program was used for the hierarchical Pathway analysis. The molecular pathways of the DEGs identified by mRNA sequencing were analyzed using Pathway Studio version 12.0 (Elsevier, USA). Based on a text-mining algorithm and the scientific literature, this web-based software provides a schematic output consisting of the biological relationships among the imported gene lists as well as an expanded analysis of relevant cell processes and diseases (PMID: 14594725). All relationships among genes, cell processes and diseases were trimmed based on the number of references (> 10) and the reliability of the reference sentences. The selection of key genes in the analyzed pathways was performed based on the local connectivity; Pathway Studio provides the number of relationships with neighbors as local connectivity, which is considered as a parameter for scoring the significance of the entity in the pathway. The top six genes that showed high local connectivity values were selected as the key genes in the pathway. ). The primers designed for these assays are shown in SI. Total RNA was extracted from each cell group using the RNeasy Mini Kit (Qiagen, Hilden, Germany), and DNase I (Qiagen) was used to remove genomic DNA from the RNA samples during RNA preparation according to the manufacturer's instructions. cDNA was synthesized with Superscript II reverse transcriptase (Bioneer, Inc., Daejeon, Korea). RNA yield and purity were determined with a Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA). The quantitation of mRNA expression was carried out using the ABI Prism 7300 real-time polymerase chain reaction (PCR) system (Applied Biosystems, Foster, CA, USA). PCR amplification was performed in a 20-µL final reaction mixture containing 1 µL of 5 pmol forward and reverse primers, 1 µL of cDNA, and Power SYBR™ Green PCR Master Mix (Life Technologies, Carlsbad, CA, USA). The PCR mixtures were incubated at 95 °C for 15 s and 60 °C for 1 min, followed by amplification for 45 cycles. The target mRNA expression levels were normalized to those of the gene encoding endogenous glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and relative gene expression in the experimental groups was calculated via the 2 (−∆∆Ct) method.

Statistical analysis.
All values are represented as the mean ± SD. For data analysis, we used the SPSS 24.0 statistical program. The Kruskal-Wallis test was employed for comparisons between two groups, and ANOVA was used to compare multiple groups in the cell viability assay. A p value of < 0.05 was considered statistically significant. For multiple comparisons, Bonferroni correction was performed.