Global methylation profiles in buccal cells of long-term smokers and moist snuff consumers

Abstract

Purpose: Alternations in gene methylation and other epigenetic changes regulate normal development as well as drive disease progression. The aim of this study is to investigate global methylation changes in the buccal cells of smokers and smokeless tobacco users.

Materials and methods: Generally healthy adult male subjects were recruited into smoker (SMK), moist snuff consumer (MSC) and non-tobacco consumer (NTC) cohorts (40 subjects/cohort) (ClinicalTrials.gov Identifier: NCT01923402). Global methylation profiling was performed on Illumina 450 K methylation arrays using buccal cell DNAs.

Results: The SMK cohort exhibited larger qualitative and quantitative changes relative to MSC. Approximately half of the differentially methylated 1252 gene loci were grouped as combustible tobacco-related (CTR) signatures and a third of the changes, tobacco-related (TR) signatures, were associated with smoking. Very few (41) differentially methylated gene loci were exclusively associated with moist snuff use and were designated as moist snuff-related (MSR) signature. Pathway enrichment analyses revealed that developmental and immune response pathways, among others, were impacted due to tobacco use.

Conclusions: Chronic cigarette smoking causes hyper- and hypo-methylation of genes that could contribute to smoking-related diseases. These results help place combustible and non-combustible tobacco products along a risk continuum and provide additional insights into the effects of tobacco consumption.

Keywords:

• Buccal cells
• methylation profiles
• tobacco

Introduction

Cigarette smoking is a major risk factor for several diseases, such as lung cancer, oral cancer, chronic obstructive pulmonary disease (COPD) and cardiovascular disease (CVD) (U.S. Department of Health and Human Services 2010). Chronic smoking has been reported to impact functioning of multiple organs and hence, normal physiology. Numerous biochemical changes at the molecular level that could collectively drive physiological perturbations, and ultimately progression to disease states have been identified (Buro-Auriemma et al. 2013, Besingi and Johansson 2014, Gao et al. 2015, Prasad et al. 2015, Titz et al. 2015). These molecular changes encompass epigenetic, transcriptomic, proteomic and metabolomic perturbations, and map to distinct biochemical pathways. DNA methylation is a key regulatory mechanism in normal development and disease pathways for several diseases, including cancer (Tsaprouni et al. 2014, Zhang et al. 2016). Alterations in DNA methylation status have been reported to regulate chromatin structure, and potentially influence transcription (Sneller and Gunter 1987, Kang et al. 2015, Liu et al. 2016b).

In well-documented examples, hyper- or hypo-methylation of CpG sites in the promoter regions of oncogenes and tumour suppressor genes have been shown to silence or activate those genes, respectively, consequently regulating oncogenic pathways (Ohtani-Fujita et al. 1993, Kanaya et al. 2003, Djos et al. 2012, Liu et al. 2016a). Changes in DNA methylation status have been associated with lifestyle, including chronic smoking (Lee and Pausova 2013). Healthy smokers are reported to experience changes in the methylation of genes in lung (Buro-Auriemma et al. 2013), buccal cells (Teschendorff et al. 2015) and peripheral blood mononuclear cells (PBMCs) (Gao et al. 2015). Previous work has identified differential methylation in smokers in genes involved in several biological pathways including xenobiotic metabolism (AHRR) (Tsaprouni et al. 2014).

Combustion-related toxicants are the primary drivers of biological processes that culminate in smoking-related diseases, and combustible tobacco products are recognized as the most harmful form of tobacco use. Clearly, non-use of tobacco or smoking cessation are better options for harm reduction, and use of any tobacco product presents a certain degree of risk. In fact, a relative risk between the consumption of combustible, non-combustible tobacco products and pharmaceutical nicotine therapy has been proposed. In the tobacco product risk continuum, cigarette smoking is regarded as the most harmful form of tobacco consumption followed by non-combustible (smokeless) tobacco products, which in turn are more hazardous than the pharmaceutical nicotine products (Hatsukami et al. 2007, Zeller and Hatsukami 2009).

Smokeless tobacco (ST) consists of diverse non-combustible oral tobacco products. STs are consumed globally but significantly vary in the type, form and chemical composition. For example, moist snuff is the most widely consumed form of ST in the United States (US) chemically, moist snuff contains nicotine and a number of toxicants; for example tobacco-specific nitrosamines (TSNAs), cadmium and benzo(a)pyrene (Stepanov et al. 2008, Borgerding et al. 2012, Song et al. 2016). Generally, the epidemiology relating to Swedish and US STs suggest that ST consumers experience lower risk for lung cancer, COPD, and oral cancer compared to smokers (Luo et al. 2007, Rodu 2007, Colilla 2010, Lee 2013).

Considering that the tobacco use-related diseases develop over a long period of product usage, short-term clinical measurements including biomarkers of exposure and effect offer insights into the effects of product usage prior to manifestation of clinical symptoms. In this context, the early biological effects of chronic cigarette smoking are fairly well understood (U.S. Department of Health and Human Services 2010), however, such information is very limited regarding the use of smokeless tobacco products, including moist snuff.

Several biomarker studies, including those from R.J. Reynolds Tobacco (RJRT), have shown that moist snuff consumers (MSC), although exposed to nicotine and TSNAs, exhibit significantly lower levels of combustion-related biomarkers than smokers, with biomarker levels comparable to those found in the non-tobacco consumers (NTC) (Campbell et al. 2015, Prasad et al. 2016). Further, several biomarkers that are indicative of the effect of tobacco consumption on the consumers (known as biomarkers of effect) across different biological pathways were not significantly different from NTC. Smokers (SMKs), however, exhibited significant differences in these biomarkers, relative to the consumers of STs and NTC (Campbell et al. 2015, Prasad et al. 2016).

In efforts to identify novel candidate biomarkers of tobacco effect and gain insights into the molecular changes resulting from long-term consumption of different tobacco product categories, we utilized ‘omic’ technologies. Recently, we showed that global metabolomic profiles from NTC and MSC are similar, but the smokers’ profiles are distinct (Prasad et al. 2015). Since cigarette smoking is known to be associated with changes in DNA methylation status of select loci and globally (Shenker et al. 2013, Sun et al. 2013, Tsaprouni et al. 2014, Teschendorff et al. 2015, Gao et al. 2016), here we sought to characterize global methylation changes in the buccal cells of chronic smokers and MSC.

Clinical significance

• Cigarette smoking is a major risk factor for several diseases and altered DNA methylation is believed to contribute to the pathophysiology of smoking-related diseases.

• This work for the first time comparatively assesses buccal cell DNA methylation profiles in smokers and users of non-combustible tobacco.

• Cigarette smoking appears to be a potent modifier of altered methylation compared to smokeless tobacco, which is consistent with the existing epidemiology.

• This work provides an understanding of the potential physiological impact of different tobacco products.

• The differentially methylated loci, including AHRR, serve as potential biomarkers that distinguish the effects of tobacco product use.

Methods

Study population and clinical conduct

This cross-sectional study was conducted between June 2010 and January 2011 in accordance with applicable sections of the US Code of Federal Regulations (21 CFR 50, 54, 56) and the International Conference on Harmonization Guidance on Good Clinical Practice, which is consistent with Declaration of Helsinki). All pertinent study documents (i.e. protocol, amendments, informed consent forms and associated subject materials) were reviewed and approved by Independent IRB, Inc. (currently Schulman IRB, Cincinnati, OH). The study was registered with the Clinical Trials Registry (Clinical Trials.gov) with the identifier NCT01923402 as described previously (Prasad et al. 2016). Adult male subjects, aged 35–60 years, were enrolled into three cohorts of 40 subjects each: exclusive SMKs, MSC and NTC. The sample size for this biomarker study was based on the estimates for differentiating metabolomic profiles and those estimates were derived from historical experience with metabolomic profiling.

