DNA Methylation Abnormalities and Altered Whole Transcriptome Profiles after Switching from Combustible Tobacco Smoking to Heated Tobacco Products

Heated tobacco products (HTP) are a novel type of tobacco products that allow the inhalation of aerosols containing nicotine and various chemicals produced by electrically heated tobacco leaves (1, 2). HTPs such as Ploom by JT (3), IQOS by Philip Morris International (4), and glo by British American Tobacco (5) have been marketed since mid-2014 in several countries including Japan; however, the associated risk of developing tobacco-related diseases remains unknown.

In July 2020, the FDA authorized the marketing of the IQOS Tobacco Heating System with “Reduced Exposure” due to its significant reduction in harmful and potentially harmful chemical production. Despite insisting that switching completely from conventional cigarettes to the IQOS system could reduce the risk of tobacco‐related diseases and that presents lower risks than the continuing smoking of conventional cigarettes, the Modified Risk Tobacco Products assessment claim was denied because of a lack of sufficient and direct clinical or epidemiologic evidence (6). In recent years, analyses of heated tobacco ingredients have revealed that HTP contain nicotine, at levels similar to those of combustible tobacco, and also carcinogens, although in lower concentrations than combustible tobacco (7). Nevertheless, in Japan, the use of heated tobacco products rapidly spread, mainly because, unlike in the United States, they were approved for sale without performing any risk assessment (8). Moreover, a claim made in the media about HTP being less harmful than combustible cigarettes, despite the lack of scientific evidence, could have also contribute to its spread (8). In addition, the HTPs are sold in stylishly designed packaging to appear harmless, and despite the small print regarding the harmfulness of tobacco products, it has been found that approximately 90% of tobacco product users do not read these hazard warnings (9). Because the number of HTP users has increased dramatically in Asian countries, such as Japan (8, 10, 11) and South Korea (12, 13), as well as the United States and other Western countries, independent epidemiologic studies regarding the health risks of HTPs are urgently required.

Recent improvements in next-generation sequencing (NGS) have enabled inexpensive and rapid analysis of comprehensive molecular profiling, including the genome, DNA methylome, and transcriptome (14). Although the genomic information remains essentially unchanged throughout life, DNA methylation (DNAm) and gene expression profiles can be altered by the lifestyle and exposure to environmental chemicals. DNAm plays important roles in mammalian embryonic development and displays time-, cell-, and tissue-specific profiles. Analysis using Infinium arrays (HumanMethylation450 and MethylationEPIC BeadChips) and epigenome-wide association studies (EWAS) have shown that changes in DNAm status caused by exposure to environmental factors (e.g., combustible tobacco smoking, alcohol consumption, and stress) are associated with the development of certain lifestyle diseases including hypertension (15), type 2 diabetes (16), dyslipidemia (17), and various types of cancer (18–21). Tobacco-associated reliable DNAm biomarkers include AHRR (22–24), ALPPL2 (22), C14orf48 (25), F2RL3 (21, 24, 26), GFI1 (22), GNG12 (22), GPR15 (26, 27), HUS1 (28), LRRN3 (24, 26), MGAT3 (29), MYO1G (22), NCF4 (30), PRSS23 (25), RARA (22, 25), SLAMF7 (25), TMEM51 (22), and TNXB (31). However, how smoking-associated DNAm biomarkers are affected by HTP usage has yet to be elucidated.

Transcriptome profiles also display cell- and tissue-specific patterns that vary depending on lifestyle and environmental factors. For instance, LRRN3 and GPR15 are consistently hypomethylated and upregulated by cigarette smoke exposure (32). Although transcriptomic studies have examined tobacco smoke and HTPs exposure using three-dimensional airway tissue cultures (33) and laboratory animals (34), the effects of HTPs use on human transcriptomic profiles have not been fully investigated.

Although HTPs are associated with lower levels of hazardous chemical exposure than combustible cigarettes, their health effects remain largely unknown, as HTPs have been only recently available. There are concerns that the distribution of heated cigarettes may expand globally as the FDA approves them for sale in the United States in 2020. Because a long-term epidemiologic assessment of the health effects of HTPs would be time consuming, it is necessary to elucidate the health effects using trans-omics surrogate biomarker approaches with short-term observations.

