Epigenome-450K-wide methylation signatures of active cigarette smoking : The Young Finns Study

Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 *These authors contributed equally to this work. Received: 03 March 2020 Revised: 11 June 2020 Accepted: 23 June 2020 Accepted Manuscript online: 25 June 2020 Version of Record published: 07 July 2020 Research Article Epigenome-450K-wide methylation signatures of active cigarette smoking: The Young Finns Study Pashupati P. Mishra1,2,3,*, Ismo H¨ anninen1,2,3,*, Emma Raitoharju1,2,3, Saara Marttila1,2,3,4, Binisha H. Mishra1,2,3, Nina Mononen1,2,3,MikaK¨ ah¨ onen2,5, Mikko Hurme4,6, Olli Raitakari7,8,9, Petri T ¨ or¨ onen10, Liisa Holm10,11 and Terho Lehtim¨ aki1,2,3 1Department of Clinical Chemistry, Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland; 2Finnish Cardiovascular Research Center-Tampere, Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland; 3Department of Clinical Chemistry, Fimlab Laboratories, Tampere, Finland; 4Gerontology Research Center (GEREC), Tampere University, Tampere, Finland; 5Department of Clinical Physiology, Tampere University Hospital, Tampere, Finland; 6Department of Microbiology and Immunology, Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland; 7Centre for Population Health Research, University of Turku and Turku University Hospital, Turku, Finland; 8Research Centre of Applied and Preventive Cardiovascular Medicine, University of Turku, Turku, Finland; 9Department of Clinical Physiology and Nuclear Medicine, Turku University Hospital, Turku, Finland; 10Institute of Biotechnology, Helsinki Institute of Life Sciences (HiLife), University of Helsinki, Helsinki, Finland; 11Organismal and Evolutionary Biology Research Program, Faculty of Biological and Environmental Sciences, University of Helsinki, Helsinki, Finland Correspondence: Pashupati P. Mishra (pashupati.mishra@tuni.fi) Smoking as a major risk factor for morbidity affects numerous regulatory systems of the human body including DNA methylation. Most of the previous studies with genome-wide methylation data are based on conventional association analysis and earliest threshold-based gene set analysis that lacks sensitivity to be able to reveal all the relevant effects of smoking. The aim of the present study was to investigate the impact of active smoking on DNA methylation at three biological levels: 5-C-phosphate-G-3(CpG) sites, genes and functionally related genes (gene sets). Gene set analysis was done with mGSZ, a modern threshold-free method previously developed by us that utilizes all the genes in the experiment and their differential methylation scores. Application of such method in DNA methylation study is novel. Epigenome-wide methylation levels were profiled from Young Finns Study (YFS) participants’ whole blood from 2011 follow-up using Illumina Infinium HumanMethylation450 BeadChips. We identified three novel smoking related CpG sites and replicated 57 of the previously identified ones. We found that smoking is associated with hypomethylation in shore (genomic regions 0–2 kilobases from CpG island). We identified smoking related methylation changes in 13 gene sets with false discovery rate (FDR) ≤0.05, among which is olfactory receptor activity, the flagship novel finding of the present study. Overall, we extended the current knowledge by identifying: (i) three novel smoking related CpG sites, (ii) similar effects as aging on average methylation in shore, and (iii) a novel finding that olfactory receptor activity pathway responds to tobacco smoke and toxin exposure through epigenetic mechanisms. Background Tobacco use has been estimated to cause approximately 7 million yearly deaths globally through general morbidity, disability, non-communicable chronic diseases, and increased risk of communicable diseases [1]. The chronic diseases caused by smoking are in nature multifactorial, mediated by both genetic and environmental factors [2,3]. Despite the well-established link between smoking and associated diseases in literature, there is a knowledge gap in the causative mechanisms of the diseases. Previous studies have led to identification of several transcriptomic signatures and their role in diseases related to smoking [4]. Epigenetic mechanisms have been understood to be modulated by environmental ©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). 1 Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 factors such as smoking and play a crucial role in disease development process by regulating gene expression [5]. Several studies have reported smoking-related methylation changes at CpG site, gene as well as gene set level [6–11]. Zeilinger et al. reported 972 CpG sites significantly associated with smoking from discovery data with 1793 participants and replicated 187 out of 972 CpG sites in independent cohort of 479 participants [6]. The three most significant CpG sites and corresponding genes were cg05575921:AHRR, cg21566642: ALPP/ALPPL2, and cg03636183:F2RL3. McCartney et al. [7] and Zeilinger et al. [6] suggested that smoking-induced DNA methylation changes are reversible. Shenker et al. developed methylation index by combining four smoking related CpG sites (one in AHRR, two intergenic in chromosome 2 and one intergenic in chromosome 6) that can predict former smoking status [8]. The study by Joehanes et al., one of the largest study of smoking-related DNA methylation changes, identified 18760 active smoking-related CpG sites annotated to 7201 genes with FDR <0.05 [9]. Prince et al. investigated association of 2620 previously reported smoking-related CpG sites with different smoking behaviors (such as ever smoking, current weekly smoking, ever weekly smoking, and blood cotinine levels) and reported smoking behavior related methylation patterns in 11 CpG sites mapping to seven genes [10]. Several studies on impact of aging on DNA methylation, as summarized by Ciccarone et al., have established age-induced DNA methylation changes as hallmark of aging [12]. Recent studies such as by Yang et al. [13] and Lei et al. [14] have reported association between cigarette smoking and age-induced DNA methylation changes. Joehanes et al. also performed Gene Ontology (GO) [15] based gene set analysis of the list of statistically significant CpG sites and reported 99 gene sets that broadly included molecular processes such as signal transduction, protein metabolic process, and transcription pathways [9]. Another study performed gene set analysis of list of CpG sites with altered methylation levels due to exposure to maternal smoking during pregnancy and reported cell cycle, cancer, white blood cell differentiation, genotoxicity as major findings [16]. Bakulski et al. performed gene set analysis of lung-specific smoking related CpG sites and the most significant pathways included mRNA catabolic processes, protein targeting, angiogenesis, and mRNA translation [17]. Traditional epigenome-wide association studies might fail to identify relevant CpG sites due to noise and lack of statistical power. Further, long list of CpG sites discovered through such studies are harder to interpret in terms of underlying biological theme. Gene set analysis, also referred to as pathway analysis, addresses the problems by shifting the analysis from individual CpG sites to groups of related CpG sites. CpG sites are mapped to respective genes and grouped together based on shared biological features. The grouping of corresponding genes can be based on reference knowledge bases such as Reactome Pathway Database [18], BioCarta pathways [19], Kyoto Encyclopedia of Genes and Genomes (KEGG) [20], and GO [15]. The grouping approach averages out errors at single CpG site level and increases statistical power. Additionally, the outcome of gene set analysis is meaningful terms that provide insights into biological theme underlying the list of CpG sites. Studies that have reported smoking-related methylation changes at biological pathway level have done so mostly with the earliest gene set analysis approach, over-representation analysis (ORA) [9,17,21]. ORA is a type of gene set analysis that is based on statistical tests like hypergeometric and chi-square. The approach requires users to provide a list of genes selected based on arbitrary threshold (for example, P-value<0.05). The arbitrariness of the threshold leads to unstable results that are difficult to replicate and massive reduction in sensitivity. Furthermore, the approach assumes independence between genes, which is violated in most cases resulting to false positive results [22]. Over a decade of development of gene set analysis approaches have provided several robust gene set analysis methodsthatarelikelytoprovideimprovedinsightsintoepigeneticconsequencesofsmoking.Aclassofsuchmethods includes threshold-free competitive gene set analysis methods that test whether genes in a given gene set are more differentially methylated than the other genes in the dataset. The approach can be used to test a large number of gene sets, such as those provided by Molecular Signatures Database (MSigDB) [23,24] in order to profile smoking-related epigenetic alterations at biological pathways or processes level. An example of one of such methods is mGSZ (modified Gene Set Z-score) that have been developed and improved over time by us [25–27] and have been shown to be robust and efficient by us [20–22] as well as another independent study [28]. Unlike over-representation analysis methods, mGSZ is threshold-free and thus analyzes all the genes in the data resulting into sensitive, comprehensive, andreproducibleresults.ThemGSZ is different from other methods in the class in that it is based on asymptotic P-value estimation that is significantly more efficient than empirical P-value estimation [26,29]. The objective of the present study was to characterize DNA methylation differences between current smokers and non-smokers at different biological levels such as: (1) CpG site level with differential methylation analysis; (2) biological pathway level using a threshold-free gene set analysis method, mGSZ; and (3) average methylation level at different genomic regions with Wilcoxon rank-sum test. We show that application of threshold-free gene set analysis method for DNA methylation data has improved benefit over earliest threshold-based ORA methods used currently by most studies as illustrated by our novel findings such as smoking-induced altered methylation in olfactory system. 2©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Table 1 Population characteristics of Young Finns Study participants; data are mean + −SD or proportions Characteristics Active smokers Non-smokers Number of subjects 21 104 Sex (%men) 57% 34% Age, years 45 (+ −3.95) 44.1 (+ −3.1) Body mass index, kg/m226.59 (+ −4.38) 25.7 (+ −4.36) Even though olfactory ability is known to be adversely affected by active cigarette smoking [30], we for the first time show the underlying epigenetic mechanism. Materials and methods Study population Young Finns Study (YFS) is one of the largest existing follow-up studies into cardiovascular health from childhood to adulthood, running in a longitudinal prospective setup with regular follow-ups from 1980 onwards [31]. The study began in 1980 with 3596 children and adolescents aged 3 to 18 years randomly selected from the areas of five university hospitals in Finland (Turku, Tampere, Helsinki, Kuopio, and Oulu). The participants have been followed up for over 40 years. The methylation measurements were performed on a subset of 192 individuals from whole-blood samples from 2011 follow up. The smoking history of the subjects was self-reported and belonged to six categories based on smokingfrequency(1.activesmokeroratleastonceaday,2.onceaweekormoreoften,howevernotdaily,3.lessoften thanonceaweek,4.attemptstoquit,5.hasquit,6.hasneversmoked).Thepresentstudyisbasedonasubsampleof 125 participants, 40–49 years of age who were either active smoker (n=21) or have never smoked (n=104) [Table 1]. Participants in the middle conditions such as those who smoked less frequently or has quit smoking were eliminated in order to obtain maximum possible smoking-related biological signal in DNA methylation. DNA methylation assessment DNA was extracted from EDTA-blood samples using Wizard®Genomic DNA Purification Kit (Promega Corporation, Madison, WI, U.S.A.), according to the manufacturer’s instructions. Genome-wide quantification of DNA methylation levels were done using Illumina Infinium HumanMethylation450 BeadChips [32] in the Core Facility at the Institute of Molecular Medicine Finland, University of Helsinki, following manufacturer’s protocols. The HumanMethylation450 BeadChip measures DNA methylation at more than 485,000 CpG cites across the genome. Thearrayswereimagedwithahigh-precisionscanner(iScansystem,IlluminaInc.),andthesignalintensitieswere extracted using a GenomeStudio software package (Illumina Inc.). The methylation data are available in Gene Expression Omnibus (GEO) under accession number GSE69270. Data filtering and normalization Data were obtained and processed from raw methylation image files using the minfi package in R/Bioconductor [33].AlltheanalyzedsampleshavesumofdetectionP-values across all the probes less than 0.05. Logged (log2) median of methylated and unmethylated intensities of the analyzed samples clustered within the default threshold (10.5) of getQC function in minfi. Further, samples for which real sex did not match the predicted sex obtained with getSex function in minfi were excluded. Background subtraction and dye-bias normalization were performed via noob method [34] implemented in minfi. Stratified quantile normalization was performed using preprocessQuantile function in minfi. Probes with detection P-value more than 0.01 in 99% of the samples were filtered out. CpG loci on sex chromosomes were excluded from the analysis to avoid gender-based methylation bias. Also, cross-reactive probes and probes with SNPs were excluded from the analysis. After quality control, the total number of autosomal CpGs was 429,773 in 125 samples (21 active smokers and 104 non-smokers). Batch correction was done by adding the first two principal components of control probes as covariates in multiple linear regression model. Differential methylation analysis All statistical analyses were performed using R statistical software (v.3.5.1) [35]. M-values, calculated as the log2 ratio of the intensities of methylated probe versus unmethylated probe were used as measures of methylation level. Differentially methylated CpG loci for smoking status were identified using CpGassoc packageinR[36].Inorder to keep the present study exploratory in nature, analyses were based on two different models adjusted with varying ©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). 3 Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 number of covariates. Model 1 involved adjustment for age, sex, body mass index, cell type proportions, and technical covariates (chip and array). Cell type proportions consisted of proportions of CD8T, CD4T, natural killer cells, B cells, monocytes, and granulocytes in white blood cells. Selection of covariates for model 1 was motivated from study by [9]. Model 2 involved adjustment with two additional covariates, alcohol usage and socioeconomic status, in addition to those in model 1. The addition of the two covariates in model 2 is based on previous reports suggesting effects of alcohol usage [37] and socioeconomic status [38] on DNA methylation. There was reduction in sample size for model 2 because 16 out of 125 participants did not have information on alcohol usage and socioeconomic status. The number of active smokers in the analysis was reduced from 21 to 16 because of missing information on socioeconomic status. Alcohol consumption was measured by asking participants to report their alcohol consumption during the previous week. One unit is equivalent to 14 g of alcohol [39]. Socioeconomic status was based on occupation and was categorized as manual, lower non-manual, and upper non-manual. Potential batch affects were addressed by including the first two principal components of array control probes into the regression models. The cell type proportions of white blood cells were estimated through the reference-based Houseman method [40] using the estimateCellCounts function in the minfi BioconductorpackageinR.CpGsitesweremappedtogenesbyusing annotation database provided by Illumina [41]. Differentially methylated genes were identified by utilizing as a proxy the CpG site with maximum absolute t-score from any location in the genomic region of the gene, since the mechanism how the methylation influences gene expression is not completely understood. Statistical significance level was set to false discovery rate (FDR) of 0.05 in the CpG site and gene level analysis. Gene set analysis Biological relevance of the differentially methylated genes based on both model 1 and 2 was investigated using mGSZ method implemented in mGSZ R package [26]. mGSZ is a gene set analysis method based on robust gene set scoring function and efficient P-value estimation method. Unlike over-representation-based pathway analysis methods [42,43], this approach is threshold free and thus includes all the genes in the analysis irrespective of their effect size or significance level. This is particularly important in the present study as our goal is to identify set of genes contributing to a biological pathway that have milder but coherent smoking-related changes in methylation level. Analyzed gene sets were downloaded from MSigDB version 7.0. The database contained 22569 gene sets (as of September 21, 2019) divided into eight major collections: Hallmark gene sets, positional gene sets, curated gene sets, motif gene sets, computational gene sets, gene ontology (GO) gene sets, oncogenic signatures, and immunogenic signatures. Hallmark gene sets are generated computationally and represent well-defined biological processes [44]. Positional gene sets are generated based on genomic locations of genes. Curated gene sets are generated from knowledge sources such as pathway databases, biomedical literature, and domain experts. The gene sets are divided into two sub-categories: (i) canonical pathways curated from online databases [18–20] and (ii) chemical and genetic perturbations (CGP). CGP-based gene sets represent expression signatures of chemical and genetic perturbations mostly curated from biomedical literature [45]. Motif gene sets represent target genes for transcription factors or micro RNAs. Computational gene sets are generated from cancer related microarray gene expression data with data mining techniques. GO gene sets are derived from GO annotations and have three sub-categories—biological process (BP), molecular function (MF), and cellular component (CC). Average methylation level at different genomic regions We analyzed effect of smoking on the average methylation level at different genomic regions in relation to both gene region (TSS200, 0–200 bases upstream of the transcriptional start site; TSS1500, 200–1500 bases upstream of the TSS; 5UTR, within the 5untranslatedregion,betweentheTSSandtheATGstartsite;Body,betweentheATGandstop codon irrespective of the presence of introns, exons, TSS, or promoters; 3UTR, between the stop codon and poly A signal) and CpG islands (Shores, 0–2 kb from CpG island; Shelves, 2–4 kb from CpG island). The information about the different genomic regions of CpG sites on HumanMethylation450K array was obtained from Illumina Inc. Mean of methylation β-values of all CpGs belonging to the genomic regions for each of the study participants was calculated. Smoking-related statistical difference in average methylation specific to each of the regions was tested with Wilcoxon rank-sum test. Results Study participants’ characteristics CharacteristicsoftheYFScohortparticipantsinthepresentstudyissummarizedinTable1. 4©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Figure 1. Manhattan plot and quantile-quantile (Q-Q) plot of epigenome-wide association analysis of smoking habit (A) Manhattan plot showing the P-values of genome-wide CpG sites. X-axis represents position of the CpG sites on each chromosome. Y-axis represents negative log10 of the P-values for the association. The dotted line indicates false discovery rate (FDR)-corrected significance threshold and the solid horizontal line represents Bonferroni-corrected significance threshold (experiment-wide significance). The annotations for the three novel CpG sites and corresponding genes (if annotation is known) are shown. (B) Quantile-quantile plot showing genomic inflation factor (lambda =1.11) of the epigenome-wide association study. The genomic inflation factor (ratio of the median of the empirically observed distribution of the test statistic to the expected median) represents the extent of inflation and false positive rate in the results. Differential methylation analysis We identified 60 differentially methylated CpG sites with respect to smoking habit with FDR <0.05 from statistical model 1 [Table 2]. There were 26602 nominally significant (P-value <0.05) CpG sites (Figure 1). The results identify three novel smoking-associated CpG sites: cg26038589 (CCDC55, Coiled-Coil Domain-Containing Protein 55), cg10385208 (CWC25, Spliceosome Associated Protein Homolog), and cg09355027 (no known gene) and replicate the findings of the previous studies, indicating that our results are technically robust. In addition, a site found in only one previous study, cg13898430 (RUNX3, Runt-related transcription factor 3) [46] was also replicated in the present study. Model 2 identified only 18 differentially methylated CpG sites with respect to smoking habit with FDR <0.05. Seventeen of the CpG sites from model 2 that were also reported by model 1 are highlighted in Table 2 with bold font. The three novel sites identified in model 1 were not recovered with model 2. Gene set analysis Using the gene sets from MSigDB and mGSZ, we identified smoking induced altered methylation in a total of 13 gene sets with FDR ≤0.05. The results are presented in Table 3 and Figure 2. The figure serves the purpose of observing overall methylation pattern and we do not expect significant statistical differences with such approach due to sample size limitation. This is where our proposed state-of-the-art gene set analysis method, mGSZ,becomesusefulasthe approach increases the statistical power and identifies biologically relevant results even with smaller sample size. Among eight major gene set collections of MSigDB, significant results (FDR ≤0.05) were obtained for three—Hallmark, Curated, and GO gene sets [Table 3]. The gene set, Genes down-regulated in response to ultraviolet (UV) radiation (FDR =0.05) identified from Hallmark collection suggests that smoking induces hypermethylation of genes that are known to be down-regulated by ultraviolet radiation [Table 3 and Figure 2]. We identified three canonical pathways with smoking related altered methylation (hypomethylation): two from reactome database and one from BioCarta. Gene sets based on reactome pathways are aryl hydrocarbon receptor signaling (FDR = 0.003) and roraactivates gene expression (FDR =0.03). Similarly, BioCarta based results suggest that smoking alters ©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). 5 Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Table 2 List of differentially methylated CpG sites with respect to smoking habit with false discovery rate (FDR) <0.05 identified with model 1, their corresponding effect size, standard error, P-values, FDR, and genomic regions in relation to both gene region (TSS200, 0–200 bases upstream of the transcriptional start site; TSS1500, 200–1500 bases upstream of the TSS; 5UTR, within the 5untranslated region, between the TSS and the ATG start site; Body, between the ATG and stop codon irrespective of the presence of introns, exons, TSS, or promoters; 3UTR, between the stop codon and poly A signal) and CpG islands (Shores, 0–2 kb from CpG island; Shelves, 2–4 kb from CpG island) CpG Sites Effect size Standard error P-value FDR Gene/Island name Gene region/Relation to island Chromosome 1 cg20295214 0.03 0.005 9.53 ×10-9 1.84 ×10-4 AVPR1B Body 1 2 cg08709672 0.02 0.004 3.21 ×10-6 2.82 ×10-2 AVPR1B 5UTR 1stExon 1 3 cg09935388 0.06 0.009 9.87 ×10-10 2.65 ×10-5 GFI1 Body 1 4 cg18316974 0.03 0.006 7.64 ×10-7 8.87 ×10-3 GFI1 Body 1 5 cg18146737 0.05 0.01 5.55 ×10-6 4.18 ×10-2 GFI1 Body 1 6 cg25189904 0.03 0.006 3.85 ×10-6 3.18 ×10-2 GNG12 TSS1500 1 7cg21408581 0.02 0.005 2.18 ×10-6 2.03 ×10-2 RAB3GAP2 Body 1 8 cg08869700 0.01 0.003 2.75 ×10-6 2.46 ×10-2 RFX5 Body 1 9 cg13898430 0.008 0.002 5.07 ×10-6 3.96 ×10-2 RUNX3 TSS1500 1 10 cg04885881 0.02 0.004 1.52 ×10-7 2.43 ×10-3 chr1:11119030 -11120634 Shelf 1 11 cg27537125 0.02 0.003 9.85 ×10-9 1.84 ×10-4 – OpenSea 1 12 cg12547807 0.01 0.002 6.36 ×10-7 7.84 ×10-3 – OpenSea 1 13 cg11555067 0.02 0.003 3.42 ×10-6 2.94 ×10-2 INPP4A 5UTR 2 14 cg21566642 0.03 0.002 1.02 ×10-23 2.19 ×10-18 chr2:233283397 -233285959 Island 2 15 cg01940273 0.05 0.004 2.46 ×10-21 3.53 ×10-16 chr2:233283397 -233285959 Island 2 16 cg05951221 0.05 0.004 8.74 ×10-19 7.51 ×10-14 chr2:233283397 -233285959 Island 2 17 cg03329539 0.02 0.003 2.71 ×10-9 6.86 ×10-5 chr2:233283397 -233285959 Shore 2 18 cg06644428 0.06 0.01 9.64 ×10-7 1.09 ×10-2 chr2:233283397 -233285959 Island 2 19 cg21949194 0.02 0.003 6.06 ×10-6 4.49 ×10-2 chr2:39351355 -39351733 Shelf 2 20 cg23079012 0.02 0.003 5.43 ×10-10 1.67 ×10-5 – OpenSea 2 21 cg09355027 (novel) 0.02 0.004 1.43 ×10-6 1.46 ×10-2 – OpenSea 2 22 cg19859270 0.03 0.004 9.59 ×10-12 3.75 ×10-7 GPR15 1st Exon 3 23 cg24719910 0.02 0.003 5.41 ×10-6 4.15 ×10-2 TGFBR2 Body 3 24 cg05575921 0.10 0.006 8.53 ×10-28 3.66 ×10-22 AHRR Body 5 25 cg21161138 0.04 0.004 1.09 ×10-19 1.17 ×10-14 AHRR Body 5 26 cg25648203 0.04 0.004 5.82 ×10-13 3.57 ×10-8 AHRR Body 5 27 cg14817490 0.05 0.005 2.01 ×10-12 9.60 ×10-8 AHRR Body 5 28 cg26703534 0.03 0.003 4.53 ×10-12 1.95 ×10-7 AHRR Body 5 29 cg24090911 0.03 0.005 1.38 ×10-10 4.56 ×10-6 AHRR Body 5 30 cg12806681 0.04 0.006 1.44 ×10-8 2.59 ×10-4 AHRR Body 5 31 cg03991871 0.05 0.009 1.12 ×10-7 1.85 ×10-3 AHRR Body 5 32 cg04551776 0.02 0.004 1.68 ×10-7 2.58 ×10-3 AHRR Body 5 33 cg23916896 0.05 0.009 3.73 ×10-7 5.17 ×10-3 AHRR Body 5 34 cg11554391 0.02 0.004 1.21 ×10-6 1.27 ×10-2 AHRR Body 5 35 cg01899089 0.02 0.003 4.16 ×10-6 3.31 ×10-2 AHRR Body 5 36 cg11902777 0.05 0.009 6.94 ×10-6 4.97 ×10-2 AHRR Body 5 37 cg06126421 0.06 0.007 1.04 ×10-12 5.61 ×10-8 – OpenSea 6 38 cg24859433 0.03 0.004 6.12 ×10-11 2.19 ×10-6 – OpenSea 6 39 cg14753356 0.02 0.003 3.15 ×10-9 7.52 ×10-5 – OpenSea 6 40 cg15342087 0.02 0.003 6.49 ×10-9 1.39 ×10-4 – OpenSea 6 Continued over 6©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Table 2 List of differentially methylated CpG sites with respect to smoking habit with false discovery rate (FDR) <0.05 identified with model 1, their corresponding effect size, standard error, P-values, FDR, and genomic regions in relation to both gene region (TSS200, 0–200 bases upstream of the transcriptional start site; TSS1500, 200–1500 bases upstream of the TSS; 5UTR, within the 5untranslated region, between the TSS and the ATG start site; Body, between the ATG and stop codon irrespective of the presence of introns, exons, TSS, or promoters; 3UTR, between the stop codon and poly A signal) and CpG islands (Shores, 0–2 kb from CpG island; Shelves, 2–4 kb from CpG island) (Continued) CpG Sites Effect size Standard error P-value FDR Gene/Island name Gene region/Relation to island Chromosome 41 cg12803068 -0.09 0.01 2.41 ×10-7 3.45 ×10-3 MYO1G Body 7 42 cg04180046 -0.03 0.005 6.39 ×10-7 7.84 ×10-3 MYO1G Body 7 43 cg22132788 -0.05 0.007 1.58 ×10-6 1.58 ×10-2 MYO1G Body 7 44 cg07826859 0.01 0.003 1.94 ×10-6 1.89 ×10-2 MYO1G TSS1500 7 45 cg21322436 0.01 0.003 4.06 ×10-6 3.29 ×10-2 CNTNAP2 TSS1500 7 46 cg03450842 0.01 0.003 4.42 ×10-7 5.93 ×10-3 ZMIZ1 5UTR 10 47 cg01901332 0.02 0.004 1.19 ×10-6 1.27 ×10-2 ARRB1 Body 11 48 cg21611682 0.02 0.003 7.56 ×10-10 2.17 ×10-5 LRP5 Body 11 49 cg14624207 0.01 0.003 6.91 ×10-6 4.97 ×10-2 LRP5 Body 11 50 cg07986378 0.02 0.005 2.49 ×10-6 2.28 ×10-2 ETV6 Body 12 51 cg22851561 0.02 0.003 3.90 ×10-8 6.71 ×10-4 C14orf43 5UTR 14 52 cg05284742 0.02 0.003 4.53 ×10-9 1.03 ×10-4 ITPK1 Body 14 53 cg07069636 0.01 0.002 4.82 ×10-7 6.28 ×10-3 chr16:30669107 -30671155 Shore 16 54 cg26038589 (novel) 0.03 0.006 1.17 ×10-6 1.27 ×10-2 CCDC55 Body 17 55 cg10385208 (novel) 0.01 0.003 3.81 ×10-6 3.18 ×10-2 CWC25 TSS1500 17 56 cg19572487 0.02 0.004 8.07 ×10-9 1.65 ×10-4 RARA 5UTR 17 57 cg03636183 0.05 0.004 1.00 ×10-17 7.17 ×10-13 F2RL3 Body 19 58 cg03707168 0.02 0.004 7.10 ×10-7 8.47 ×10-3 PPP1R15A Body 19 59 cg17566560 0.02 0.004 2.13 ×10-6 2.03 ×10-2 DLGAP4 Body, 5UTR 20 60 cg23110422 0.03 0.006 1.85 ×10-7 2.74 ×10-3 ETS2 Body 21 Intergenic chromatin regions are called OpenSea. CpG sites retained in model 2 are highlighted with bold font. Statistical models: Model 1 involved adjustment for age, sex, body mass index, cell type proportions, and technical covariates (chip and array). Model 2 involved adjustment with two additional covariates, alcohol usage and socioeconomic status, in addition to those in model 1. methylation in genomic region responsible for RARRXR pathway (FDR =0.05), which is linked to cancer. CGP-based gene set identified in the present study suggests altered methylation in genes that are up-regulated in glioma cell lines after knockdown of SPARC gene by RNAi (FDR =0.003). GO annotation-based gene sets were analyzed separately for the three GO categories. We identified one BP and 7 MF related gene sets with FDR ≤0.05 [Table 3]. Olfactory receptor activity, a novel finding of this study, belongs to GO MF molecular. Gene set analysis of differentially methylated CpG sites identified with model 2 reported seven gene sets with FDR <0.25 [Supplementary Table S2]. The novel finding of the present study, smoking-induced altered methylation in olfactory system, was retained with FDR of 0.14. Also, the well-established effect of smoking on aryl hydrocarbon receptor signaling was retained with FDR of 0.004. Similarly, effect of smoking on UV sensitive genes as well as four other gene sets with gene regulatory effect was retained. Average methylation level analysis specific to different genomic regions We found that active smoking is significantly (Wilcoxon rank-sum test, P-value <0.05) associated with hypomethylation in shore region among active smokers [Supplementary Table S1 and Figure 3]. Although the conventional P-value threshold of 0.05 was not reached for other CpG regions, perhaps due to lack of power, smoking seems to be associated with hypomethylation in most of the regions [Supplementary Table S1]. These results suggest that effect of active smoking on average methylation in genomic regions such as shore resembles to that of aging [12]. ©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). 7 Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Table 3 Molecular Signature Database based gene sets identified with false discovery rate (FDR) ≤0.05 in the present study Significant gene sets Size (number of genes) of gene sets GSZ-score P-value FDR Gene Ontology (GO) based gene sets 1. Regulation of smoothened signaling pathway involved in dorsal ventral neural tube (biological process) 8 4.89 5.12 ×10-6 0.04 2. Thrombin-activated receptor activity (molecular function) 5 4.27 2.3 ×10-5 0.03 3. Protein domain specific binding (molecular function) 657 4.75 5.7 ×10-5 0.03 4. Olfactory receptor activity (molecular function) 321 9.96 6.09 ×10-5 0.03 5. Semaphorin receptor activity (molecular function) 9 4 6.85 ×10-5 0.03 6. AP-2 adaptor complex binding (molecular function) 7 4.22 1.15 ×10-4 0.03 7. Nuclear receptor activity (molecular function) 45 4.81 1.18 ×10-4 0.03 8. G-protein coupled receptor activity (molecular function) 725 7.21 2.28 ×10-4 0.05 Curated gene sets 9. Aryl hydrocarbon receptor signaling (canonical, reactome) 7 5.76 6.27 ×10-7 0.003 10. SHI SPARC targets up (chemical and genetic perturbations) 22 5.86 1.22 ×10-6 0.003 11. RORA activates gene expression (canonical, reactome) 17 4.99 1.68 ×10-5 0.03 12. RARRXR pathway (canonical, biocarta) 7 4.38 3.41 ×10-5 0.05 Hallmark gene sets 13. UV response down 138 4.50 0.001 0.05 Discussion To the best of our knowledge, this is the first study that implements a modern threshold free gene set analysis method to study smoking-induced alterations in methylation in genomic regions at biological pathway level. Studies of biological implications of CpG sites with altered methylation due to smoking have been limited to traditional over-representation-based pathway analysis [9,17,21,47]. We performed differential methylation analysis to first identifyCpGsitesthataresignificantlyassociatedwithsmokinghabit.Wethenimplementedourrobustgenesetanalysis method, mGSZ, on the whole methylation data at gene level to identify biological processes that might be affected due to smoking-induced altered methylation. The present study extended the current knowledge by identifying three novel smoking-associated CpG sites and altered methylation in genomic regions that play role in olfactory sensing system, immune response, cardiovascular disease, and cancer development. Furthermore, the results suggest that the global effect of smoking is similar to aging process that is known to be associated with hypomethylation in genomic regions with high frequency of CpG sites such as shores [12]. Hypomethylation in the genomic regions of active smokers is in line with findings at CpG site and gene set level. All except three CpG sites in the region of MYO1G gene are hypomethylated among active smokers. Similarly six out of seven significant gene sets based on model 2 are hypomethylated among active smokers. Previously, unidentified site cg26038589 maps to gene CCDC55. The protein product of the gene plays role in integrated regulation of gene expression in various ways and contribute to the development of cancers, viral infections, and neurological disorders [48]. CpG site cg10385208 is located in the genomic region of gene CWC25. This gene has been annotated to be overexpressed in lung and peripheral blood mononuclear cells in GeneCards—the human gene database [49]. The annotation is based on analysis of 69 proteomics datasets from Human Integrated Protein Expression Database resided within GeneCards. The gene has also been shown to play role in cognitive impairment [50]. Furthermore, the genes CCDC55 and CWC25 are linked to alternative splicing suggesting a role of DNA methylation on alternative splicing. The genes, therefore, should be interesting targets to follow up in the laboratory. The 8©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020 Bioscience Reports (2020) 40 BSR20200596 https://doi.org/10.1042/BSR20200596 Figure 2. Box plots showing differences in average methylation of all CpG sites belonging to genes of the analyzed gene sets in smokers and non-smokers X-axis represents smoking habit (active smokers and non-smokers). Y-axis represents mean methylation β-values of CpG sites belonging to genes of the analyzed gene sets. third novel site, cg09355027 has not been annotated to any gene and its biological significance is unknown. However, the CpG site is one of the numerous sites whose methylation is modulated by methyl-CpG-binding domain protein 2 (MBD2) [51]. MBD2 participates in gene silencing and has active role in shaping the cancer methylome. Site cg13898430 identified in our study has been identified previously to be associated with smoking by only one study [46]. The CpG site maps to gene RUNX3, a runt domain-containing transcription factor and a known tumor suppressor, frequently deleted or transcriptionally silenced in cancer according to the NCBI Gene database [52]. The gene is associated with numerous cancers, including lung cancer and several autoimmune diseases, including the inflammatory bowel disease. Overall, both the replicated and the discovered methylation sites map to genes that play roles in biological processes and pathways plausibly connected to smoking-induced disease. Olfactory receptor activity, one of the significant gene sets identified in the present study, is responsible for initiating cell activity in response to smell detection. Surprisingly, despite the long-established knowledge that smoking is associated with olfactory and gustatory dysfunction, relevant pathways have not been reported in previous studies. Decline in the sense of smell has been proposed to originate from several mechanisms operating in parallel, including the tobacco-induced sinonasal inflammation and squamous cell metaplasia in the airways [53]. The olfactory epithelium is a highly dynamic structure going through constant renewal through the exfoliation of the aged cells and the generation of new ones from the stem cells, thus constantly regenerating the sensory neurons [54]. The division and ©2020 The Author(s). This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY). 9 Downloaded from https://portlandpress.com/bioscirep/article-pdf/40/7/BSR20200596/887717/bsr-2020-0596.pdf by guest on 02 September 2020