Briefly, the inclusion criteria were: (1) exclusive cigarette SMKs of any brand ≥6 mg ‘tar’ (measured by the Cambridge Filter Pad method), who self-reported smoking at least 10 cigarettes per day for at least 3 years and had an expired carbon monoxide (ECO) level of 10–100 ppm; (2) exclusive MSC of any brand, who self-reported using ≥2 cans of moist snuff/week for at least 3 years and had an ECO of 0–5 ppm; (3) NTC, who self-reported not using any tobacco or nicotine-containing products for at least 5 years and had an ECO of 0–5 ppm. The exclusive smokers smoked for 25.1 ± 9.6 years and 21.5 ± 5.3 cigarettes per day. The MSC have reported 20.6 ± 8.5 years of product use and 6.3 ± 3.5 cans of moist snuff per week (Prasad et al. 2016).

The Consort Flow Diagram below illustrates the subject flow in the study. Buccal cells were collected by swishing method following a 2-h fasting from food and tobacco. Subjected rinsed their mouth with Scope mouthwash followed by a water rinse. Buccal cells were collected in water by vigorously swishing. The cell pellet was washed in phosphate buffered saline and used for DNA extraction.

Global methylation profiling

DNA extraction and global methylation profiling of 485,577 CpG sites were performed by Expression Analysis, Inc., on Illumina Infinium HumanMethylation450 BeadChip arrays. Statistical analysis was performed by Covance Genomics Laboratories, Seattle, WA. All data were normalized using an Illumina recommended normalization method to reduce the effects of experimental variation, which involved background estimation using negative control probes and normalization of intensities to the housekeeping genes with no CpG sites. Beta values (β), which were used to estimate the methylation level of the CpG locus using the ratio of intensities between methylated and un-methylated alleles, were exported from Illumina Genome Studio software and all analyses were performed using the SAS statistical software package (SAS Inc., Cary, NC). A stringent filter based on the standard deviation (SD), which is sensitive to extreme observations and interquartile range, which is more resistant to outliers, was used to exclude loci with outlier beta values and low biological variability. For normally distributed beta values, we analysed the data by modelling each individual methylation locus by linear regression with age and race as independent variables (Model 1).

Model 1: Methylation (β) ∼ Product + Race + Age.

All genome-wide comparisons were corrected for multiple comparisons using the method of Benjamini and Hochberg (1995). Differentially, methylated loci associated with tobacco products were identified as those that resulted below the 5% false discovery rate (FDR) thresholds. We identified 1252 loci that were differentially methylated between combustible, non-combustible and non-tobacco users. In order to characterize these selected 1252 loci, we empirically classified them into five main signature classes using the pair-wise contrast p values (< 0.05) between the three cohorts. The classes and rules for the classification are shown in Figure 1(A).

Figure 1. Classification of differentially methylated gene loci. A. Rules for classifying loci with FDR p value <0.05 (contrast p values). B. Major categories of differentially methylated gene loci and number of loci per category.

Methylation clustering

Raw methylation values were exported from Illumina Genome Studio and imported into R. Methylation values were colour-adjusted and quantile normalized using the R Bioconductor lumi package. The top 50 loci heat map was created with the R gplots package (R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org).

Determination of transcriptional changes based on methylation status

A promoter region was assigned to the target ID if the Illumina HumanMethylation BeadChip annotation or UCSC RefGene Group referenced a ‘TSS’ (transcriptional start site) or ‘5'UTR.’ A gene region was assigned to the target ID if the annotation for UCSC RefGene Group referenced ‘Body,’ ‘Exon’ or ‘3'UTR’ and did not include ‘TSS’ or ‘5'UTR.’ Entrez gene IDs that consisted of loci both hyper- and hypo-methylated were assigned to a promoter or gene region first based on the larger number of loci and, if equal, by the most significant loci.

Network analysis

In addition to evaluating the functional enrichment of three major group gene signatures, we also performed network analysis using two different algorithms in Thomson Reuters MetaCore (http://thomsonreuters.com/metacore) and evaluated the functional enrichment of the resulting network(s). 1) Dijkstra's shortest path algorithm was used to calculate the smallest possible number of directed one-step interactions between pairs of initial objects in each direction. 2) An analyse network algorithm generated 50-node sub-networks enriched with seed nodes. We prioritized the sub-networks according to their p-score and evaluated the top four networks that met an empirically determined criteria (i.e. contained 10 or more seed nodes).

Enrichment analysis

We performed functional enrichment analysis for Gene Ontology Biological Process (GO:BP) using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) (Huang da et al. 2009a, 2009b). Enrichment analysis often identifies multiple pathways/processes down the same branch of the Gene Ontology graph; pathways/processes further down the graph typically provide more information but are less significant and have fewer associated genes than higher pathways/processes. To ensure judicious, systematic selection of pathways/processes and to minimize over-representation, the top 120 terms from each enrichment were visualized and manually reviewed to select those that best balance informative content and gene enrichment. As a means to partition pathways/processes for review and selection, prior to visualization long lists of Gene Ontology terms were summarized by finding a representative subset of terms using a simple clustering algorithm that relies on semantic similarity measures (Supek et al. 2011). All six enrichments were pooled and integrated, and the pathways/processes were then categorized and prioritized according to the percentage of signature genes in each category.

Results

In a single-site cross-sectional study in NC, USA, healthy subject populations were enrolled into three subject cohorts: exclusive MSC, exclusive cigarette SMK and never-tobacco consumers (NTC) (Supplementary Table S1a). Each cohort consisted of 40 generally healthy males between the ages of 35–60 (inclusive). Descriptive statistics show age distributions were similar for all three groups (Supplementary Table S1b). A majority of subjects (77%) were of Caucasian background followed by African American (Prasad et al. 2016). The average age was 45.0 (MSC), 46.9 (SMK) and 47.2 (NTC), with years of tobacco product use 20.6 ± 8.5 for MSC and 25.1 ± 9.6 for SMK.

Buccal cell DNAs were prepared and subjected to methylation profiling using the Illumina 450 K chip array. For each sample, at least 99.4% of total CpG loci were detected. Unsupervised analysis of the beta-values, which are used to estimate the methylation level of the CpG locus using the ratio of intensities between methylated and un-methylated alleles, showed that there are two major groups of different methylation status. Clustering analysis of the 3000 loci with most variable beta values (i.e. highest SDs) showed two major clusters that appeared to correspond to the two major race groups in the samples (data not shown).

DNA methylation associated with tobacco use

DNA methylation status was measured in 120 buccal cell samples (MSC, n = 40; SMK, n = 40; NTC, n = 40) collected from study subjects on day 1 following a 2-h fasting from tobacco and food. We identified 1252 CpG sites with differential methylation levels related to tobacco use after stringent filtering and correction for multiple testing (FDR < 0.05). To characterize these loci, we empirically classified them into three major groups using the pair-wise contrast p values (<0.05) (Figure 1(A)). The combustible tobacco-related (CTR) signature consists of loci whose methylation changes are only observed in the SMK group. Similarly, the moist snuff-related (MSR) signature consists of loci whose methylation changes are only detected in the MSC group. The tobacco-related (TR) signature is characterized by common methylation aberrations observed in both the SMK and MSC cohorts. We analysed each of these tobacco use methylation signatures by performing enrichment analysis below.