Herein, we performed a trans-omics, molecular epidemiologic population-based cohort study of DNAm and transcriptome analysis, and reported the effects of HTP exposure on smoking-related DNAm biomarkers as well as gene expression alterations in peripheral blood mononuclear cells (PBMC).

This study was conducted according to the Japanese Ethical Guidelines for Medical and Health Research Involving Human Subjects and was approved by the Ethics Committees of the School of Medicine, Keio University (Tokyo, Japan; #2011-0264 and #2018-0336) and Iwate Medical University (Iwate, Japan; #HG2018-006). Written informed consent was obtained from all participants.

The study workflow is shown in Fig. 1. Participants were enrolled from the Tsuruoka Metabolome Cohort Study (TMCS) managed by Keio University, Tsuruoka City, Yamagata Prefecture, Japan. The baseline TMCS survey was conducted during fiscal years (FY; April–March) 2012–2015 with 11,002 registered individuals ages 35–74 years, who had either participated in a municipal or workplace health checkup, as reported previously (35–37). Moreover, a follow-up study is currently underway for the 2018–2021 FYs. The study described here, was conducted with participants of the 2018 FY and the workplace population was supplemented with additional recruited individuals under 40 years of age from the same workplace (until June 2019), to investigate the smoking habits of young adults after the market introduction of HTPs. On the day of the health checkup, health information, including smoking habits, was obtained face-to-face and blood samples were collected using a BD Vacutainer CPT mononuclear cell isolation vacuum collection tube (8 mL; Becton, Dickinson and Company). The tubes were centrifuged immediately and the PBMCs were separated for DNA methylation and transcriptome analysis using a previously established high-quality PBMC isolation method (38). PBMCs were stored at −80°C until further analysis.

HTP users were defined as individuals who smoked combustible cigarettes in the baseline survey but switched to HTPs by the follow-up survey. Dual users were excluded to focus on the health effects of HTPs alone. The individuals in the three smoking habit groups [never smokers (NS), combustible tobacco smokers (CS) at time of blood sampling, and past smokers (PS)], were matched as close as possible with the HTP users considering parameters as gender, age, drinking habits (current drinker or not) and propensity score. The propensity score was calculated taking into account the following variables: (i) NS group—body mass index (BMI), systolic blood pressure (SBP), non-high-density lipoprotein (non-HDL), hemoglobin A1c (HbA1c), hypertension medication, dyslipidemia medication, and diabetes treatment were used as matching factors; (ii) CS group—number of tobacco products smoked/used per day; (iii) PS group—number of years since quitting smoking and, in the case of HTP users, the number of years since switching to HTP.

The exclusion criteria were as follows: (i) patients with a history of malignant neoplasms or (ii) myocardial infarction, (iii) angina patients treated with catheterization, (iv) patients with a history of stroke, (v) individuals who changed their smoking habits due to illness, and (vi) those who lacked smoking habit information from either of the surveys.

Combustible tobacco smoking–associated DNAm biomarkers located near 17 genes (AHRR, ALPPL2, C14orf48, F2RL3, GFI1, GNG12, GPR15, HUS1, LRRN3, MGAT3, MYO1G, NCF4, PRSS23, RARA, SLAMF7, TMEM51, and TNXB) were selected for DNAm analysis in this study (22–31). The primer sequences are shown in Supplementary Table S1.

Cryopreserved PBMC suspensions were thawed at 22°C–24°C and genomic DNA (gDNA) was extracted using a Maxwell 16 Blood DNA Purification Kit and Maxwell 16 Instrument (Promega). DNA absorbance was measured at 260/280 nm (A_260/280) and 260/230 nm (A_260/230) using a Nanodrop 2000/2000c Spectrophotometer (Thermo Fisher Scientific); and its yield was calculated using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). Fragment size and DNA integrity (DIN) were confirmed using a TapeStation with Genomic DNA ScreenTape and Reagents (Agilent Technologies).

