Abstract

Background

DNA methylation alteration extensively associates with smoking and is a plausible link between smoking and adverse health. We examined the association between epigenome-wide DNA methylation and serum cotinine levels as a proxy of nicotine exposure and smoking quantity, assessed the role of SNPs in these associations, and evaluated molecular mediation by methylation in a sample of biochemically verified current smokers (N = 310).

Results

DNA methylation at 50 CpG sites was associated (FDR < 0.05) with cotinine levels, 17 of which are novel associations. As cotinine levels are influenced not only by nicotine intake but also by CYP2A6-mediated nicotine metabolism rate, we performed secondary analyses adjusting for genetic risk score of nicotine metabolism rate and identified five additional novel associations. We further assessed the potential role of genetic variants in the detected association between methylation and cotinine levels observing 124 cis and 3898 trans methylation quantitative trait loci (meQTLs). Nineteen of these SNPs were also associated with cotinine levels (FDR < 0.05). Further, at seven CpG sites, we observed a trend (P < 0.05) that altered DNA methylation mediates the effect of SNPs on nicotine exposure rather than a direct consequence of smoking. Finally, we performed replication of our findings in two independent cohorts of biochemically verified smokers (N = 450 and N = 79).

Conclusions

Using cotinine, a biomarker of nicotine exposure, we replicated and extended identification of novel epigenetic associations in smoking-related genes. We also demonstrated that DNA methylation in some of the identified loci is driven by the underlying genotype and may mediate the causal effect of genotype on cotinine levels.

Electronic supplementary material

Background