The largest signature – the CTR signature – consists of 690 loci (55.1% of the 1252 loci identified) and is associated with the SMK group (Figure 1(B)). It is characterized by decreased levels of methylation for the SMK cohort compared to both the SMK and NTC cohorts (Figure 2(A)) The TR signature consists of the next highest number of loci (391; 31.2% of the 1252 loci identified) and is characterized by decreased levels of methylation in the NTC cohort, with the SMK cohort having the highest levels, followed by the SMK cohort with intermediate levels and trailed by the NTC cohort with the lowest levels (Figures 1(B) and 2(B)). Thus, the TR signature is composed of methylation changes in both SMK and MSC cohorts. CTR and TR signatures with ‘hits’ in more than one locus are provided in Table 1; each loci identified in the MSR signature was hit once.

Figure 2. Methylation group example signatures based on beta values. A. Combustible tobacco-related (CTR) signature for cg05575921 (AHRR). B. Tobacco-related (TR) signature for cg04018738 (VARS). C. Moist snuff-related (MSR) signature for cg14895183 (AGAP1).

Table 1. List of CTR and TR top hits with more than one locus altered for each signature.

CTR signature top hits
Gene symbol	Markers	Regulation	Gene annotation
AHRR	7	Hypo-methylated	Aryl-hydrocarbon receptor repressor
GNG12	5	Hypo-methylated	Guanine Nucleotide Binding Protein (G Protein), Gamma 12
TCN2	5	Hyper-methylated	Transcobalamin II
LOC404266	5	Hypo-methylated
DLEU7	4	Hypo-methylated	Deleted in lymphocytic leukaemia, 7
GAPDHS	4	Hypo-methylated	Glyceraldehyde-3-phosphate dehydrogenase, spermatogenic
GFI1	4	Hypo-methylated	Growth factor independent 1 transcription repressor
PRR15	4	Hypo-methylated	Proline rich 15
CACNA1D	3	Hyper-methylated	Calcium channel, voltage-dependent, L type, alpha 1D subunit
NR2F2	3	Hyper-methylated	Nuclear receptor subfamily 2, group F, member 2
C10orf41	3	Hypo-methylated
CYP1A1	3	Hypo-methylated	Cytochrome P450, family 1, subfamily A, polypeptide 1
HS3ST3B1	3	Hypo-methylated	Heparan sulphate (glucosamine) 3-O-sulfotransferase 3B1
CREB5	2	Hyper-methylated	cAMP responsive element binding protein 5
GFAP	2	Hyper-methylated	Glial fibrillary acidic protein
ID3	2	Hyper-methylated	Inhibitor of DNA binding 3, dominant negative helix-loop-helix protein
IGF1	2	Hyper-methylated	Insulin-like growth factor 1 (somatomedin C)
KCNK2	2	Hyper-methylated	Potassium channel, subfamily K, member 2
PKNOX2	2	Hyper-methylated	PBX/knotted 1 homeobox 2
PRKAR1B	2	Hyper-methylated	Protein kinase, cAMP-dependent, regulatory, type I, beta
SPSB4	2	Hyper-methylated	splA/ryanodine receptor domain and SOCS box containing 4
TSLP	2	Hyper-methylated	Thymic stromal lymphopoietin
TR signature top hits
SEPT9	4	Hyper-methylated	Septin 9
RASSF1	4	Hyper-methylated	Ras Association (RalGDS/AF-6) Domain Family Member 1
BST1	3	Hyper-methylated	Bone marrow stromal cell antigen 1
PITX2	3	Hyper-methylated	Paired-like homeodomain 2
RAD52	3	Hyper-methylated	RAD52 homolog
ARHGEF6	2	Hyper-methylated	Rac/Cdc42 guanine nucleotide exchange factor (GEF) 6
CARS	2	Hyper-methylated	Cysteinyl-tRNA synthetase
GPRC5C	2	Hyper-methylated	G protein-coupled receptor, family C, group 5, member C
NKX3-2	2	hyper-methylated	NK3 homeobox 2
PAX1	2	Hyper-methylated	Paired box 1
PDXK	2	Hyper-methylated	Pyridoxal (pyridoxine, vitamin B6) kinase
PER1	2	Hyper-methylated	Period circadian clock 1

Very few loci were detected in the MSR Signature, consisting of 41 loci (3.3% of the 1252 loci identified) associated with MSC. It is characterized by a decrease in levels of methylation for the MSC cohort compared to both the SMK and NTC cohorts (Figures 1(B) and 2(C)). Hierarchical clustering of the top 50 significant differential methylation gene loci co-cluster MSC and NTC cohorts, and a majority of SMK appear to be distinct from the MSC and NTC (Figure 3), consistent with studies that have evaluated combustion-related biomarkers (Campbell et al. 2015, Prasad et al. 2016).

Figure 3. Six of seven AHRR loci are significantly hypo-methylated in the SMK cohort.

Predicted gene expression from differential methylation

In order to provide a more complete functional interpretation of each signature, we calculated the predicted transcriptional impact of each differentially methylated site. DNA methylation in gene promoters is associated with down-regulation of gene expression while methylation in gene bodies is associated with enhanced DNA transcription, i.e. up-regulation (Lee and Pausova 2013). Thus, we predicted 256 up- and 201 down-regulated genes in the CTR signature, 130 up- and 102 down-regulated genes in the TR signature and 8 up- and 27 down-regulated genes in the MSR signature (Supplementary Tables S2–S4).

Candidate biomarker identification

For this study, three criteria were used to identify candidate biomarkers in smokers and/or MSC: differential methylation, biological relevancy and status. Specifically, each biomarker candidate met the following: (1) the gene was hypo- or hyper-methylated; (2) the gene or gene product was associated with a specific biological process or pathway; and (3) the candidate was currently considered a biochemical, genomic or proteomic biomarker in some context as catalogued in the Thomson Reuters Integrity database.

Combustible tobacco-related (CTR) methylation signature

The CTR methylation signature consisted of 510 statistically significant differentially methylated target IDs that mapped to 147 hyper-methylated and 305 hypo-methylated genes. The most hypo-methylated site was AHRR, comprising 5 of the top 10 most significant genes (Supplementary Table S5). In the SMK cohort, AHRR was hypo-methylated at 6 of 7 loci (Figure 4)

Figure 4. Top 50 loci with the most significant adjusted p values cluster into two groups. The majority of NTC (green; 80%) and MSC (blue; 63%) samples co-cluster to the right side of the heat map (circled in black), while the majority of SMK (red; 68%) samples cluster to the left. Coloured heat map regions indicate methylation changes: hypo-methylation (green) and hyper-methylation (red).

Enrichment analysis of the 452 CTR genes identified nine significantly associated (p value = 5.89 × 10⁻⁶) with response to the retinoid vitamin A: AQP3, BMP4, CD44, CYP1A1, HSD17B2, IGFBP7, PDGFA, RARA and RARG (Table 2). In rats, cigarette smoke-induced depletion of vitamin A is associated with the development of emphysema (Li et al. 2003). CYP1A1, one of the main cytochrome P450 enzymes, is one of four genes predicted to be up-regulated in the response to vitamin A association. CYP1A1 polymorphisms are associated with risk of oral, oesophageal and lung cancers. Ten genes were significantly associated (p value = 1.12 × 10⁻⁴) with inner ear development: CDH23, CHD7, CUX1, DLL1, GFI1, PDGFA, PDGFRA, TCAP, TGFB3 and WNT3A; seven genes were significantly associated (p value = 3.62 × 10⁻²) with respiratory tract development: BMP4, CUX1, FGF18, HS6ST1, PDGFA, PDGFRA and TGFB3.