Genomic DNA (500 ng) was treated with bisulfate using the EZ DNA Methylation Kit (Zymo Research) according to the manufacturer's instructions. Then, it was amplified using a PyroMark PCR kit (Qiagen) under the following conditions: denaturation at 95°C for 15 minutes, 45 cycles of 94°C for 30 seconds, 56°C for 30 seconds, and 72°C for 30 seconds, and 72°C for 10 minutes. The PCR products were prepared for pyrosequencing according to the PyroMark Q96 Vacuum Prep Workstation protocol (Qiagen, PyroMark Q96 ID User Manual v5). Pyrosequencing was conducted using the PyroMark Gold Q96 reagent kit (Qiagen) with the PyroMark Q96 ID (Qiagen). The DNAm ratios in combustible tobacco smoking biomarkers were calculated using PyroMark Q96 (v2.5.10) and DNAm levels between different smoking habit users were compared using paired t tests. In addition, statistical analysis of the DNAm analysis was performed with two different correction methods: one for gender, age, BMI, and alcohol consumption, and the other for urinary cotinine only.

Cryopreserved PBMC suspensions were thawed 22°C–24°C and total RNA was extracted using a Maxwell 16 Instrument with a Maxwell 16 simplyRNA Blood Kit (Promega). RNA absorbance was measured at 260/280 nm and 260/230 nm using a Nanodrop 2000/2000c Spectrophotometer (Thermo Fisher Scientific) and the RNA yield with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). RNA integrity number (RIN) values were confirmed using an Agilent 2100 Bioanalyzer Instrument and Agilent RNA 6000 Nano Kit (Agilent Technologies). RNA sequencing (RNA-seq) libraries were constructed using total RNA (200 ng) with a SureSelect Strand-Specific RNA Library Preparation Kit (Agilent Technologies) and Agilent Bravo NGS automation system (Agilent Technologies) according to the manufacturer's instructions.

RNA-seq alignment and gene expression quantification were performed using the GTEx pipeline V8 (39), with slight modifications. Briefly, sequence reads were aligned to a human reference genome (GRCh37) using STAR v2.5.0 (40) with GENCODE gene annotation release 19. Gene expression was quantified with RSEM v1.3.1 software (41). Differentially expressed genes (DEG) were identified using the “DESeq2” v1.26.0 package in R v3.6.1 (42, 43). All possible pairwise comparisons were carried out and count data were normalized count and analyzed, using the log likelihood ratio test. Additional analyses were performed adjusting for sex, age, BMI, and alcohol-drinking status. Comparisons between CS, PS, and HTP were further adjusted for smoking duration (year), and urinary cotinine concentration was adjusted for creatinine. The resultant P values were corrected using Benjamini and Hochberg FDR correction. Genes with FDR-corrected P < 0.05 and fold change >1.5 or <0.67 were defined as DEGs (33).

To visualize the similarity of expression profile of smoking status-related genes between groups, principal component analysis (PCA) was performed on the basis of the differentially expressed genes identified in the aforementioned analyses. For visualization of PCA result, the “factoextra” package in R v3.6.1 (44) was used.

In this study, we collected PBMCs from 2,789 of the 11,002 TMCS participants, including 53 HTP; 1,657 NS; 257 CS; and 822 PS users. One HTP user was excluded from the study because the individual smoked one cigarette per day. The remaining 52 HTP users were matched with 52 individuals from each of the NS, CS, and PS groups, whose basic characteristics are shown in Table 1. The P values of the continuous variables were calculated with ANOVA, and those of the categorical variables with the χ² test. Propensity score matching revealed no significant differences between the HTP users and the other groups, suggesting that they were appropriately matched.

High-quality gDNA and total RNA were extracted from the PBMCs of 208 individuals (n = 52 per group) as described previously (38); however, gDNA extraction failed in several samples due to mechanical problems. Consequently, gDNA was obtained from 204 individuals (n = 51 per group) and RNA from 208 individuals (n = 52 per group), and used for subsequent analysis.

The gDNA and RNA quality assessment results are shown in Supplementary Table S2. All extracted gDNA solutions had almost no RNA, protein, or organic contamination and the mean A_260/280 and A_260/230 ratios were 1.88 and 1.83–2.00, respectively, indicating high quality. The DIN ranged from 8.0 to 8.1, which is sufficient for pyrosequencing. Total RNA contained little impurity contamination, with average A_260/280 and A_260/230 ratios of 2.07 and 2.09–2.14, respectively. The RIN was over 9.8, indicating sufficient quality for sequencing analyses.