Smoking remains a major preventable cause of morbidity and mortality worldwide and has been shown to associate extensively with DNA methylation changes across the genome as evidenced by several epigenome-wide association studies (EWAS) [1]. Almost all EWAS so far, including a recent large meta-analysis (N = 15,907) [2], assessed the association between DNA methylation and self-reported smoking status or smoking quantity [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Self-reported smoking status and quantity, however, are prone to inaccuracies due to reporting bias (usually under reporting or recall bias) [24]. Cotinine, the primary metabolite of nicotine, is a reliable measure of nicotine exposure among current smokers [25] and provides higher statistical power compared to self-reported smoking quantity, as shown by Ware et al. [26] in a genome-wide association study (GWAS) meta-analysis of cotinine levels.

The cotinine GWAS meta-analysis [26] identified a locus on chromosome 4, within UDP glucuronosyltransferase family 2 member B10 (UGT2B10) gene, a key enzyme in nicotine and cotinine metabolism [27]. However, genetic variants in this gene did not associate with self-reported smoking quantity suggesting that cotinine may capture more than mere nicotine exposure [26]. Cotinine level is not only dependent on nicotine intake, but is also affected by the rate of formation from nicotine and the rate of cotinine metabolism to 3-hydroxycotinine, both are mediated by one highly genetically polymorphic enzyme, CYP2A6 [28]. The rate of nicotine metabolism also influences smoking behavior; for instance, fast metabolizers tend to smoke more cigarettes per day and more intensely than slow metabolizers, likely due to the longer retention of nicotine among slow metabolizers from a given intake [29, 30]. Nicotine metabolite ratio (NMR; 3-hydroxycotinine/cotinine) is a reliable proxy for CYP2A6-mediated nicotine and cotinine metabolism [29]. NMR is highly heritable with just three independent SNPs accounting for ~ 30% of variance in NMR in Finnish population [31].

In this study, we performed an EWAS to examine the association of serum cotinine levels with DNA methylation in a sample of current (daily) smokers from the Finnish Twin Cohort (N = 310). To reduce the impact of variation in cotinine levels due to CYP2A6-mediated metabolism, we utilized the genetic risk score (GRS) of NMR. As many of the genes identified in our EWAS had genetic variants that were linked to smoking-related traits previously, we further investigated the effects of genetic variants on methylation (methylation quantitative trait loci (meQTLs)) and cotinine levels. Finally, at loci where the genetic variants were associated with both methylation and cotinine levels, we examined whether altered methylation maybe a molecular mediator in the association observed between the genotype (single nucleotide polymorphism; SNP) and cotinine levels rather than being altered directly as a consequence of smoking.

Results

The overall study design and summary of our main findings is shown in Fig. 1.

Manhattan and QQ plots showing epigenome-wide associations from discovery analysis. a QQ plot showing observed versus expected − log10(P) for association at all loci. b Manhattan plot showing chromosomal locations of − log10(P) for association at each locus. All CpG sites with FDR < 0.05 are highlighted in green and the top gene for each of the highlighted loci is labeled

As cotinine levels are influenced not only by nicotine intake but also by nicotine and cotinine clearance rate (both mediated by CYP2A6), we performed a secondary analysis, wherein we included the GRS for NMR as an additional covariate in the model. In this secondary EWAS, altogether 30 CpG sites (in 20 genes) were significantly (FDR < 0.05) associated with cotinine levels. Five (in five genes) of these 30 CpG sites were novel and not genome-wide significant in the discovery EWAS (cg04776445 in MSX1, cg05335388 in ABTB2, cg10961758 in EDIL3, cg12817959 in SSTR3, and cg02306995 in SEMA5B), and the remaining 25 CpG sites (in 15 genes) overlapped with the discovery results. Half of the discovery EWAS hits (25 CpG sites) were no longer significant after controlling for nicotine and cotinine clearance rate. Manhattan and QQ plots for the secondary analysis are presented in Additional file 1: Figure S1.

Altogether, 55 unique CpG sites in 40 genes (22 of which were novel in 22 genes) were identified in our discovery and secondary EWAS analyses.

Mediation analysis

At loci where genotype was associated with both methylation (meQTL) as well as cotinine levels (19 SNPs overlap), we performed mediation analysis implemented with causal inference test (CIT) to investigate the role of DNA methylation as a molecular mediator in the observed association between genotype and cotinine levels (Additional file 7: Figure S3). CIT is a widely used approach to infer causal indirect effects of a genetic variant on an outcome using a series of statistical tests (see the “Methods” section and Additional file 7: Figure S3) evaluating conditional independencies between covariates in order to distinguish a mediated effect of the genetic variant (G) on an outcome (cotinine; C) through an intermediate (methylation, M) from a reverse cause and a common cause (pleiotropic) effect [33, 34].

Mediation analysis of 19 SNPs that associated with both cotinine levels and methylation revealed a trend (P < 0.05) for methylation at seven CpG sites (in LSM6, TNRC18, MTSS1, CASC8, PTGDR2, THSD4, and TMEM220-AS1) being a molecular mediator between the effects of SNPs (in GNG12, CNTNAP2, CYP2C18, ABTB2, and THSD4) and cotinine levels (Additional file 8: Table S5). With a cis effect observed between rs532712997 and cg26589665 in THSD4, the remaining six CpG sites showed mediating effects with multiple SNPs in trans. These results indicate that alteration of methylation at these CpG sites may be a mediator in the causal pathway between genotype and observed cotinine levels, rather than a direct consequence of nicotine exposure/smoking. These seven CpG sites (Additional file 8: Table S5) bind multiple transcription factors (TF), include DNase I hypersensitive site (DHS), and overlap with enhancers (Additional file 9: Table S6 and Additional file 10: Table S7).

Replication

To replicate our findings, we assessed the association between serum cotinine levels and methylation at the 55 CpG sites identified in our discovery and secondary EWAS in a Dutch population-based sample from the Netherlands Twin Register (NTR) as well as in an independent Finnish population cohort DIetary, Lifestyle, and Genetic determinants of Obesity and Metabolic syndrome (DILGOM). Among the 55 CpGs identified in our EWAS, 23 CpG sites were nominally associated (P < 0.05) with cotinine levels in NTR (N = 450) and 11 CpG sites in DILGOM (N = 79) with consistent direction of association, i.e., higher methylation associated with lower cotinine levels (Additional file 10: Table S7).

To replicate the trend observed in mediation analysis, we followed the same procedure as in FTC. For the 11 CpG sites where methylation levels were nominally associated (P < 0.05) with cotinine levels in the DILGOM sample (Additional file 10: Table S7), we performed genetic association and meQTL analysis. Altogether, eight SNPs were observed as meQTLs and were also nominally associated with cotinine levels (P < 0.05) in the DILGOM sample. We employed CIT to assess mediation between these eight SNPs and 11 CpG sites and observed a trend (P < 0.05) for mediation at seven CpG sites in AHRR, SLC44A1, NOTCH1, F2RL3, and CASC8 (Additional file 11: Table S8) with cg25305703 in CASC8 consistent with the mediation trend observed in FTC.

Discussion

Smoking has a major influence on methylation changes across the genome, as evidenced by numerous EWAS conducted in last few years and extensively elaborated in previous reviews [1, 35]. Cotinine is a reliable measure of recent nicotine exposure compared to self-reported smoking status or smoking quantity that are prone to reporting bias [24, 25]. In this study, we assessed the association of serum cotinine level with genome-wide DNA methylation in a sample of biochemically verified current smokers (serum cotinine > 4.85 ng/ml) from the Finnish Twin Cohort. With a half-life of 15–20 h, cotinine is a reliable indicator of recent nicotine exposure [25]. In our relatively small sample (N = 310), we replicated several previously reported findings (33 CpG sites), many of which were originally identified in studies utilizing self-reported smoking measures in much larger samples (N ~ 1000 or more) [2, 3, 4], thus showing that a biomarker with biological proximity to the phenotype provides high statistical power.

Aside from the consistently reported association with CpG sites in AHRR, F2RL3, ALPPL2, IER3, and MTSS1, which have previously been discussed extensively in the context of smoking [3, 9, 10, 22], we identified novel loci pointing to smoking-related genes such as THSD4, LSM6, CACNA2D4, and CYP2C18. An interesting novel candidate identified in our EWAS is THSD4 (thrombospondin type 1 domain containing 4). Genetic variants in this gene have been associated with several smoking-related phenotypes including nicotine dependence [36], smoking cessation success in clinical trials [37], and lung function affected by smoking [38]. Other genes highlighted in our EWAS with genetic variants associated with smoking-related phenotypes include LSM6 and CNTNAP2 associated with nicotine dependence [36] and CACNA2D4 with pack years (indicator of lifelong accumulated smoking exposure) [39]. In our data, SNPs in CNTNAP2 and THSD4 genes were significantly associated with cotinine levels, in line with these previous findings. Some of the pathways highlighted among the genes identified in our EWAS have been implicated in smoking-associated aberration in cancer pathogenesis [40, 41] and pulmonary disorders [42, 43, 44], highlighting the clinical significance of our findings.

The level of cotinine is influenced not only by intake (such as smoking pattern and smoking quantity) but also by the rate of clearance of cotinine [28]. Therefore, we utilized a GRS for NMR (a proxy for nicotine and cotinine clearance rate) to account for the impact of this in our secondary analysis and identified five additional novel loci but also observed that half of the associations from our discovery EWAS were no longer significant. A plausible explanation for this observation could be that adding the GRS adjusts for the genetic variation of nicotine metabolism rate (normal vs faster metabolizers) among the individuals analyzed. Associations which disappear when the NMR GRS is added to the model may be related to the impact of nicotine metabolism rate on nicotine intake or nicotine/cotinine metabolism. Among such associations (no longer genome-wide significant after adjusting for the NMR GRS) are CpG sites in genes AHRR and CYP2C18, members of the metabolic pathways involved in nicotine and bupropion (prescription medication for smoking cessation [45]) degradation and xenobiotic metabolism [46]. CYP2C18, a member of the cytochrome P450 (CYP) family, and AHRR, a regulator of CYP genes, are both involved in xenobiotic (which includes components of cigarette smoke) metabolism. Altered expression of CYP2C18 in the context of smoking has been reported previously [47]; however, the role of methylation as a plausible regulator in such observed alterations warrants further experimental investigation. It should be noted that although just three SNPs used in the NMR GRS account for a relatively high proportion (~ 30%) of variance in NMR, a large portion of variance remains unaccounted for and this is a major limitation of our study. The three independent SNPs constituting the GRS were identified in the Finnish population sample and may not capture the same effects across other populations (as we observed with the NTR replication sample).

As some of the CpGs identified in our EWAS were annotated to genes that have previously been linked to smoking-related phenotypes (36–39), we examined the role of genetic variants in the genes identified in our EWAS. It is noteworthy that a plethora of associations between the methylation levels of highlighted CpG sites and SNPs were observed in both cis and trans. We can speculate a cross-genome interplay between nicotine exposure-associated DNA methylation and genetic variation in the highlighted genes [48, 49]. Some of the cis meQTLs also act in trans with other CpG sites. Most of the CpG sites bind multiple TFs, overlapping enhancers or DHS sites hinting at global-directed networks at play [48, 50]; however, these may be independent effects occurring simultaneously. It should also be noted that we tested meQTLs only among genes highlighted in our analysis potentially missing other cis and trans meQTLs in the whole genome. Interestingly, for the 19 meQTLs that were also associated with cotinine levels, we observed a trend of molecular mediation by methylation at seven CpG sites. These results suggest that alterations in methylation of these nicotine exposure-associated CpG sites might not always be a consequence of smoking as commonly suggested, but rather methylation at such CpG sites may act as a causal mediator (molecular mechanism) in regulating the effect of genetic variation (in genes that modulate nicotine intake and/or nicotine and cotinine metabolism rate) on cotinine levels or smoking behavior. We observed support for these findings in an independent replication sample (DILGOM). These findings are particularly interesting as these SNPs do not overlap SNPs in the probes. We observed cis mediation between CpG sites and SNPs in THSD4 (FTC) and AHRR (DILGOM) both involved in smoking [36, 37, 38, 51], while trans mediation in CpG sites that bind multiple TFs or overlap DHS sites and enhancers providing plausible action mechanism for molecular mediation. These findings are crucial and provide caution to interpretation of smoking-related EWAS but require further experimental evidence.

Cotinine, even though a reliable measure of nicotine exposure, is not specific to tobacco smoking. Other potential sources of nicotine, such as smokeless tobacco (e.g., snus), e-cigarettes, and/or nicotine replacement therapy, could also contribute to cotinine levels; however, the likelihood of such alternative sources in our sample was low (see the “Methods” section). There have been only a few EWAS so far using cotinine as a continuous phenotype [51, 52], both of which were performed in samples including smokers and non-smokers. When considering replication of prior findings, it should be noted that all the previous EWAS we used to assess the novelty of our findings (see “Methods”) were conducted in samples including both smokers and non-smokers. We cannot rule out that our novel findings might be specific to current or heavy smoking as we limited our analyses to self-reported current smokers with elevated cotinine levels only. Having included only current smokers, the range of variability is limited (in contrast to studies including both smokers and non-smokers) resulting in smaller effect sizes observed. Similar results (small effects) were noted by Zhang et al. [52] and Su et al. [22] when restricting analysis to current smokers only. Although we observed a great overlap with previous findings, non-overlapping associations could be due to differences in study design and sample (smokers and non-smokers in previous publications versus current smokers only in our study), population specificity of methylation levels [12], and precision of cotinine measurement technique (mass spectrometry versus immunoassay) [53]. The lack of replication of the novel associations identified in our EWAS in the NTR and DILGOM samples may be due to several aforementioned factors including population specificity, differences in precision of cotinine measurement (mass spectrometry (FTC and DILGOM) and immunoassay (NTR)) and small sample size (DILGOM N = 79), or small and inconsistent effects of novel associations. Even though we replicated several prior associations between smoking and methylation, the novel associations identified in our study, while biologically meaningful, would require replication in other larger samples.

Conclusions

In conclusion, we examined the epigenetic signature of nicotine exposure in a sample of biochemically verified current smokers and identified several novel loci, including genes involved in nicotine degradation and metabolism. Our study is the first cotinine EWAS in which the rate of nicotine clearance is accounted for, allowing us to discover smoking-pertinent novel associations which may be specific to regular/heavy smoking. We further expose genetic influences on methylation and provide suggestive evidence for the role of methylation as a molecular mediator between genetic variation and cotinine levels, as opposed to a direct consequence of nicotine exposure in some of the genes.

Methods

Study sample

The study sample of cotinine verified current smokers comes from the Finnish Twin Cohort (FTC) [54], a population-based longitudinal study designed to examine genetic and environmental determinants of health-related behaviors. A total of 310 individuals (51 full monozygotic pairs, 44 full dizygotic pairs, and 120 twin individuals without a co-twin) were included in our discovery analysis. Genome-wide DNA methylation was assessed with Infinium HumanMethylation450 BeadChip in peripheral blood using standard protocols [55]. The average age in the sample was 29.5 years (SD 14.2, range 21–69), and 52% were females. Genotype data imputed to 1000Genomes phase 3 (processing and imputation of the genotype data has been described in detail previously [31]) was available for 304 of these individuals, which were consequently included in the secondary analysis.

Phenotype data

Current smokers were selected based on the threshold of serum cotinine above 4.85 ng/ml, as suggested by Benowitz et al. [56]. Selecting current smokers with a threshold of cotinine > 4.85 ng/ml instead of 10 ng/ml was done to maximize sample size available for discovery analysis, and sensitivity analysis (results not shown) was performed to ensure that using a higher threshold did not affect the results and subsequent inference. Cotinine was measured from frozen serum samples using high-precision LC-MS/MS at University of Toronto, Canada (N = 242) [57] and at the Metabolomics Unit, Institute for Molecular Medicine Finland (FIMM), University of Helsinki, Finland (N = 68). The protocols and instrument parameters have been described previously [31, 58, 59]. Average cotinine in our discovery sample (N = 310) was 192.7 ng/ml (median = 178.8, SD = 148.4, range = 5.1–820.5). In our discovery sample, two individuals (with cotinine levels above 100 ng/ml) reported using nicotine replacement therapy and none declared use of smokeless tobacco. E-cigarette use was non-existent in Finland at the time of data collection.

Replication samples

DILGOM

To replicate our findings, we utilized a Finnish population cohort DILGOM (DIetary, Lifestyle, and Genetic determinants of Obesity and Metabolic syndrome) [60, 61]. Briefly, DILGOM originates from the population-based national FINRISK 2007 study and includes a total of 631 unrelated Finnish individuals aged 25–74 years from the Helsinki area. A total of N = 79 individuals with methylation data from the peripheral blood, genome-wide genotype data imputed to 1000Genome phase I, and serum cotinine > 4.85 ng/ml (current smokers) were analyzed. The average age of the sample was 51 years (SD 13.5, range 25–72), and 44% were females. Cotinine was measured with gas chromatograph mass spectrometry at the Laboratory of Analytical Biochemistry at the Institute of Health and Welfare, Helsinki, Finland, as described earlier [57]. All 79 individuals were self-reported current smokers. Average cotinine levels were 134.2 ng/ml (median = 129.6, SD = 98.0, range 5.0–408.8).

DNA methylation data processing and analysis

DNA methylation data was preprocessed and normalized using the pipeline suggested by Lehne et al. [67] implemented using R packages “minfi” [68] and “limma” [69]. We modified the pipeline to accommodate the relatedness in our sample. Altogether the quality control involved exclusion of probes with detection P > 1 × 10− 16, bead count < 3 (estimated with “wateRmelon” R package [70]) and call rate < 98%. Samples with sample call rate < 98% and sex mismatch were further excluded. After quality control, 418,302 probes and N = 310 samples remained (FTC data). Quantile normalization stratified on probe type, color channel, and probe subtypes was performed to obtain a methylation beta matrix. The first 30 principal components (PCs) from the control probes were calculated to adjust for technical bias. White blood cell subtypes (CD8T, CD4T, Natural Killer cells, B cell, monocyte, and granulocyte) were estimated using “FlowSorted.Blood.450 k” [71] within the “minfi” R package based on a modified version of the Houseman algorithm [72] and included in the regression model. Intermediary residuals were estimated by regressing the methylation beta values on age, sex, body mass index (BMI), 30 control probe PCs, and white blood cell types. PCs from these intermediary residuals were further included in the final regression model to adjust for any unaccounted global covariation. To analyze the association between methylation (dependent variable) at each CpG site and serum cotinine levels (explanatory variable), we used linear mixed effects model with “lmer” function implemented in “lme4” R package [73]. Age, sex, BMI, batch variable (site of cotinine measurement), 30 control probe PCs, and 10 intermediary residual PCs were included as covariates (fixed effects), while family identifier and zygosity were included as random effects to account for the relatedness in our sample. To avoid spurious associations, we further applied ad hoc exclusion of probes reported by Zhou et al. [74] to affect beta methylation at the CpG (probes with non-unique mapping, inconsistent extension base, SNP in the extension base, and overlap with any known SNPs with global minor allele frequency (MAF) and MAF in a Finnish population > 1%).

We applied false discovery rate (FDR) adjustment to the p values for multiple testing correction and considered FDR-adjusted P < 0.05 as statistically significant. Manhattan and QQ plots were produced using the R package “qqman” [75]. Lambda value estimates for the QQ plot were calculated using the “estlambda” function in R package “GenABEL” [76]. For replication of our results, analysis pipeline was identical for the replication sample NTR. For the DILGOM sample, analysis was identical except that instead of a mixed effects model, a linear model was used as all individuals were unrelated.

Secondary analysis

To account for the individual rate of nicotine metabolism, we included the GRS for NMR in the regression model (described above) as an additional covariate. We extracted the genotypes for the three independent SNPs (rs56113850, rs113288603, and rs12461964) and calculated the weighted mean of minor allele counts to calculate the GRS, as described previously [31]. One of the top variants (indel esv2663194) identified in the NMR GWAS meta-analysis [31] was not available in our samples (discovery and replication). GRS calculation and analysis pipeline was identical for the replication samples.

Literature search

We used the same literature search strategy as Joehanes et al. [2]. Briefly, PubMed literature database was queried using medical subject heading (MeSH) terms of DNA methylation and smoking. To limit duplication of efforts, we used the supplementary data on literature search used by Joehanes et al. which included articles until 2015 and applied additional filter on year (starting 2015) along with a filter for species (human). We reviewed the abstracts to determine if the studies were as follows: (1) performed in healthy adult human populations (also excluding maternal smoking), (2) sample type analyzed was peripheral blood, (3) studies assessed genome-wide methylation only, and (4) public reporting of P values and gene annotations for the CpGs identified was available. A total of 24 publications met all the inclusion criteria and are listed in the Additional file 12: Table S9. A compiled list of CpGs (19,208 CpG sites) and associated gene names (8582 genes) from all 24 studies were utilized to assess novelty of our findings.

Annotation and pathway analysis

All statistically significant CpG sites were annotated using the data aggregated by Zhou et al. [74]. For the probes that did not have a gene name listed, the name of the nearest gene was fetched from Ensembl [77]. This gene list was analyzed with the Ingenuity Pathway Analysis software (IPA; Ingenuity® Systems, https://www.qiagenbioinformatics.com/) to identify gene networks at play. We applied multiple testing correction and an enrichment for genes functioning in signaling and metabolic pathways [78] was considered significant when FDR < 0.05. A total of 38 of the 40 highlighted genes were recognized by IPA. We also performed pathway analysis (Gene Ontology) to assess enrichment of biological processes among the top 55 CpG sites while taking into account the non-uniform CpG probe distribution per gene [79] on the 450 k array with R package “missMethyl” [80] and corrected for multiple testing and considered FDR < 0.05 significant.

Association of cotinine with genetic variants

As the genes identified in our EWAS harbor genetic variants associated with smoking behavior phenotypes, such as nicotine dependence [36], we examined whether SNPs in the genes identified in our EWAS analysis associate with cotinine levels. We performed association analysis using mixed effects model implemented as “lmer” function in “lme4” R package. Cotinine levels were treated as dependent variable, while genotype dosage (coded as 0, 1 or 2 for copies of effect allele), age, sex, and BMI were included as fixed effects in the model. Family and zygosity were included as random effects to account for the relatedness in our sample. Only polymorphic SNPs (N = 46,780) were tested for association with cotinine levels. Analysis was otherwise identical for the DILGOM sample except that linear model was used instead of mixed effects model, as all individuals were unrelated. To estimate proportion of variance explained by individual SNPs, we subtracted the R2 of the model including only covariates (age, sex, BMI) from the R2 of the model also including the SNP.

Methylation quantitative trait loci analysis

To examine the association between methylation and SNPs within the 40 genes (with 50 kb flanking regions; gene boundaries are available in Additional file 13: Table S10), we used R package “MatrixEQTL” [81]. Only polymorphic SNPs (N = 46,780) were tested using linear model setting with a cis distance of 2.5 Mb. The longest gene was approximately 2.4 Mb (Additional file 13: Table S10); thus, a cis distance of 2.5 Mb was chosen to ensure that all possible combinations of SNPs, and the CpG within a gene are tested. All association tests between remaining SNP and CpGs were considered trans. Genotypes were coded as copies of effect allele (0, 1, or 2) and methylation data was extracted for the 55 probes from the CPACOR normalized data set used in the EWAS. We included, age, sex, BMI, white blood cell counts, control probe PCs, and cotinine as covariates in the model. In addition, owing to the relatedness in our sample, we included a covariance matrix based on CpG sites being tested in the model. Results of meQTL analysis were visualized using R package “RCircos” [82].

Mediation analysis

For the SNPs that were significantly associated with both cotinine levels and methylation, we investigated the potential of methylation (M) being a mediator between the effects of genetic variant (G) on outcome cotinine levels (C) using causal inference test (CIT) [83]. In CIT, four conditions are tested to assess consistent causal mediation as described by Millstein et al. [83]. Conditions for CIT are (1) G is associated with C, (2) G is associated with M conditional on C, (3) M is associated with C conditional on G, and (4) G is independent of C conditional on M. We implemented CIT analysis with R package “cit” and function “cit.cp” with 50 permutations for conditional analysis and included, age, sex, BMI, family, and zygosity as covariates. We report the overall p value (omnibus p value) based on collective conditional analysis and considered P < 0.05 as nominally significant.

Notes

Acknowledgements

Funding

Data collection of the Finnish Twin Cohort samples has been supported by the Academy of Finland Center of Excellence in Complex Disease Genetics (grants 213506, 129680), the Academy of Finland (grants 265240 and 263278 to JK, and 297908 to MO), Sigrid Juselius Foundation (to JK and MO), University of Helsinki Research Funds (to MO), and Global Research Award for Nicotine Dependence, Pfizer Inc. and Finnish Foundation for Cardiovascular Research (to JK), Canada Research Chair in Pharmacogenomics and Canadian Institutes of Health Research (CIHR) grant FDN-154294 (to RFT). Data from the Netherlands Twin Register were funded by the BBRMI-NL-financed BIOS Consortium (NWO 184.021.007), and Genetics of Mental Illness, a lifespan approach to the genetics of childhood and adult neuropsychiatric disorders and comorbid conditions (ERC-230374).

Availability of data and materials

Due to the consent given by study participants, FTC data cannot be made publicly available. Data are available through the Institute for Molecular Medicine Finland (FIMM) Data Access Committee (DAC) for authorized researchers who have IRB/ethics approval and an institutionally approved study plan. For more details, please contact the FIMM DAC (fimm-dac@helsinki.fi). Data from DILGOM cohort can be accessed by application to the THL Biobank (thl.fi/web/thl-biobank). The HumanMethylation450 BeadChip data from the NTR are available as part of the Biobank-based Integrative Omics Studies (BIOS) Consortium in the European Genome-phenome Archive (EGA), under the accession code EGAD00010000887.

Authors’ contributions

The study design was conceived and evolved by MO, AL, JK, and RG. RG analyzed and interpreted the data and wrote the paper. JVD, AA, and YF analyzed the replication data. JK, RFT, VV, and DIB generated the phenotype data. JK, AL, and DIB acquired the genotype data. JK, MO, and DIB acquired the DNA methylation data. AL, MO, JK, RFT, TK, and JVD critically revised the manuscript for important intellectual content. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Written informed consent was obtained from all subjects who gave DNA samples in accordance to the current edition of the Declaration of Helsinki and the collection of blood samples followed the recommendations given in the Declaration of Helsinki and its amendments. Data collection for the Finnish Twin Cohort has been approved by the Hospital District of Helsinki and Uusimaa, the ethics committee for epidemiology and public health (HUS-113-E3–01, HUS-346-E0–05, HUS 136/E3/01, HUS 270/13/03/01/2008, HUS 154/13/03/00/11). DILGOM participants provided written informed consent and the Ethics committee of Helsinki University Central Hospital has provided the approval (#46 20.2.2007, #90 3.4.2007, #332 13.03.2013). The protocol was designed and performed according to the principles of the Helsinki Declaration and was approved by the coordinating ethics committee of the hospital district of Helsinki and Uusimaa, Finland. Informed consent was obtained from all participants of NTR study. The study was approved by the Central Ethics Committee on Research Involving Human Subjects of the VU University Medical Centre, Amsterdam, an Institutional Review Board certified by the U.S. Office of Human Research Protections (IRB number IRB00002991 under Federal-wide Assurance FWA00017598; IRB/institute codes, NTR 03–180).

Consent for publication

Not applicable.

Competing interests

JK and TK have provided consultation to Pfizer on nicotine dependence and its treatment. RFT has consulted for Apotex, Quinn Emmanuel, and Ethismos and is a member of several scientific advisory boards (e.g., Canadian Centre for Substance Abuse, Quitta, Health Canada (Vaping), and Brain Canada). The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary material

Additional file 1:Figure S1. Manhattan and QQ plots showing epigenome-wide association results from secondary analysis when accounting for the rate of nicotine metabolism using a GRS. (A) QQ plot showing observed versus expected − log10(P) for association at all loci. (B) Manhattan plot showing chromosomal locations of − log10(P) for association at each locus. All CpG sites with FDR < 0.05 are highlighted in green and the top gene for each of the highlighted loci is labeled. (PDF 171 kb)

Copyright information

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.