Table 2. CTR network analysis summary.

Enriched category/pathway/process	Analysis enrichment	Signature gene enrichment (%)	Hyper-methylated genes	Hypo-methylated genes	Derived up-regulated signature gene expression	Derived down-regulated signature gene expression
Response to drug	Shortest paths and 2 of 4 networks	24	COL18A1 and NGF	CYP1A1, LCK, PDGFA, RELA, SEMA3C and SOCS3	COL18A1, CYP1A1 and SEMA3C	LCK, NGF, PDGFA, RELA and SOCS3
Response to vitamin A	Signature	100	None	AQP3, BMP4, CD44, CYP1A1, HSD17B2, IGFBP7, PDGFA, RARA and RARG	BMP4, CYP1A1, HSD17B2 and RARA	AQP3, CD44, IGFBP7, PDGFA and RARG
Mesoderm development	Signature, Shortest paths and 1 of 4 networks	35	SECTM1 and WNT3A	BMP4, EXT1, PPP2CA and TXNRD1	BMP4, SECTM1 and WNT3A	EXT1, PPP2CA and TXNRD1
Respiratory system development
Respiratory tube development	Signature, Shortest paths and 1 of 5 networks	37	CUX1	BMP4, FGF18, HS6ST1, PDGFA, PDGFRA and TGFB3	BMP4, CUX1, HS6ST1 and TGFB3	FGF18, PDGFA and PDGFRA
Lung development	Shortest paths and 1 of 4 networks	33	CUX1	BMP4, FGF18, PDGFA, PDGFRA and TGFB3	BMP4, CUX1 and TGFB3	FGF18, PDGFA and PDGFRA
Respiratory system development	Shortest paths and 1 of 4 networks	37	CHD7 and CUX1	BMP4, FGF18, PDGFA, PDGFRA and TGFB3	BMP4, CHD7, CUX1 and TGFB3	FGF18, PDGFA and PDGFRA
Circulatory system development
Heart development	Signature, Shortest paths and 1 of 4 networks	36	BCOR, CHD7, ID3, KCNJ8, NFATC4, NRP2, SMARCD3 and WNT3A	ADAM19, BMP4, MED1, RXRB, SEMA3C, TCAP and ZMIZ1	BMP4, CHD7, ID3, KCNJ8, MED1, NFATC4, NRP2, RXRB, SEMA3C, SMARCD3, TCAP, WNT3A and ZMIZ1	ADAM19 and BCOR
Vasculature development	Shortest paths	50	CHD7, COL18A1, ELK3, NR2F2, NRP2, SH2D2A and TCF7L2	BMP4, CD44, FGF18, PDGFA, PTK2 and SEMA3C	BMP4, CHD7, COL18A1, ELK3, NRP2, SEMA3C and TCF7L2	CD44, FGF18, NR2F2, PDGFA, PTK2 and SH2D2A
Blood vessel development	Signature	100	CHD7, COL18A1, ELK3, NR2F2, NRP2, SH2D2A and TCF7L2	BMP4, CD44, FGF18, HS6ST1, PDGFA, PTK2, SEMA3C and ZMIZ1	BMP4, CHD7, COL18A1, ELK3, HS6ST1, NRP2, SEMA3C, TCF7L2 and ZMIZ1	CD44, FGF18, NR2F2, PDGFA, PTK2 and SH2D2A
Oral/digestive system development
Palate development	2 of 5 networks	45	BCOR and CHD7	MEN1, PDGFRA and TGFB3	CHD7 and TGFB3	BCOR, MEN1 and PDGFRA
Inner ear development	Signature and 1 of 5 networks	71	CHD7, CUX1 and WNT3A	CDH23, DLL1, GFI1, PDGFA, PDGFRA, TCAP and TGFB3	CHD7, CUX1, TCAP, TGFB3 and WNT3A	CDH23, DLL1, GFI1, PDGFA and PDGFRA
Immune response
Cytokine-mediated signalling	2 of 5 networks	20	CLCF1 and JAK3	RELA	None	CLCF1, JAK3 and RELA
Positive regulation of T cell activation	2 of 5 networks	20	IL2	LCK and RARA	IL2 and RARA	LCK
Response to hypoxia	2 of 5 networks	29	None	CREBBP, CYP1A1, PDGFA, SOCS3 and TGFB3	CYP1A1 and TGFB3	CREBBP, PDGFA and SOCS3
Regulation of transforming growth factor beta receptor signalling	Signature and 1 of 5 networks	56	None	FURIN, MEN1, PRDM16, SKI and TGFB3	TGFB3	FURIN, MEN1, PRDM16 and SKI
Regulation of cell death	Shortest paths and 1 of 5 networks	33	CLCF1, COL18A1, ID3, IGF1, IL2, IL3, NGF, NTRK1, SOCS2, TCF7L2 and TRAF3	AHRR, BCL6, BIK, BMP4, CD44, CLU, FURIN, LCK, MCF2, MEN1, NQO1, NUAK2, PPP2CA, PRKCA, RARG, RELA, SOCS3 and TGFB3	BCL6, BIK, BMP4, CLU, COL18A1, ID3, IL2, MCF2, SOCS2, TCF7L2, TGFB3 and TRAF3	AHRR, CD44, CLCF1, FURIN, IGF1, IL3, LCK, MEN1, NGF, NQO1, NTRK1, NUAK2, PPP2CA, PRKCA, RARG, RELA and SOCS3
Negative regulation of cell cycle	Signature	100	NGF and DGKZ	BCL6, BMP4, MEN1 and MYO16	BCL6, BMP4 and DGKZ	MEN1, MYO16 and NGF

The number of analyses in which each enriched category/pathway/process was identified is indicated. Percent signature gene enrichment (i.e. detected methylation changes) is shown, along with a breakdown of hyper- or hypo-methylation and derived up- or down-regulated gene expression.

Moist snuff-related (MSR) methylation signature

The MSR methylation signature consisted of 32 statistically significant differentially methylated target IDs that mapped to 10 hyper-methylated and 25 hypo-methylated genes. Each gene was detected once. The gene most significantly hypo-methylated and predicted to be down-regulated was ST6GALNAC5 (Supplementary Table S6). ST6GALNAC5 belongs to a family of sialyltransferases that modifies proteins and ceramides on the cell surface to alter cell–cell or cell–extracellular matrix interactions. A recent study suggested that ST6GALNAC5 mutations can cause coronary artery disease (InanlooRahatloo et al. 2014). Enrichment analysis of the 35 MSR genes identified four significantly associated (p value = 0.0097) with microtubules: TUBGCP2, CLASP2, KIF15 and FSD1. All four genes are predicted to be down-regulated (Table 3). A majority (21 of 35) of MSR genes are known to encode proteins for which at least two isoforms exist due to distinct pre-mRNA splicing events (Swiss-Prot Protein Information Resource Keywords, p value = 5.87 × 10⁻³). These include three of the four genes associated with microtubules: CLASP2, KIF15 and TUBGCP2. Between the CTR and the MSR signatures, only one common gene MCF2 is found, indicating that cigarette smoking and moist snuff usage induce markedly different epigenetic modulation in the buccal cells of the consumers, with cigarette smoking being the predominant driver of perturbations in the gene methylation.

Table 3. MSR network analysis summary table.