To elucidate the DNAm profiles of smoking-related DNAm biomarkers in HTP users, we compared DNAm levels between HTP users and individuals with other smoking statuses. We analyzed AHRR (cg23576855, chr5:373,299; Supplementary Table S3; Fig. 2) and found that its DNAm levels were significantly higher in NS (90.1%) than in HTP users (82.4%, paired t test P = 1.73 × 10⁻¹³). Conversely, they were significantly lower in CS (76.1%) than in HTP users (P = 0.0001), but were similar to those of PS users (81.8%, P = 0.9870). The mean DNAm of chr5:373,315 (AHRR), which was amplified using the same PCR primer set, was significantly higher in NS (87.3%) than in HTP users (77.9%, P = 3.79 × 10⁻¹⁰) and significantly lower in CS (72.0%, P = 8.43 × 10⁻⁵), but not significantly difference was observe between PS and HTP users (77.9%, P = 0.6895; Supplementary Table S3; Fig. 2). Similarly, the mean DNAm of cg05575921 (chr5:373,378) in AHRR was significantly higher in NS (82.6%, P = 1.68 × 10⁻¹⁰) and significantly lower in CS (64.9%, P = 0.0002) than in HTP users (71.7%), with no significant difference between PS (70.3%) and HTP users (P = 0.3728; Supplementary Table S3; Fig. 2). Although NS and CS were not matched in this study, our finding that AHRR hypomethylation was 10% lower in CS than in NS is consistent with previous reports (22, 23), even after excluding outlying DNAm data for CS in sensitivity analysis.

Consistently, F2RL3 (cg03636183), MGAT3 (cg05086879), and RARA (cg17739917) DNAm were significantly higher (P = 9.06 × 10⁻⁶, 0.041, and 2.25 × 10⁻⁸) in NS (79.0%, 85.5%, and 22.7%) and significantly lower (P = 0.0013, 0.0028, and 0.0175) in CS (71.9%, 82.8%, and 15.3%) than in HTP users (75.2%, 84.4% and 17.4%), but similar between PS (74.6%, 84.3%, and 18.2%) and HTP users (P = 0.5043, 0.8548, and 0.3926; Supplementary Table S3; Fig. 2). MYO1G (cg12803068) DNAm was significantly lower (P = 0.0021) in NS (86.2%) and significantly higher (P = 8.94 × 10⁻⁴) in CS (91.5%) when compared with HTP users (89.4%), but it was similar between PS (89.7%) and HTP users (P = 0.6448; Supplementary Table S3; Fig. 2). In LRRN3 (chr7:110733941), the DNAm levels were significantly higher in NS (70.6%, P = 1.59 × 10⁻⁶) than in HTP users (63.3%), whereas those in CS (63.5%, P = 0.8936) and PS (65.9%, P = 0.0759) were similar to those of HTP users (Supplementary Table S3; Fig. 2). The DNAm level of cg11556164, another CpG of the LRRN3 gene analyzed in this study, was significantly (76.4%) lower in CS (74.3%, P = 0.0455) than HTP. In contrast, there was no significant difference in DNAm levels between NS (78.0%, P = 0.1618) and HTP or between PS (75.2%, P = 0.3346) and HTP users (Supplementary Table S3; Fig. 2).

DNAm levels in ALPPL2 (cg21566642), GNG12 (cg25189904), GPR15 (cg19859270), and PRSS23 (cg14391737) were approximately 5% higher in NS (16.3%, 36.5%, 90.2%, and 35.8%; P = 4.51 × 10⁻⁴, 9.40 × 10⁻⁸, 2.90 × 10⁻¹¹, and 0.0334, respectively) than in HTP users, with little variation among individuals (Supplementary Table S3; Fig. 2), but did not differ significantly between CS (11.8%, 29.7%, 84.9%, and, 33.6%; P = 0.4324, 0.1261, 0.1771, and 0.8964, respectively) or PS (13.2%, 30.9%, and 86.2%; P = 0.1664, 0.6070, and 0.5210) and HTP users (12.2%, 31.6%, and 85.7%). In PRSS23, DNAm levels were significantly lower in PS (31.5%, P = 0.0267) than in HTP (33.5%).