Enriched category/pathway/process	Analysis enrichment	Signature gene enrichment (%)	Hyper-methylated genes	Hypo-methylated genes	Derived up-regulated signature gene expression	Derived down-regulated signature gene expression
Microtubules	Signature	100	KIF15 and TUBGCP2	CLASP2 and FSD1	None	CLASP2, FSD1, KIF15 and TUBGCP2\
Alternative splicing	Signature	100	C6orf48, DHX36, KIF15, SNORD52, TUBGCP2, VAX1, ZNF236 and ZNF511	AGAP1, CLASP2, DIP2A, HUWE1, MCF2, MINPP1, PHF21A, RIMS1, SBF2, SKIL, SLC9A3R2, TMEM126B, ZDHHC16 and ZFYVE28	MCF2, RIMS1, SKIL, TMEM126B and ZDHHC16	AGAP1, C6orf48, CLASP2, DHX36, DIP2A, HUWE1, KIF15, MINPP1, PHF21A, SBF2, SLC9A3R2, SNORD52, TUBGCP2, VAX1, ZFYVE28, ZNF236 and ZNF511
Wnt receptor signalling pathway	Shortest paths and 3 of 4 networks	4	None	NLK	NLK	None
Regulation of cell development	Shortest paths	4	VAX1	None	None	VAX1
Regulation of programmed cell death	Shortest paths	3	None	LAMTOR5 and MCF2	MCF2	LAMTOR5
Mitosis	Shortest paths and 1 of 4 networks	13	KIF15	CLASP2 and FSD1	None	CLASP2, FSD1 and KIF15

The number of analyses in which each enriched category/pathway/process was identified is indicated. Percent signature gene enrichment (i.e. detected methylation changes) is shown, along with a breakdown of hyper- or hypo-methylation and derived up- or down-regulated gene expression.

Tobacco-related (TR) methylation signature

The TR methylation signature consisted of 255 statistically significant differentially methylated target IDs that mapped to 214 hyper-methylated and 18 hypo-methylated genes. The gene most significantly hyper-methylated and predicted to be up-regulated was VARS. Two genes that are known tumour suppressors had the largest number of hyper-methylated loci: RASSF1 and SEPT9 (Table 1; Supplementary Table S7). Both genes were methylated at four different loci in the promoter region and predicted to be down-regulated. Enrichment analysis of the 232 TR genes identified 13 significantly associated (p value = 7.32 × 10⁻⁵) with hematopoietic or lymphoid organ development: BCL11A, BCL11B, BCL2L11, FOXP1, HDAC4, HOXA9, ITPKB, JARID2, MAEA, NKX3-2, PAX1, TGFBR2 and WNT3A (Table 4).

Table 4. TR network analysis summary table.

Enriched category/pathway/process		Analysis enrichment	Signature gene enrichment (%)	Hyper-methylated genes	Hypo-methylated genes	Derived up-regulated signature gene expression	Derived down-regulated signature gene expression
Hemopoietic or lymphoid organ development		Signature	100	BCL11A, BCL11B, BCL2L11, FOXP1, HDAC4, HOXA9, ITPKB, MAEA, NKX3-2, PAX1, TGFBR2, WNT3A	JARID2	BCL11A, BCL11B, BCL2L11, HDAC4, ITPKB, MAEA, NKX3-2, PAX1, TGFBR2, WNT3A	FOXP1, HOXA9, JARID2
Skeletal system development		Signature	100	HDAC4, HOXA9, HOXC4, HOXC5, MMP2, NKX3-2, NODAL, PAX1, SIX2, TGFBR2	None	HDAC4, NKX3-2, NODAL, PAX1, SIX2, TGFBR2	HOXA9, HOXC4, HOXC5, MMP2
Ectoderm development		Shortest paths, 3 of 5 networks	23	COL5A3, PPARD, ZBTB17	None	COL5A3	PPARD, ZBTB17
Respiratory system development
Lung development		Signature, Shortest paths, 2 of 5 networks	38	CRISPLD2, CUX1, FOXP1, NODAL, TGFBR2	None	CRISPLD2, CUX1, NODAL, TGFBR2	FOXP1
Respiratory system  development		Signature, 2 of 5 networks	38	CRISPLD2, CUX1, FOXP1, NODAL, TGFBR2	None	CRISPLD2, CUX1, NODAL, TGFBR2	FOXP1
Circulatory system development
Heart development		Signature, Shortest paths, 2 of 5 networks	50	ALPK3, FOXP1, KIF3A, NFATC4, NODAL, TGFBR2, WNT3A	None	ALPK3, NFATC4, NODAL, TGFBR2, WNT3A	FOXP1, KIF3A
Vasculature  development		Shortest paths	50	MMP2, NODAL, TGFBR2	None	NODAL, TGFBR2	MMP2
Oral/digestive system development
Inner ear  morphogenesis		Shortest paths, 1 of 5 networks	25	WNT3A	None	WNT3A	None
Ear development		Signature, 1 of 5 networks	67	BCL2L11, CUX1, NKX3-2, SIX2, WNT3A, ZIC1	None	BCL2L11, CUX1, NKX3-2, SIX2, WNT3A	ZIC1
Immune response
T cell differentiation		Signature	100	BCL11A, BCL11B, ITPKB, PAX1	None	BCL11A, BCL11B, ITPKB, PAX1	None
Response to  cytokine stimulus		Shortest paths, 1 of 5 networks	20	HDAC4	None	HDAC4	None
Response to wounding		Signature, 1 of 5 networks	68	ADORA3, C3AR1, ELK3, F2RL1, HDAC4, NFATC4, PPARD, RTN4RL1, TGFBR2, TLR1, TNFAIP6, TOLLIP, VWF	None	ADORA3, F2RL1, HDAC4, NFATC4, RTN4RL1, TGFBR2, TOLLIP, VWF	C3AR1, ELK3, PPARD, TLR1, TNFAIP6
Transforming growth  factor beta receptor signaling pathway		2 of 5 networks	17	TGFBR2	None	TGFBR2	None
Cell death
Negative regulation   of apoptosis		Shortest paths, 3 of 5 networks	9	IL3, NKX3-2	None	NKX3-2	IL3
Induction of apoptosis  by extracellular signals		Signature	100	ARHGEF6, BCL2L11, DAPK1, MCF2L, OBSCN	None	ARHGEF6, BCL2L11, DAPK1, MCF2L, OBSCN	None
Cell cycle arrest		Shortest paths, 2 of 5 networks	30	CUL1, DST, RASSF1	None	None	CUL1, DST, RASSF1

The number of analyses in which each enriched category/pathway/process was identified is indicated. Percent signature gene enrichment (i.e. detected methylation changes) is shown, along with a breakdown of hyper- or hypo-methylation and derived up- or down-regulated gene expression.

Several additional TR associations are linked to these processes: skeletal system development (p value = 0.0152), response to wounding (p value = 0.0267), T cell differentiation (p value = 0.0436) and heart development (p value = 0.0454).

Integrated systems biology analysis

Pathways involved in the response to moist snuff, cigarette or general tobacco use were determined by performing a series of network analyses on the CTR, TR and MSR methylation signature gene lists using the Thomson Reuters systems biology functional analysis platform MetaCore. To interpret epigenetic changes in the context of their impact on pathways, we expanded each list by identifying secondary genes/proteins that are known to interact with the hyper- or hypo-methylated signature genes. Smoking-associated DNA methylation changes in buccal cells have been shown to correlate with DNA methylation changes in epithelial cancers (Teschendorff et al. 2015). Based on MetaCore filtering, CTR and TR network genes were restricted to those expressed in epithelial cells; MSR network genes were restricted to those expressed in oral tissues.