The DNAm levels of the C140rf48 (cg16398761), HUS1 (cg10190813), and TMEM51 genes (cg090690072 and cg21913886) were not significantly different between HTP (9.2%, 10.1%, 69.6%, and 87.8%) and the other three smoking status groups. Nevertheless, the DNAm levels of the NS user (10.0%, 10.5%, 71.1%, and 88.5%) were significantly higher than those of the PS (8.9%, 9.8%, 87.4%, and 69.2%) and CS users (8.3%, 9.6%, 68.9%, and 87.6%). While the differences might appear small, they are statistically significant. Furthermore, there was no significant difference in the DNAm levels of GFI1 (cg04535902), NCF4 (cg025322700), and TNXB (cg25596754) among any of the groups. This was true even after adjusting DNAm levels for sex, age, BMI, and alcohol-drinking status (Supplementary Table S3; Fig. 2).

In addition, we corrected the DNAm levels for sex, age, BMI, and alcohol drinking status, which are known to influence DNAm profiles, but the corrected results did not significantly differ from the uncorrected results (Supplementary Fig. S1). On the other hand, the results corrected only for urinary cotinine showed that there were no CpG sites with significantly different methylation levels between the CS and HTP groups. When NS group with little or no cotinine in urine was compared with CS group, only 5 of 29 CpG sites showed significant difference in DNAm levels (Supplementary Fig. S2).

To clarify the differences in transcriptome profiles between individuals from the four smoking-habit groups, we conducted a whole transcriptome analysis using total RNA derived from their PBMCs (n = 208, 52 per group). On average, >21-million reads were obtained per group and the percentage of uniquely mapped reads was >91% (Supplementary Table S4). Because the transcriptomic data obtained were of sufficient quantity and quality, we performed a transcriptome-wide association study (TWAS) of smoking habits. Six pairwise comparisons between user groups (NS vs. HTP, PS vs. HTP, CS vs. HTP, NS vs. PS, PS vs. CS, and NS vs. CS) revealed significant variations in the expression of 95 genes in the smoking habits analyzed (FDR-adjusted P value < 0.05, and fold change >1.5 or <0.67). PCA, without correction by any variables, revealed that the gene expression profile of HTP users was similar to that of PS and in an intermediate position between NS and CS (Fig. 3). The details of the six pairwise analyses are shown in Supplementary Tables S5–S16. In the PCA corrected for the variables sex, age, BMI, drinking weeks, smoking duration, and urinary cotinine concentration adjusted for creatinine (Supplementary Fig. S3), the four groups were plotted in close proximity, and there were no evident differences between the gene expression profiles of the four groups.

A comparison of the gene expression profiles of the NS and HTP users revealed that GPR15 and NAV2 were significantly upregulated in HTP users (FDR-adjusted P value < 0.05, fold change > 1.5; Fig. 4A; Table 2; Supplementary Table S5). Notably, GPR15 expression was significantly higher in PS (Supplementary Fig. S1A; Supplementary Table S11) and CS (Supplementary Fig. S1C; Supplementary Table S15) when compared with the NS group; however, 13 genes (FEZ1, FCGR3B, NXF3, LRRC43, PI3, DAB2IP, CXCR1, DGKK, CMTM2, PODN, VIPR2, CNTNAP3, and CNTNAP3B) were significantly downregulated (FDR-adjusted P value < 0.05, fold change < 0.67; Fig. 4A; Table 2; Supplementary Table S6) in HTP users when compared with those in the NS group. None of these 13 genes were significantly downregulated in the PS group when compared with the NS group (Fig. 4A; Supplementary Table S12), but three of them (NXF3, DGKK, and PODN) were significantly downregulated in CS (Fig. 4E; Supplementary Table S16).