Network models were built with MetaCore using two different algorithms: ‘Shortest Paths’ and ‘Analyze Network’. The Shortest Paths algorithm generates a single network; for the Analyse Network algorithm, those networks that contained more than 10 signature genes/proteins as well as the network with the most significant p value that included pathways were selected for analysis. Once created, each network was analysed for pathway/gene set enrichment. All enrichments for a given signature were then integrated and pathways/processes categorized. This integrative systems biology approach identified a total of 413 genes (CTR), 196 genes (TR) and 362 genes (MSR) significantly associated with at least one function or pathway. Below, we discuss analysis limited to pathways relevant to smoking-induced diseases that were identified using the ‘Shortest Paths’ algorithm or in multiple ‘Analyze Network’ analyses and which contained at least one differentially methylated gene (Tables 2–4).

Developmental pathways

The CTR signature and two of five CTR network analyses are significantly enriched for genes associated with mesoderm development (17 genes; 35% (6) are hyper- or hypo-methylated). The mesoderm is one of the three primary germ layers in the very early embryo that forms a variety of tissues, including mesenchyme (connective tissue), mesothelium, serous membranes, muscle (smooth and striated), bone, cartilage, adipose tissue, circulatory system, lymphatic system, dermis, genitourinary system and notochord. The related to mesoderm development consists of several that encode nuclear/cytoplasmic (5/6) and membrane/extracellular (2/4) proteins. Candidate biomarkers include BMP4, PPP2CA, TXNRD1 and WNT3A. Consistent with these data, others have shown that in vivo exposure to cigarette smoke influenced gene expression in umbilical cord tissue, predominantly organ systems derived from mesoderm (haematological, musculoskeletal, reproductive and cardiovascular) (Liszewski et al. 2012).

In contrast, the TR signature was enriched with genes associated with ectoderm development, identified in four of six TR network analyses (13 genes; 23% (3) are hyper-methylated). The ectoderm is one of the three primary germ layers in the very early embryo that differentiates to form, among other tissues, the nervous system. The ectoderm development signature contains an equal number of genes that encode nuclear (4), cytoplasmic (4) and extracellular (4) proteins; 1 gene encodes a membrane protein. Candidate biomarkers include COL5A3, PPARD and ZBTB17; all three genes are hyper-methylated.

Consistent with the known mechanisms of action and effects of combustible tobacco in the body, pathways and/or processes related to respiratory and circulatory development were also significantly impacted. Three pathways or processes identified in the CTR signature and two of five CTR network analyses are involved with respiratory system development (20 total genes, 40% (8) are hyper- or hypo-methylated). The CTR respiratory system development signature contains slightly more genes that encode membrane/extracellular (5/7) than nuclear/cytoplasmic (4/4) proteins. Features identified include the indicated candidate biomarkers hyper- or hypo-methylated: respiratory tube development (BMP4, CUX1, FGF18, PDGFA, PDGFRA and TGFB3), lung development (BMP4, CUX1, FGF18, PDGFA, PDGFRA and TGFB3) and respiratory system development (BMP4, CHD7, CUX1, FGF18, PDGFA, PDGFRA and TGFB3). All candidate biomarkers for CTR respiratory system development – BMP4, FGF18, PDGFA, PDGFRA and TGFB3 – were hypo-methylated except CHD7 and CUX1. In the TR signature and three of six TR network analyses, three pathways or processes consisting of 12 total genes (42% (5) are hyper- or hypo-methylated) associated with respiratory system development were identified. Similar to the CTR respiratory signature, the TR respiratory signature contains slightly more genes that encode membrane/extracellular (2/5) than nuclear/cytoplasmic (4/1) proteins. Features identified include the indicated methylated candidate biomarkers: lung development (CUX1, FOXP1, NODAL and TGFBR2) and respiratory system development (CUX1, FOXP1, NODAL and TGFBR2). In contrast to the CTR respiratory system development signature, all four candidate biomarkers for TR respiratory system development – CUX1, FOXP1, NODAL and TGFBR2 – were methylated.

Eight pathways or processes identified in the CTR signature and three of five CTR network analyses are involved with circulatory system development (70 total genes, 35.7% (25) are hyper- or hypo-methylated). The CTR circulatory system development signature contains a high number of genes that encode membrane proteins (27) relative to cytoplasmic (17), nuclear (19) or extracellular (7) proteins. Features identified in the category include the indicated candidate biomarkers hyper- or hypo-methylated: heart development (ADAM19, BCOR, BMP4, CHD7, ID3, KCNJ8, MED1, NFATC4, NRP2, RXRB, SEMA3C, SEMA3C, SMARCD3 and WNT3A), vasculature development (BMP4, CD44, CHD7, COL18A1, ELK3, FGF18, NR2F2, PDGFA, PTK2, SH2D2A and TCF7L2) and blood vessel development (BMP4, CD44, CHD7, COL18A1, ELK3, FGF18, NR2F2, NRP2, PDGFA, PTK2, SEMA3C, SH2D2A and TCF7L2). Tobacco has known effects on smooth muscle contraction, vasoconstriction and regulation of heart rate, all of which were identified by pathway analysis (Chiba et al. 2005, Xu et al. 2010).

In the TR signature and four of six TR network analyses, five pathways or processes associated with circulatory system development were identified (21 total genes; 38.1% (8) are hyper- or hypo-methylated). The TR circulatory system development signature contains almost equally divided numbers of genes that encode nuclear/cytoplasmic (5/6) and membrane/extracellular (3/7) proteins. Features identified in the category include the indicated hyper-methylated candidate biomarkers: heart development (FOXP1, NFATC4, NODAL, TGFBR2 and WNT3A) and vasculature development (MMP2, NODAL and TGFBR2). Genes related to oral/digestive system development are significantly enriched in both CTR and TR analyses. In the CTR signature and three of five CTR network analyses, five pathways or processes consisting of 28 total genes (50% (14) are hyper- or hypo-methylated) were identified. The CTR oral/digestive system development signature contains almost equally divided numbers of genes that encode nuclear/cytoplasmic (8/7) and membrane/extracellular (8/5) proteins. No membrane receptors were identified. Prominent features of this category identified in multiple analyses and associated candidate biomarkers that were hyper- or hypo-methylated include palate development (BCOR, CHD7, MEN1, PDGFRA and TGFB3) and inner ear development (CHD7, CUX1, DLL1, GFI1, PDGFA, PDGFRA, TGFB3 and WNT3A). Both processes are highly enriched with signature genes: palate development contains 5 of 11 (45%) hyper- or hypo-methylated genes, and strikingly, inner ear development contains 10 of 14 (71%) hyper- or hypo-methylated genes. In the TR signature and three TR network analyses, six pathways or processes consisting of 17 total genes (41% (7) are hyper- or hypo-methylated) were identified. The TR oral/digestive system development signature contains a high number of genes that encode nuclear proteins (9) relative to membrane (5), extracellular (2) or cytoplasmic (1) proteins. Prominent features of this category and associated methylated candidate biomarkers include inner ear morphogenesis (WNT3A) and ear development (CUX1, NKX3-2 and WNT3A). All three candidate biomarkers for the oral/digestive system – CUX1, NKX3-2 and WNT3A – are methylated.