When the gene expression of CS and HTP users was compared, we found that three genes (AHRR, LINC00599, and RP11-713M15.1) were significantly downregulated in HTP users (Fig. 4B, Table 4; Supplementary Table S8). Increased AHRR and LINC00599 expression was observed in the CS group when compared with the PS (Fig. 4E; Supplementary Table S13) and NS groups (Fig. 4F; Supplementary Table S15). Conversely, CYP4F22 and PARD3 were significantly upregulated in HTP users when compared with CS (Fig. 4B; Table 2; Supplementary Table S7), and CYP4F22 was significantly downregulated in CS compared with PS and NS (Fig. 4E and F; Supplementary Tables S14 and S16).

In PS and HTP users' transcriptome profile comparison, no genes showed significant differences. For instance, F2RL3 (also used for DNAm analysis) showed no significant differences in gene expression in any pairwise comparison among the four groups (Figs. 4A–F). Although the expression of some genes changed by approximately 9-fold between PS and HTP users, none had an FDR-adjusted P value < 0.05 (Fig. 4B; Supplementary Table S7 and S8), indicating considerable individual variation in the expression profiles of each gene.

The volcano plots of transcriptome profiles corrected for the variables sex, age, BMI, drinking weeks, smoking duration, and urinary cotinine concentration adjusted for creatinine; and that corrected only for urinary cotinine concentration adjusted for urinary creatinine are shown in Supplementary Fig. S4 and S5. As a result of the corrections performed, and in addition to the fact that no known smoking-related markers were detected, the plot shows a distorted shape that is different from that of the original volcano plot. Therefore, these results suggest that the transcriptome data may have been over-corrected. On the other hand, the results of correction by urinary cotinine concentration adjusted for urinary creatinine were almost the same as the results without the correction, although some genes were out of the DEGs.

HTP use has increased exponentially in Japan since 2016; however, their effects on health remain unknown. Therefore, we conducted a cohort-base study in a Japanese population and analyzed their DNAm and whole transcriptome. This research is one of the first to reveal the molecular genetic effects of HTP exposure on 17 smoking-associated DNAm markers as well as in whole transcriptome profiles of a Japanese population. Among the 29 CpGs present in the 17 genes examined, 17 CpGs in 10 genes (AHRR, ALPPL2, F2RL3, GNG12, GPR15, LRRN3, MGAT3, MYO1G, PRSS23, RARA) were significantly hypomethylated in CS and HTP users, who had switched from combustible tobacco to HTP for an average of 1.7 years, when compared with those in NS; whereas AHRR, F2RL3, and RARA were less hypomethylated in HTP users than in CS. In addition, the transcriptome profiles of HTP users and PS differed from those of the NS and CS groups, being positioned in an intermediate position between the NS and CS profiles. Although AHRR expression was the highest in the CS group, and no difference between HTP, PS, and NS users were detected. Moreover, GPR15 expression was the lowest in NS and no difference between HTP, PS, and CS users were detected. Altogether, these results indicate that HTP users display more DNAm abnormalities than the NS, but to a lesser degree than CS, at two DNAm sites in AHRR and F2RL3. Interestingly, of the 17 genes that were subjected to DNAm analysis, only three (AHRR, GPR15, and LRRN3), showed profiles in the CS and HTP groups that differed from that of the NS group. In recent studies, cis-expression quantitative trait methylation analysis (cis-eQTM) has shown that changes in DNAm levels can affect the expression of adjacent genes (45). Therefore, to investigate how changes in DNAm levels are associated with the expression of adjacent genes, comprehensive DNAm analysis using microarrays and targeted bisulfite sequencing and EWASs, together with transcriptome data, are needed to perform an integrated analysis.

Guida and colleagues observed that some CpG sites revert back to normal within 30 years of smoking cessation, whereas other established DNAm markers of combustible tobacco smoking persist, such as AHRR and F2RL3 (24). However, we found that DNAm levels were significantly higher in HTP users than in CS, even just a few years after switching to HTP. Our results may differ from those of Guida and colleagues for the following reasons: (i) gDNA was obtained from samples with different cell composition (whole blood vs. PBMC); (ii) DNAm analysis was performed in different platforms (microarray vs. pyrosequencing); and (iii) the sample matching methods used were different (meta-analysis vs. propensity score matching analysis). These results suggest that the interpretation of the results of DNAm analysis needs to consider the type of sample, platform, and matching.

We also found that the DNAm status of LRRN3 and GPR15 between HTP and CS users was similar, whereas AHRR and F2RL3 were less hypomethylated in HTP than in CS users (Fig. 2). These findings indicate that DNAm is more likely to change in some regions than in others, supporting a previous report (46). The AHRR protein has a similar molecular structure to the aryl hydrocarbon receptor (AhR), which is intimately involved in dioxin toxicity, and inhibits AhR function by forming a complex with the AhR nuclear translocator (ARNT) that binds to AhR (47). Toxic substances in combustible cigarette smoke, such as polycyclic aromatic hydrocarbons and dioxin, induce AHRR expression and suppress its toxicity (48), while it has been suggested that increased AHRR expression induced by toxic substances in combustible cigarette smoke is caused by the promotion of AHRR hypomethylation (49). Bojesen and colleagues reported that the multifactor-adjusted HRs for the lowest versus. highest AHRR methylation quintile were 4.58 [95% confidence interval (CI), 2.83–7.42] for chronic obstructive pulmonary disease (COPD) exacerbation and 4.87 (95% CI, 2.31–10.3) for lung cancer, indicating that AHRR hypomethylation and gene expression are strongly associated with COPD and lung cancer development (50). The methylation status of three CpG sites in AHRR (chr5: 373,299; 373,315; 373,378) is strongly related to AHRR gene expression, which is in agreement with our iMETHYL database (51).

In addition to demonstrating GPR15 hypomethylation at cg19859270, we found that GPR15 gene expression was significantly upregulated in PS, CS, and HTP users compared with that in NS (Fig. 4A; Supplementary Fig. S1A and S1C). These results were consistent with previous studies in which smoking-induced GPR15 hypomethylation correlated with increased gene expression (27). GPR15 expression is regulated by the direct binding of AhR to its open chromatin region, and by the interaction with the transcription factors FOXP3 (master transcription factor of certain regulatory T cells), and RORγt (retinoic acid receptor-associated orphan receptor γt) (52). Therefore, AHRR and GPR15 hypomethylation and elevated gene expression may be the result of a series of responses to immune activation by combustible tobacco smoking. Here, we found that HTP users displayed greater AHRR and GPR15 hypomethylation levels and lower AHRR expression than CS, which may reflect a tendency for subsiding inflammation in HTP users as a result of switching from combustible tobacco to HTP. However, GPR15 expression remained high in HTP and PS users, suggesting that the immune response pathway activated by combustible tobacco smoking may not recover immediately after reducing tobacco smoke exposure, but instead remains active for several years (53).

CYP4F22 (cytochrome P450, family 4, superfamily F, polypeptide 22), which encodes a fatty acid metabolizing omega-hydroxylase, displayed specifically reduced gene expression in CS. Although its function is poorly understood, CYP4F22 has been reported as a causative gene for ichthyosis (54) and was recently identified as a missing fatty acid ω-hydroxylase required for the synthesis of acylceramide, an important skin barrier lipid, suggesting that it may play an important role in skin barrier formation (55). Because no relationship between CYP4F22 and tobacco smoking has yet been reported, this study is the first to suggest its association with tobacco smoking and further studies are required to determine whether this gene is involved in the development of smoking-associated skin diseases (e.g., palmoplantar pustulosis and plaque psoriasis).

Of the 13 genes significantly downregulated in HTP users compared with those in NS (Fig. 4A), three (NXF3, DGKK, and PODN) were also significantly downregulated in CS users and have not yet been associated with tobacco smoking. Although their subcellular locations differ and the genes share no reported common pathways, they have recently been associated with hepatocellular carcinoma (56) and gastric cancer (57). Furthermore, because smoking is considered a risk for the development of those cancer types (58, 59), the use of HTP may also constitute a cancer risk factor. However, the functions of each gene have not yet been fully elucidated and further analysis is required.