The MSR signature, different from CTR and TR signatures, was enriched with two developmental pathways/processes that included one or more signature genes that were differentially methylated. Four of five MSR network analyses identified the WNT signalling pathway. Aberrant WNT signalling is a common feature of tumours and plays important roles in tumour progression and metastasis of many cancer types (Shiah et al. 2016). Nemo-like kinase (NLK) was hypo-methylated and predicted to be upregulated. Aberrant expression of NLK is significantly associated with the initiation and progression of various types of human cancers, as well as clinicopathologic features and survival rate (Huang et al. 2015). Additionally, regulation of cell development was identified in one of five MSR network analyses; Ventral Anterior Homeobox 1 (VAX1) was hyper-methylated and predicted to be down-regulated. Analysis of CpG islands in squamous cell carcinomas and adenocarcinomas of the lung by others has previously identified VAX1 methylated in >80% of adenocarcinomas (Rauch et al. 2012).

Immune response pathways

Pathways and processes associated with an immune response were identified in the CTR signature and two of five CTR network analyses (86 total genes; 43% (37) are hyper- or hypo-methylated) were identified. The CTR immune response signature contains almost equally divided numbers of genes that encode nuclear/cytoplasmic (18/24) and membrane/extracellular (26/18) proteins. Prominent features of this category identified in multiple analyses and associated candidate biomarkers that were hyper- or hypo-methylated include cytokine-mediated signalling (CLCF1, JAK3 and RELA), positive regulation of T cell activation (BCL6, HLA-DOA, IL2, LCK and RARA) and response to hypoxia (CREBBP, CYP1A1, PDGFA, SOCS3 and TGFB3). Notably, all derived expression levels for genes associated with cytokine-mediated signalling were down-regulated while four of five derived expression levels for positive regulation of T cell activation were up-regulated.

In the TR signature, immune response was the most enriched functional category, consisting of seventeen pathways or processes (42 total genes, 45% (19) are hyper- or hypo-methylated). Almost half of the TR immune response signature encodes nuclear proteins (20) while another 26% (11) encode cytoplasmic proteins; eight genes encode membranes and three encode extracellular proteins. The process of T cell differentiation was identified in the TR signature; all four genes were hyper-methylated and predicted to be up-regulated (Table 4). This is consistent with the predicted up-regulation of the CTR genes involved in T cell activation (Table 2). With the exception of one gene – JARID2 – all 18 other signature genes were hyper-methylated. Prominent features of this category and associated candidate biomarkers differentially methylated include response to cytokine stimulus (HDAC4) and response to wounding (HDAC4, NFATC4, PPARD, TGFBR2, TLR1, TNFAIP6, TOLLIP and VWF).

Nine genes involved in the regulation of the transforming growth factor beta receptor signalling pathway were identified in the CTR signature and one of five CTR network analysis (56% (5) are hypo-methylated). Most genes encode nuclear (6) or cytoplasmic (2) proteins. All five associated candidate biomarkers are hypo-methylated: FURIN, MEN1, PRDM16, SKI and TGFB3. Other studies have implicated TGF-b1 in both the inflammatory and remodelling components of early and late persistent respiratory disease (Aubert et al. 1994, Vignola et al. 1997, Takizawa et al. 2001, Wang et al. 2005). In the TR signature, the transforming growth factor beta receptor signalling pathway was enriched in two of six TR network analyses (6 total genes; 17% (1) are hyper-methylated). Five of the six genes encode nuclear proteins; one is a membrane protein and candidate biomarker, transforming growth factor, beta receptor II (TGFBR2). TGFBR2 is hyper-methylated and predicted to be up-regulated in this study. In oesophageal squamous cell carcinoma, TGFBR2 promoter methylation is thought to play an important role in TGFBR2 gene silencing (Dong et al. 2012).

Cell death pathways

Two processes identified in the CTR signature and two of five CTR network analyses are involved with cell death (94 total genes, 36% (34) are hyper- or hypo-methylated). The cell death signature contains a high number of genes that encode cytoplasmic proteins (39) relative to nuclear (26), membrane (17) or extracellular (12) proteins. In particular, regulation of cell death was identified, consisting of 89 genes, 29 of which were differentially methylated and are candidate biomarkers: AHRR, BCL6, BIK, BMP4, CD44, CLCF1, CLU, COL18A1, FURIN, ID3, IGF1, IL2, IL3, LCK, MCF2, MEN1, NGF, NQO1, NTRK1, NUAK2, PPP2CA, PRKCA, RARG, RELA, SOCS2, SOCS3, TCF7L2, TGFB3 and TRAF3.

Three processes identified in the TR signature and five TR network analyses are involved with cell death (38 total genes, 18% (7) are hyper- or hypo-methylated). The TR cell death signature contains a high number of genes that encode cytoplasmic (19) or nuclear (13) relative to membrane (3) or extracellular (3) proteins. Prominent features of this category and associated methylated candidate biomarkers include negative regulation of apoptosis (IL3 and NKX3-2) and induction of apoptosis by extracellular signals (ARHGEF6, BCL2L11, DAPK1, MCF2L and OBSCN).

In the MSR signature, one of five network analyses identified regulation of programmed cell death (64 total genes 2 (3%) are hyper- or hypo-methylated). MCF2 encodes an oncogenic protein that is a guanine nucleotide exchange factor (GEF), which exerts control over some members of the Rho family of small GTPases. Similarly, as part of the Ragulator complex, LAMTOR5 is involved in amino acid sensing and activation of mTORC1. Ragulator functions as a GEF activating the small GTPase Rag.

Cell cycle pathways

One process identified in the CTR signature is involved with negative regulation of the cell cycle (6 total genes, 100% (6) are hyper- or hypo-methylated), consisting of BCL6, BMP4, DGKZ, MEN1, MYO16 and NGF. In the TR signature, three of six TR network analyses identified cell cycle arrest (10 total genes, 3 (30%) are hyper- or hypo-methylated). All three genes are hyper-methylated and predicted to be down-regulated: CUL1, DST and RASSF1. An association was recently shown between 13 genes involved in cell cycle gene transcription, including CUL1 and tumour size in oral squamous cell carcinoma (Diniz et al. 2015). Another recent study found that RASSF1 inhibits YAP activation through the GEF-H1/RhoB pathway, suppressing the invasion and metastatic potential of non-small cell lung cancer cells (Dubois et al. 2016).

All five MSR network analyses-identified processes are involved in cell cycle regulation (60 total genes, 5% (3) are hyper- or hypo-methylated). A prominent feature of this category is mitosis. The gene encoding cytoplasmic linker associated protein 2 (CLASP2), which promotes the stabilization of dynamic microtubules and fibronectin type III and SPRY domain containing 1 (FSD1), which associates with a subset of microtubules and may be involved in the stability and organization of microtubules during cytokinesis – were hypo-methylated and predicted to be down-regulated. A third mitosis gene, kinesin family member 15 (KIF15), maintains bipolar microtubule spindle apparatus in dividing cells and was methylated and predicted to be down-regulated.

Discussion

The main objective of this work is to identify methylation changes in consumers of combustible and non-combustible tobacco, particularly MSC. Cigarette smoking has been known to induce global methylation changes, which have been suggested to contribute to the smoking-related diseases, including COPD, cancer and CVD (Lee and Pausova 2013). While the existing US epidemiology indicates that smokeless/MSC are at significantly lower risk, compared to cigarette smokers, for the smoking-related diseases, the molecular effects, such as gene methylation differences, due to smokeless tobacco product is incompletely understood. In this manuscript, we show that cigarette smoking induces marked DNA methylation changes in the buccal cells, compared to those detected in MSC.

Consistent with the published data (mostly derived from blood), we find that buccal cells from SMK exhibit prominent changes in DNA methylation. The top 20 differentially methylated gene loci were from SMK. Hierarchical clustering of the top 50 differentially methylated signatures revealed that SMK tended to segregate from the MSC and NTC, whereas NTC and MSC co-cluster (Figure 3). We have grouped the differentially methylated genes into three signatures: CTR, MSR and TR. Among these, 22 genes exhibited multiple methylation changes (>2) in CTR, and 12 genes were differentially methylated in more than 2 loci. The few MSR genes found to be differentially methylated were altered at only one locus.

The CTR signature, which appears to be exclusively derived from the SMK, was found to be the largest and consisted of gene loci, previously reported in smokers. For example AHRR gene was differentially methylated in the CTR at 7 loci (Table 1).

The AHRR gene has been shown to be differentially methylated in previous studies of smokers (Joubert et al. 2012, Monick et al. 2012, Philibert et al. 2012, Shenker et al. 2013, Besingi and Johansson 2014, Joehanes et al. 2016). The most significant site in our study (cg05575921) was also the most significant in previous studies (Joubert et al. 2012, Sun et al. 2013, Zeilinger et al. 2013, Besingi and Johansson 2014, Joehanes et al. 2016). AHRR is a negative regulator of the aryl hydrocarbon receptor (AHR) gene that codes for a protein that binds to a wide range of xenobiotics. AHR induces the expression of CYP1A1 and other genes involved in the metabolism of xenobiotics such as dioxin and polycyclic aromatic hydrocarbons (PAHs), as well as nicotine. AHRR has been proposed to be a tumour suppressor gene in multiple types of human cancers (Zudaire et al. 2008). Optimal functioning of the xenobiotic detoxification system requires retinoids for sensing, detoxifying, and eliminating xenobiotics (Shmarakov 2015), and we find that response to the retinoid vitamin A is impacted in the CTR signature.

A number of CTR genes with multiple hits are associated with diseases linked to combustible tobacco use, including atherosclerosis and arteriosclerosis (AHRR, CYP1A1, ID3, IGF1, NR2F2 and TSLP) and coronary heart disease (CYP1A1, DLEU7, IGF1, NR2F2 and TCN2).

The top hits in the TR gene loci were all hyper-methylated (Table 1), and the SEPT9 and RASSF1 were more prominent among them, both of which were differentially methylated at 4 loci. Consistent with our findings, RASSF1A (RASSF1) has been previously shown in a large case-control study of lung cancer to be influenced by sex, with males showing higher levels of methylation (Vaissiere et al. 2009). In head and neck squamous cell carcinoma, a trend was previously shown between SEPT9 hyper-methylation and HPV16 positivity (Bennett et al. 2010). HPV16 is associated with oral papillomas and oropharyngeal cancer, as well as other tumour types and conditions. Thus, the hyper-methylation of tumour suppressor genes in the SMK is in agreement with the known epidemiology that cigarette smokers are at higher risk for cancers, including lung and oral cancer (U.S. Department of Health and Human Services 2010). Additional work is necessary to understand the contribution of exposure to combustion-related toxicants and tobacco exposure to the methylation changes in the TR signature.

The buccal cell methylation data from MSC suggests that microtubules, cell cycle and cell death are impacted in that cohort. This result is consistent with studies using smokeless tobacco extract, where toxicity was mediated by disruption of the microtubule network (Das et al. 2013). Based on the few loci detected in the MSR methylation signature, which are exclusive to MSC, it is possible TR methylation signature is predominantly driven by changes from combustible tobacco use. A previous report indicates that no differentially methylated loci/signatures were detected in the blood of smokeless tobacco users (Besingi and Johansson 2014), potentially indicating tissue-specific differences in the methylation patterns in tobacco users.

Combustible tobacco drives profound changes in buccal cell methylation status, principally impacting cell developmental and immune response pathways (Table 2) (Besingi and Johansson 2014). This is consistent with emerging evidence supporting physiological functions of AHR in these pathways (Hao and Whitelaw 2013). In contrast, very few differentially methylated loci were detected in buccal cells from MSC and the top 50 significant differential methylation gene loci co-cluster MSC and NTC (Figure 3). Similar to other studies in both buccal cells and blood, AHRR was the most significantly hypo-methylated gene in SMK (Shenker et al. 2013, Zeilinger et al. 2013, Besingi and Johansson 2014, Teschendorff et al. 2015). However, based on the location of differential methylation of a locus (Figure 4, cg17924476), AHRR expression may be up-regulated in buccal cells of SMK, and the tissue-specific level of gene expression requires further confirmation. Separately, we have shown that AHRR gene expression is upregulated in the peripheral blood mononuclear cells of SMK, relative to MSC and NTC from the same cohort (Arimilli et al. 2017). Taken together, the AHRR pathway does not appear to be activated in MSC, suggesting that xenobiotic metabolism is differentially regulated between smokers and smokeless tobacco users.

A potential limitation of this study is that the classification of signatures is mutually exclusive. Thus, a differentially methylated gene is not present in more than one signature. However, our data from cigarette smokers is consistent with other studies and it is likely that this limitation simply reduced the overall number of loci detected for a given signature and did not influence biological interpretation. In particular, this may have impacted the number of MSR loci identified, which is limited in this analysis and warrants additional study. Similarly, additional work will be necessary to understand the nature of the TR gene loci.

A different study from Sweden, which evaluated blood global methylation profiles from smokers and snuff consumers, did not find any significant methylation differences in the snuff consumers, and suggested that the combustion-derived toxicants, rather than tobacco was the driver of epigenetic changes (Besingi and Johansson 2014). Those findings are similar to those of ours, which showed minimal epigenetic alterations in buccal cells of MSC; smokers in both studies, however, experienced marked DNA methylation alterations.

In summary, we have for the first time shown global methylation in buccal cells of MSC, and compared with those profiles with global methylation differences in smokers. Relative to SMK, only a limited number of gene loci were differentially methylated in MSC from buccal cell DNA. Our findings aid in placing combustible and non-combustible tobacco products along a risk continuum related to tobacco use and provides additional insights into the effects of tobacco consumption.

Conclusion

In this study, we show differential methylation changes in buccal cells due to cigarette smoking and the use of moist snuff (used as a representative of smokeless tobacco products). The gene methylation differences were classified into three groups: CTR, TR and MSR methylation signatures. Consistent with published studies, healthy cigarette smokers’ exhibit pronounced alterations in the methylation status of several genes, which could be linked to the development of smoking-related diseases. Relatively, fewer differences detected in the methylation of genes in MSC (as evident from MSR signature), compared to smokers. Further research is needed to better understand the methylation changes in the TR signature. Methylation status of genes, such as AHRR is useful potential biomarkers to distinguish the effects of combustible and smokeless tobacco.

Acknowledgements

The authors thank Bobbette A. Jones for managing the clinical conduct of the study and for assistance in preparing this manuscript. The authors thank Megan J. Whelen for critical review and editing of this manuscript. Drs. Borgerding and Prasad are full-time employees of RAI Services Company. RAI Services Company is a wholly owned subsidiary of Reynolds American Inc., which is a wholly owned subsidiary of British American Tobacco Plc. Dr. Jessen is a full-time employee of Laboratory Corporation of America Holdings (LabCorp).

Disclosure statement

Drs. Borgerding and Prasad are full time employees of RAI Services Company. RAI Services Company is a wholly owned subsidiary of Reynolds American Inc., which is a wholly owned subsidiary of British American Tobacco plc. Dr. Jessen is a full time employee of Laboratory Corporation of America Holdings (LabCorp).