Herein, we determined the smoking habits of each participant as rigorously as possible by conducting face-to-face interviews about their smoking habits and confirming their smoking exposure using cotinine concentration tests. We also standardized the conditions from sample collection to omics analysis with the help of trained technicians to collect PBMC samples using an established protocol, resulting in high-quality DNA and RNA that yielded very reliable results. Because this study was conducted in a single rural area of Japan, the results need to be carefully generalized. Nevertheless, the study also has several limitations. Because HTPs have only been on the market for a short time and the participants enrolled here only switched to HTPs for 1.7 years on average, a limitation of this study is that we have not been able to eliminate the effects of combustible tobacco from the DNAm and gene expression profiles of HTP users. Therefore, a continuous follow-up study is required to clarify the health effects of long-term HTP use. Although we defined HTP users as those who use only HTPs to clarify the molecular genetic effects of HTP use, it has been reported that in fact, more than 10% of HTP users use both HTPs and combustible cigarettes (60). Therefore, future studies should clarify the effects of concomitant use. The detailed effects of HTP use on DNAm should therefore be investigated through comprehensive DNAm analysis using microarrays or NGS. Racial differences in DNAm and transcriptome have not been discussed as much as genomic studies at present. However, because genomic differences between races may affect DNAm and transcriptome, these results need to be verified by analyzing Japanese people in other parts of Japan and other races in other countries.

To investigate the molecular genetic effects of HTP use, we analyzed 17 smoking-associated DNAm markers and whole transcriptome in a cohort study composed of Japanese participants. The DNAm status of the biomarkers and the total transcriptome profile of HTP users (switched from CS <2 years) differed from those of the NS group, albeit to a lesser extent from those of the CS. Because the molecular genetic effects of long-term HTP use remain unknown, further long-term molecular epidemiologic studies are necessary.

H. Ohmomo reports grants from Japan Agency for Medical Research and Development during the conduct of the study. S. Harada reports grants from Japan Agency for Medical Research and Development, Yamagata Prefectural Government, and the city of Tsuruoka during the conduct of the study. K. Ono reports grants from Japan Agency for Medical Research and Development during the conduct of the study. Y. Sutoh reports grants from Japan Agency for Medical Research and Development during the conduct of the study. R. Otomo reports grants from Japan Agency for Medical Research and Development during the conduct of the study as well as personal fees from ARKRAY, Inc. outside the submitted work. S. Umekage reports grants from Japan Agency for Medical Research and Development during the conduct of the study. T. Hachiya reports grants from Japan Agency for Medical Research and Development during the conduct of the study as well as personal fees from Genome Analytics Japan Inc. outside the submitted work. K. Katanoda reports grants from Japan Society for Menopause and Women's Health outside the submitted work. T. Takebayashi reports grants from Japan Agency for Medical Research and Development and Yamagata Prefectural Government and the City of Tsuruoka during the conduct of the study. A. Shimizu reports grants from Japan Agency for Medical Research and Development during the conduct of the study. No disclosures were reported by the other authors.

H. Ohmomo: Conceptualization, resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Harada: Conceptualization, resources, data curation, funding acquisition, writing–review and editing. S. Komaki: Formal analysis, validation, investigation, visualization, methodology, writing–review and editing. K. Ono: Resources, data curation, formal analysis, visualization, writing–review and editing. Y. Sutoh: Formal analysis, writing–review and editing. R. Otomo: Formal analysis, writing–review and editing. S. Umekage: Formal analysis, writing–review and editing. T. Hachiya: Formal analysis, writing–review and editing. K. Katanoda: Formal analysis, writing–review and editing. T. Takebayashi: Conceptualization, resources, data curation, supervision, funding acquisition, project administration, writing–review and editing. A. Shimizu: Conceptualization, formal analysis, supervision, methodology, project administration, writing–review and editing.

This study was supported by the Yamagata Prefectural Government and the city of Tsuruoka, and the Japan Agency for Medical Research and Development (AMED) under grant number JP19ek0210102 (T. Takebayashi). We thank Yuki Koshiba and Shizuka Saito for the collection of the healthy information, and Atsuki Kawai and Noriko Komatsu for blood collection and pre-treatment, as well as Kumi Furusawa and Takako Aoyama for their technical assistance in DNA/RNA assessment, library preparation, and sequencing.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.