A blood- and brain-based EWAS of smoking

Cigarette smoking remains a leading cause of preventable death and disease, accounting for approximately 8 million global deaths annually¹. It is a major risk factor of more than 50 diseases including cardiovascular disease, lung cancer and dementia². As smoking history is often used in clinical risk stratification assessments, enhancing the accuracy of cumulative tobacco consumption measurements has the potential to improve prevention and treatment of smoking-related diseases.

Traditionally, tobacco use has often been quantified using self-report questionnaire-response data, such as pack years or indication of current smoking status (i.e. current, former, never smoker), which are prone to recall bias and do not account for passive smoke exposure³. A more objective approach to assess smoking is by measuring the concentration of tobacco-related chemicals. However, the commonly used nicotine biomarker, cotinine, has an average half-life of 15–20 hours⁴. Consequently, the concentration of serum cotinine does not inform about time since cessation in recent quitters and it cannot help to distinguish former smokers from never smokers. This limitation is especially pertinent when estimating the risk of diseases that take many years to develop, such as cardiovascular disease.

Blood-based DNA methylation patterns show great promise as a long-term biomarker of smoking⁵. DNA methylation (DNAm) is a cell- and tissue-specific epigenetic modification of DNA molecules that does not change the DNA sequence itself. It involves the addition of a methyl group to cytosine residues and occurs predominantly at cytosine-phosphate-guanine dinucleotides, also known as CpG sites. CpG methylation levels reflect not only smoking status but also cumulative tobacco exposure and can be informative of time since quitting in former smokers^6,7,8. In the majority of CpG sites, smoking-related DNAm changes are dose-dependent and reversible after cessation⁹.

In previous studies, blood-based DNAm biomarkers of smoking almost perfectly discriminated smokers from never smokers^7,10. However, their ability to differentiate former smokers from never smokers was relatively modest. Additionally, most of these studies relied on whole-blood DNAm assessments using arrays, which only measure a pre-selected subset of CpG sites present in the epigenome. For example, the current largest epigenome wide association study (EWAS) of smoking in adults (n = 15,907) was a meta-analysis conducted using the Illumina 450k BeadChip array (n ~ 450,000 CpG sites)⁶. The largest EWAS of smoking, in which DNAm was measured with Illumina EPIC array (n ~ 850,000 CpG sites, approximately 5% of all sites on the epigenome), considered blood samples of 15,014 individuals¹¹.

In this work, we update the existing biomarkers of smoking by analysing whole-blood DNAm levels measured with Illumina EPIC array in 17,865 individuals. For a subset of 46 individuals, we implement a high-resolution approach ( ~ 4 million CpG sites, TWIST human methylome panel and ~21 million sites, Oxford Nanopore Sequencing) aimed at measuring methylation levels at CpG sites which are currently absent from arrays. Furthermore, we develop a smoking biomarker (mCigarette) using the EPIC dataset and investigate its associations with self-reported smoking in two external cohorts (Lothian Birth Cohort – LBC1936 and The Avon Longitudinal Study of Parents and Children - ALSPAC). We also investigate variations in methylation patterns in both blood and brain, with DNA methylation measured across five post-mortem brain regions in 14 individuals. Finally, we compare the epigenetic proxy of smoking to self-reported smoking by considering the genetic loci associated with these phenotypes in genome-wide association studies (GWASs).

EWAS of smoking

First, we ran a Bayesian EWAS of smoking (Fig. 1). There were 17,865 Generation Scotland (GS) volunteers with a measure of pack years of smoking and Illumina EPIC DNAm data. The average age of the sample was 47.6 years (standard deviation [SD] = 14.9), and 59.1% of the participants were female (Supplementary Data 1). In the BayesR EWAS, DNAm explained 50.0% (95% Credible Interval 46.0 – 53.9) of the variance in the pack years phenotype.

The analysis was conducted using data from three British cohorts: Generation Scotland (GS), the Lothian Birth Cohort of 1936 (LBC1936) and the Avon Longitudinal Study of Parents and Children (ALSPAC). While blood-based DNA methylation (DNAm) data were available in all three cohorts, LBC1936 also contained information about the DNAm levels in post-mortem brain tissue, covering five brain regions from 14 individuals. A Bayesian Epigenome-Wide Association Study (EWAS) of smoking was performed in GS ( ~ 850k sites, Illumina EPIC array). For 23 pairs of age- and sex-matched smokers and non-smokers from GS, a high-resolution methylation measurement approach was implemented ( ~ 4 million sites, TWIST human methylome panel and ~21 million sites, Oxford Nanopore Technologies sequencing), followed by an EWAS analysis. An epigenetic biomarker of smoking, mCigarette, was developed in GS and tested as a predictor of self-reported smoking in LBC1936. The association between mCigarette and self-reported smoking was replicated in multiple age groups present in ALSPAC. Next, EWASs of smoking were run across five brain regions for 14 individuals using EPIC DNAm from LBC1936. Finally, Genome-Wide Association Studies (GWAS) were run in GS to compare genetic signal of self-reported and epigenetic smoking (GrimAge DNAm pack years score). EPIC – Illumina EPIC array, TWIST – TWIST Biosciences Human Methylome Panel, ONT – Oxford Nanopore Technologies Sequencing. Created in BioRender. Marioni, R. (2024) https://BioRender.com/h44z126.

Full size image

Forty-two independent CpGs were associated with smoking at posterior inclusion probability (PIP) > 80%, with 26 of these associations reaching a PIP > 95% (Supplementary Data 2). Among the associations with PIP > 80%, 33 had previously been reported in the EWAS catalogue¹² at P < 1×10⁻⁴, with 30 of these reaching P < 1×10⁻⁷. This catalogue is a resource that curates findings from published EWAS studies.

Associations unreported in the EWAS catalogue as being associated with smoking included nine sites at PIP > 80% and two CpGs with PIP > 95%. The former group included intergenic CpGs linked to neurodevelopment and addiction, such as cg22454588 (annotated to SCAMP5), cg27110277 (FGF20), and cg19404444 (SKI). The latter, high confidence associations, included cg02517189 (GRIK5) and cg00562553 (HOXA4).

High resolution EWASs of smoking

Next, we extended this analysis by running high resolution EWASs of smoking on a subset of 23 pairs (n = 46) of current vs never smokers with Illumina EPIC array ( ~ 850k CpG sites), TWIST human methylation panel ( ~ 4 million CpG sites), and Oxford Nanopore Technologies (ONT) sequencing data ( ~ 21 million CpG sites). At P < 3.6×10⁻⁸ (significance threshold set as per Saffari et al.¹³), the EPIC-based analysis revealed 15 CpG sites associated with smoking status (EWAS inflation factor, λ = 0.94), while the TWIST-based and ONT-based analyses identified 33 (λ = 1.60) and 9 (λ = 1.11) associations, respectively. At a less stringent threshold (P < 1×10⁻⁵), these counts increased to 42, 102, and 63 for the EPIC-, TWIST-, and ONT-based analyses, respectively. The overlap between the sites identified by these technologies is detailed in Supplementary Data 3. Figure 2 and Supplementary Fig. 1 (comparison of beta estimates) display the results obtained from the TWIST, ONT and EPIC EWAS.

Analyses were performed using DNA methylation data obtained using the Illumina EPIC array ( ~ 850k CpG sites), the TWIST human methylation panel (4 million CpG sites, targeted short read sequencing) and Oxford Nanopore Technologies sequencing (21 million CpG sites, long read sequencing). The X-axis represents chromosomes 1–22, while the Y-axis shows −log10(P-values). The top horizontal line marks genome-wide significant associations (P < 3.6 × 10⁻⁸, red dots), based on the multiple testing threshold estimated by Saffari et al.¹³. The bottom horizontal line denotes the suggestive significance threshold (P < 1 × 10⁻⁵, yellow dots). Dotted vertical lines highlight loci associated with smoking status at P < 3.6 × 10⁻⁸. All statistical tests were two-sided. EPIC – Illumina EPIC array, TWIST – TWIST Biosciences Human Methylome Panel, ONT – Oxford Nanopore Technologies Sequencing.

Full size image

Among the 33 associations identified as significant in the TWIST EWAS at P < 3.6x10⁻⁸, two had been previously reported in the EWAS catalog (based on DNAm profiled with array technologies). These included AHRR (chr5-373263-373264, beta = −0.35, P = 1.2x10⁻¹⁰) and an intergenic locus found on chromosome 2 (chr2-232419951-232419952, beta = −0.24, P = 3.5×10⁻⁸). The remaining 31 loci significant at this threshold included sites annotated to F2RL3 (chr19-16889741-16889742, beta = -0.31, P = 1.3x10⁻¹¹) and USP42 (chr7-6126706-6126707, beta=0.05, P = 1.8×10⁻⁸). At a less stringent threshold of P < 1×10⁻⁵, 98 uncatalogued sites were identified. They included SST (chr3-187670342-187670343, beta = −0.09, P = 2.2x10⁻⁷) and TSPAN5 (chr4-98472405-98472406, beta = -0.06, P = 6.8×10⁻⁶). Further details are provided in Supplementary Data 4.

In the ONT EWAS, 9 sites were significant at P < 3.6×10⁻⁸, of which only one had been previously listed in the EWAS catalog: a site mapping to AHRR (chr5-373263-373264, beta = -0.47, P = 2.4x10⁻¹²). The remaining eight included additional loci within the AHRR region, loci from an intergenic region on chromosome 2 (e.g., chr2-232420079-232420080, beta = -0.37, P = 3.4x10⁻¹⁰) and a site annotated to CNTNAP2 (chr7-147245588-147245589, beta=0.37, P = 1.1×10⁻⁸). At P < 1×10⁻⁵, the ONT analysis revealed 62 uncatalogued loci, such as SEPTIN9 (chr17-77351321-77351322, beta = -0.27, P = 8.4x10⁻⁶) and TERF2 (chr16-69398923-69398924, beta=0.12, P = 8.7x10⁻⁶). Complete results are available in Supplementary Data 5.

A gene set enrichment analysis of genes mapped to the 102 CpGs with P < 1×10⁻⁵ identified in the TWIST EWAS revealed 13 enriched gene sets (FDR p < 0.05; see Supplementary Data 6). These included tissue degradation through altered extracellular matrix dynamics and chronic inflammation stemming from dysregulated ATP release and immune cell recruitment. Many of the enriched pathways were driven by the presence of collagen genes, such as COL4A4 and COL4A3. In contrast, enrichment analysis based on the significant CpGs from the ONT EWAS did not identify any significantly enriched gene sets.

DNAm biomarker of cigarette consumption - mCigarette

A DNAm biomarker of cigarette consumption, mCigarette, was then developed. The biomarker was trained in GS (n = 17,865) using elastic net regression with 10-fold cross-validation. Prior to training, CpG sites were prefiltered to those associated with tobacco use at FDR < 0.05 (n_CpG=18,760) in the previous largest EWAS of smoking⁶. This meta-analysis did not include GS and was conducted using Illumina 450k BeadChip array. Following 10-fold cross validation, an optimal lambda value that minimised the mean prediction error was selected (lambda=0.012577) and fed into an elastic net model of smoking pack years. As a result, non-zero coefficients were assigned to 1,255 CpG sites (Supplementary Data 7).

The biomarker was tested in the external LBC1936 study (n = 882, mean age 69.6 years, SD = 0.8). Assessing the predictive performance using Area Under the Curve (AUC) (Fig. 3) revealed very good (AUC = 0.85) to near perfect (AUC = 0.98) ability to distinguish between current (n = 101), former (n = 368) and never smokers (n = 413). Similar results were obtained when Area Under the Precision-Recall Curve was used as a performance metric (PRAUC _{current vs never} = 0.96, PRAUC _{former vs never} = 0.85, PRAUC current vs former = 0.71).

The Areas Under the Curves (AUCs) represent the ability of the model to distinguish between the following smoking categories: blue for current smokers vs. never smokers, orange for current smokers vs. former smokers, and green for former smokers vs. never smokers. Source data are provided as a Source Data file.

Full size image

The predictive performance of mCigarette was benchmarked against four previously developed epigenetic scores for smoking, as well as a single-site biomarker based on AHRR methylation at cg05575921, in wave 1 of LBC1936 (n = 882). The results of this comparison can be found in Table 1. mCigarette yielded improved incremental R² estimates compared to EpiSmokEr, the score developed by McCartney et al.¹⁰ and GrimAge DNAm pack years.

Subsequently, the weights used to construct the studied scores (apart from the GrimAge DNAm pack years, as weights are not publicly available) were applied to methylation data in the ALSPAC cohort. Within this dataset, no single methylation score consistently outperformed others in distinguishing between current, former, and never smokers across studied age groups. While in young adults (mean age 18 and 24 years old) EpiSmokEr and the score constructed by McCartney et al.¹⁰ achieved the highest AUCs (median AUC = 0.720 [IQR:0.642-0.763]), mCigarette and EpiSmokEr showed excellent performance in older adults (mean age 29 and 50 years old, median AUC = 0.890 [IQR:0.771-0.931]). Full results of the replication study are available in Supplementary Data 8 and are visualised in Supplementary Fig. 2.

Tissue specificity

To assess whether the findings translated across different tissue types, the association between tobacco use and DNAm levels across five brain regions were explored using post-mortem samples from LBC1936 (n = 14, n_hippocampus = 13). At significance threshold P < 1x10⁻⁵, five loci in the hippocampus (BA35), one locus in dorsolateral prefrontal cortex (BA46), four loci in primary visual cortex (BA17), nine loci in anterior cingulate cortex (BA24) and three loci in ventral/lateral inferior temporal cortex (BA20/21) were associated with smoking status (Supplementary Data 9). There was no overlap between the significant loci across the studied brain regions.

Some loci demonstrated nearly perfect discrimination of smoking status in blood and brain; however, these loci did not overlap (Fig. 4 and Supplementary Fig. 3). For instance, the methylation status at cg05575921, annotated to the AHRR gene, is a well-established marker of smoking status in whole-blood DNAm. However, this marker did not discriminate smoking category in hippocampal DNAm. On the other hand, cg26381592, annotated to the PMS1 gene, did not effectively distinguish smokers in blood samples, but it exhibited a strong correlation with smoking status in hippocampus samples.

CpG methylation levels in blood (Lothian Birth Cohort 1936 baseline, n_{blood_DNAm} = 882) and brain across different self-reported smoking categories (Lothian Birth Cohort 1936, n_{hippocampus_DNAm} = 13). Box plots are defined as follows: the centre line represents the median (50th percentile). The box bounds indicate the interquartile range (IQR; 25th to 75th percentile). Whiskers extend to the smallest and largest values within 1.5×IQR. Source data are provided as a Source Data file.

Full size image

Subsequently, we tested the performance of a blood-derived DNAm biomarker of smoking (mCigarette) in the brain tissue. When applied to brain DNAm data, mCigarette did not distinguish between smoking categories in any of the studied brain regions (Supplementary Fig. 4).

GWAS of self-reported and epigenetic smoking

Finally, GWASs of self-reported and epigenetic smoking were conducted. To avoid overfitting, we used the GrimAge DNAm pack years estimator (trained externally in 1731 individuals from the Framingham Heart Cohort study¹⁴) instead of mCigarette, which was trained in GS.

In 17,105 GS individuals, there was a single nucleotide polymorphism (SNP)-based heritability of 27.3% (P = 7.0x10⁻⁴⁶) for self-reported pack years of smoking compared to 41.0% (P = 8.2×10⁻⁹⁸) for GrimAge DNAm pack years (Fig. 5). The genomic inflation factors (λ) were 1.06 and 1.10 for pack years and GrimAge DNAm pack years, respectively. At the genome-wide significance level of P < 5.0x10⁻⁸, only one locus (rs117836409 annotated to GDPD1) was associated with self-reported pack years, in contrast to 39 SNPs (three lead SNPs at two genomic risk loci) associated with GrimAge DNAm pack years (detailed in Supplementary Data 10 and 11). Of the three lead SNPs associated with GrimAge DNAm pack years, two (rs1800440 and rs6495309, annotated to CYP1B1 and CHRNA3 - CHRNB4, respectively) have been previously documented in the GWAS catalogue, which aggregates data from published GWAS studies¹⁵ (Supplementary Data 12). They have been associated with carcinogenesis and nicotine dependence, respectively^16,17. The SNP not previously annotated in the GWAS Catalog (rs114342890) maps to RMDN2:RMDN2-AS1, a long non-coding RNA previously studied in relation to eosinophil counts and melanoma^18,19. According to GeneHancer, an online database of enhancers, promoters and their inferred targets, all 27 regulatory elements which target RMDN2:RMDN2-AS1 also regulate the expression of CYP1B1²⁰. The beta coefficients of the lead loci identified in the GrimAge DNAm pack years GWAS and the smoking pack years GWAS are compared in Supplementary Fig. 5.

Studied phenotypes included: pack years (Generation Scotland, n = 17,105), pack years (largest meta-analysis to date²¹, n = 131,892), GrimAge-based DNAm estimator of smoking pack years (n = 17,105). The X-axis represents genomic positions, while the Y-axis shows −log10(P-values). The top horizontal line marks genome-wide significant associations (P < 5×10⁻⁸, red dots), based on the Bonferroni threshold for multiple testing. The bottom horizontal line denotes the suggestive significance threshold (P < 1x10⁻⁵, yellow dots). Dotted vertical lines mark genomic positions of significant associations at P < 5×10⁻⁸. All tests were two sided.

Full size image

The three lead SNPs have been characterised as methylation quantitative trait loci (mQTLs). They all act in cis on 43 CpGs (Supplementary Data 13). Only one of these CpGs (cg06264984) is featured in mCigarette, while none are included in EpiSmokEr – the list of CpGs included for the GrimAge DNAm pack years is not publicly available.

Next, the GrimAge DNAm pack years GWAS results were compared to previously published GWAS studies of tobacco use (Supplementary Data 14). At a significance level of P < 5×10⁻⁸, seven SNPs annotated to CHRNA3 and CHRNA5 aligned with the findings of the largest pack years GWAS to date (n = 131,892)²¹. Thirty-seven SNPs, mapping to CHRNA3, CHRNA5, and CHRNB4, overlapped at P < 5x10⁻⁸ with the results of a related phenotype, cigarettes per day (n = 618,489)²².

GrimAge DNAm pack years and self-reported pack years were moderately correlated (Spearman’s r = 0.65). The genetic correlation (r_g) between GrimAge DNAm pack years and pack years from Erzurumluoglu et al.²¹ was 0.62 (SE = 0.12, P = 4.4x10⁻⁷), with an LD score regression intercept of 0.99 (SE = 0.01). The r_g between self-reported pack years in the 17,105 GS and pack years from Erzurumluoglu et al.²¹ (n = 131,892) was 0.67 (SE = 0.19, P = 5.0x10⁻⁴), with an LD score regression intercept of 0.99 (SE = 0.01). Additional information on the genetic correlation between epigenetic smoking and previously studied self-reported smoking behaviours (ranging from −0.62 to 0.72)²² can be found in Supplementary Fig. 6.

This multi-tissue, multi-cohort analysis of the relationship between smoking and DNAm (assessed via arrays and sequencing) has improved both our understanding of the biological consequences of smoking and our ability to measure it objectively. The array-based study, which identified two loci not listed in the EWAS catalog as being associated with smoking, represents the largest single cohort EWAS of smoking and the largest EPIC array EWAS of smoking, to date. The updated epigenetic biomarker of tobacco-use, mCigarette, reliably predicted smoking status and was strongly correlated with self-reported pack years of tobacco use. The analysis of sites differentially methylated in the brains of smokers and non-smokers revealed evidence of tissue-specific signals. There was a partial overlap between the results of the GrimAge DNAm pack years GWAS conducted in GS (n = 17,105) and the most extensive GWAS of self-reported smoking to date (n = 131,892).

Among the loci not present in the EWAS catalog but identified in the EPIC array EWAS of smoking at PIP > 80%, five array-based CpGs are annotated to FGF20, SCAMP5, GRIK5, SKI, and HOXA4. FGF20 plays a key role in the survival and function of dopamine-producing neurons²³, which are crucial in the brain’s reward system²⁴ that nicotine stimulates. SCAMP5 is essential for dopamine release²⁵ and the pleasurable sensations associated with smoking. It reinforces smoking behaviour and potentially make individuals more susceptible to nicotine addiction. GRIK5 encodes a glutamate receptor²⁶ which modulates dopaminergic neurons within the brain’s reward pathways, influencing how strongly these pathways respond to nicotine²⁷. A proto-oncogene called SKI regulates cell growth and apoptosis²⁸. It could potentially modify neural circuitry related to addiction, which affect the risk of developing nicotine dependence or influence the severity of addiction. Lastly, HOXA4 is a homeobox ___domain gene that is normally involved in embryonic development²⁹. Homeobox genes are abnormally expressed in cancer cells and changes in the expression of HOXA4 has been specifically associated with colorectal, ovarian and lung cancer²⁹.

The findings from the next generation sequencing EWASs, which were less well-powered, also underscored the role of smoking in carcinogenesis and disrupted neurodevelopment. Significant loci not listed in the EWAS catalog (P < 1×10⁻⁵) identified in the TWIST EWAS included sites mapping to TSPAN5, which regulates tumour suppressor gene expression³⁰; USP42, involved in head and neck cancer pathogenesis³¹; and SST, encoding somatostatin, a hormone implicated in the development of pancreatic cancer³². Significant loci revealed by ONT EWAS included CpGs associated with SEPTIN9, a tumour suppressor gene³³; TERF2, a telomeric protein linked to tumour formation and progression³⁴; and CNTNAP2, a potential marker of tumour aggressiveness in oligodendrogliomas³⁵, which is also implicated in neurodevelopmental disorders. An enrichment analysis of the TWIST results indicated that smoking influences methylation patterns across genes involved in extracellular matrix interactions, chronic inflammation, platelet activation, and ATP regulation, among other pathways. Again, the genes present in enriched pathways typically included COL4A3 and COL4A4 which play crucial roles in cancer invasion and metastasis³⁶; and inflammation and remodelling of the lung extracellular matrix³⁷.

When compared to previously published epigenetic biomarkers of tobacco use, mCigarette was a better predictor of smoking pack years. It also showed an excellent performance in discriminating smoking status (current, former, never) in two external cohorts. The robust performance of mCigarette among pregnant women in ALSPAC may reflect the unique biological context of pregnancy, where hormonal changes and epigenetic plasticity make smoking-related epigenetic changes more pronounced. Loci which are most responsive to smoking during pregnancy could provide insights into changes transmitted to the foetus, and affecting the child’s health later in life. In the future, mCigarette could be used to monitor smoking cessation efforts during pregnancy, ensuring compliance with cessation programs and potentially improving their effectiveness. The AHRR single-site biomarker provided a highly practical option for smoking classification, particularly in settings where limited data points or resources are available. However, multi-site models offered greater precision, accounting for a broader smoking-related methylation signature.

While individual CpG sites offered excellent discrimination of cases and controls, these loci varied by tissue. This is consistent with the low correlations in methylation patterns between blood and brain tissue reported previously³⁸. Future work should explore if tissue-specific signals identify pathways and mechanisms by which smoking influences brain health.

In GS, the GrimAge DNAm pack years GWAS results did not align with self-reported smoking GWAS findings but did show partial overlap of lead loci with the most extensive GWAS of self-reported smoking to date. This may suggest an increased power to detect significant loci when the epigenetic score is analysed as a phenotype. However, the genetic correlation between GrimAge DNAm pack years and the meta-analysis smoking pack years was lower than the latter and self-reported pack years in GS. The moderate correlation could reflect differences in the biological pathways that phenotypic and epigenetic measures of smoking are capturing. The GrimAge-based DNAm estimator is designed to capture the cumulative biological impact of smoking, which might include broader aging-related processes beyond direct tobacco exposure. In contrast, the smoking pack years represent a more straightforward measure of cumulative smoking exposure. The shared loci between the GrimAge DNAm pack years GWAS and previous self-reported smoking GWAS included CHRNA3, CHRNA5, and CHRNB4. These genes encode subunits of the nicotinic acetylcholine receptor, responsible for neurotransmission and binding of nicotine in the brain. Variations in these genes can affect nicotine dependence and may be associated with neurological conditions as well as lung cancer^39,40.

The key strengths of this study include the high resolution of epigenetic data used in the main EWAS analysis (both sample size and the number of measured CpG sites), diverse DNAm profiling techniques and the availability of multi-tissue DNAm. Using targeted sequencing allowed for identifying associations between smoking and loci not included on the EPIC chip. Given the mean age of GS volunteers, mCigarette represents many years of exposure to cigarette smoke, and is likely to capture long term DNAm changes associated with smoking.

Limitations of this study include its potential lack of generalizability to non-Europeans and the small number of brain- and whole methylome samples. Nevertheless, given that smoking is associated with major epigenetic alternations, the effects of tobacco use are detectable, even in small datasets. The absence of serum cotinine concentrations prevented us from comparing mCigarette against a clinically established smoking biomarker. We were also unable to verify self-reported smoking status. To mitigate potential bias, we cross-referenced the reported smoking status with other variables (such as self-reported smoking initiation and cessation) and eliminated records with conflicting responses. Additionally, we acknowledge a limitation in our study regarding second-hand smoke exposure. In family-based cohorts such as GS, where related participants may live in the same household, passive smoking is an important factor that could confound the associations between direct smoking and health outcomes. Unfortunately, data on second-hand smoke exposure were not collected in this cohort.

In conclusion, this study explored methylation patterns associated with smoking in blood and brain. The blood-based analyses, using both sequencing- and array-based approaches, identified additional loci associated with tobacco use and led to the development of a highly accurate blood-based DNAm biomarker for smoking. Furthermore, the study provided insights into the differential effects of smoking across tissues. Together, these enhance our understanding of the epigenetic architecture of smoking and shed light on the molecular mechanisms by which tobacco use influences health.

Generation Scotland

Volunteer recruitment to GS has been detailed in a previous publication⁴¹. Between 2006 and 2011, patients of collaborating general medical practices in Scotland between 35 and 65 years of age were selected at random and invited to take part in the study. These individuals were then encouraged to recruit family members to volunteer to join the cohort. A total of 24,088 individuals between 18-99 years of age completed a health questionnaire. All individuals were asked to report their smoking status (classified as non-smoker, ex-smoker who stopped more than 12 months ago, ex-smoker who stopped within the past 12 months, and current smoker). In addition, current and former smokers provided information about the age they started smoking, the age they quit smoking (former smokers only), and the number of cigarettes they smoked per day. Pack years were computed by multiplying years of smoking by the number of cigarettes smoked per day divided by 20 (number of cigarettes in a pack), and assigning a value of zero to those who never smoked. Further information about the distribution of phenotypic data is available in Supplementary Data 1. DNA extracted from whole blood collected at the baseline visit was genotyped using the Illumina HumanOmniExpressExome array (8v1-2 and 8v1) for 19,992 individuals. Following quality control (QC), imputation to the Haplotype Reference Consortium (HRC) panel⁴² and post imputation QC, 7,626,922 SNPs remained for downstream analyses. Methylation levels were measured with Illumina EPIC850k array (18,413 individuals, 752,722 sites after QC). GWAS and EWAS quality control steps have been described before¹⁰ and are also documented in Supplementary Data 15 and 16. Phenotype pre-processing steps are detailed in Supplementary Fig. 7.

Lothian Birth Cohort 1936

The Lothian Birth Cohort of 1936 (n = 1091) comprises community-dwelling older adults in Scotland, most of whom completed an intelligence test aged around 11 years in 1947. Later in life, those living in the Edinburgh and Lothians region were recruited to the cohort at a mean age of ~70 years and then followed up at 3-yearly intervals. Data collected at each wave comprised cognitive test scores as well as biological measures obtained from blood samples. During a baseline interview at age 70, the participants’ self-reported smoking status (never smoker, past smoker, current smoker) and smoking behaviour (age at starting, age at stopping, average number of cigarettes smoked per day) were determined. Pack years were calculated as in GS. A brain tissue bank was established at wave 3 (from age ~76 years). Detailed information about the cohort, brain imaging and post-mortem brain samples can be found in a cohort update and brain protocol papers^43,44. DNAm from whole blood has been measured using Illumina Infinium HumanMethylation450 BeadChip array, while DNAm in five post-mortem brain tissues was profiled using the Illumina EPIC850k array⁴⁵. Quality control and processing details are provided in Supplementary Data 15. Phenotype pre-processing steps are detailed in Supplementary Fig. 8.

ALSPAC

ALSPAC is a cohort study conducted among pregnant women residing in Avon, UK, with expected delivery dates falling between 1st April 1991 and 31st December 1992^46,47. Out of the 20,248 eligible pregnancies, 14,541 were enrolled to the study, resulting in 14,062 live births, of which 13,988 children survived to age one. During pregnancy, mothers invited fathers to take part in the study. A total of 12,113 fathers completed questionnaires, with 3807 currently formally enrolled. For a subset of ALSPAC participants (mothers, fathers and children) DNAm was assayed as part of the Accessible Resource for Integrated Epigenomic Studies (ARIES) initiative^48,49. DNA was extracted from blood samples collected at various time intervals between birth and death, and methylation levels were measured using the Illumina Infinium HumanMethylation450 or MethylationEPIC BeadChip arrays. 450,838 CpG sites passed quality control and were common to these methylation arrays. This study used four collections of DNAm data. Antenatal collection includes data from the ALSPAC mothers only, the Focus on Mothers (FOM)/ Focus on Fathers (FOF) collection corresponds to the mothers/fathers at midlife ( ~ 50 years), and the F17 and F24 collections contain ALSPAC children at ages 15–17 (time-point ‘15up’) and 24 (time-point ‘F24’), respectively⁵⁰. Smoking status for F17 was based on three questionnaires administered ages 14–17. Former smokers reported having quit at the age 17 questionnaire. Current smokers reported that they had smoked weekly in one or more questionnaires but had not quit. Never smokers reported having never smoked at least one time and never reported having smoked. Smoking status for F24 was assessed using six questionnaires administered ages 14-24. Former smokers reported having smoked regularly at some point prior to age 24, but at age 24 reported not having smoked in the previous 30 days. Current smokers reported at age 24 having smoked in the last 30 days and having smoked at least 50 cigarettes in their lifetime. Never smokers reported never having smoked at age 24. Maternal antenatal smoking was assessed by questionnaires administered at 18- and 32-weeks gestation. Former smokers reported having smoked previously but have stopped smoking for the pregnancy. Current smokers reported smoking regularly in the first trimester. Never smokers reported having never smoked before or during the pregnancy. Smoking status of FOM mothers was assessed using 15 questionnaires administered at study child ages up to age 12 and a final questionnaire at age 18. Former smokers reported having smoked on at least one questionnaire but reported not smoking on the 18 y questionnaire. Current smokers reported being current regular smokers on the 18 y questionnaire. Never smokers consistently reported never having smoked on questionnaires. Smoking status of FOF fathers was assessed using 11 questionnaires administered at study child ages up to age 12 and a final questionnaire at age 20. Former smokers reported having smoked on at least one questionnaire but reported not smoking on the 20 y questionnaire. Current smokers reported being current regular smokers on the 20 y questionnaire. Never smokers consistently reported never having smoked on questionnaires.

Data for the study were gathered and administered utilizing REDCap electronic data capture tools, which are hosted at the University of Bristol. REDCap (Research Electronic Data Capture) is a secure web application specifically designed to facilitate data capture for research⁵¹. The study website provides comprehensive details on all available data, accessible through a fully searchable data dictionary located at http://www.bristol.ac.uk/alspac/researchers/our-data/.

Inclusion & ethics

This study is based on self-reported biological sex (cross-referenced with genetic data). Detailed sex distributions are provided in Supplementary Data 1. Sex/gender was not a primary focus of this study, and analyses were conducted on the full cohort to maximize statistical power. Participants in GS, LBC1936 and ALSPAC did not receive major financial compensation for participation.

All components of GS received ethical approval from the NHS Tayside Committee on Medical Research Ethics (REC Reference Number: 05/S1401/89). All participants provided broad and enduring written informed consent for biomedical research. This study was performed in accordance with the Helsinki declaration.

Ethical approval for the LBC1936 study was obtained from the Multi-Centre Research Ethics Committee for Scotland (MREC/01/0/56) and the Lothian Research Ethics committee (LREC/1998/4/183; LREC/2003/2/29). Use of human tissue for post-mortem studies has been reviewed and approved by the Edinburgh Brain Bank ethics committee and the ACCORD medical research ethics committee, AMREC (ACCORD is the Academic and Clinical Central Office for Research and Development, a joint office of the University of Edinburgh and NHS Lothian). All participants provided written informed consent. These studies were performed in accordance with the Helsinki declaration.

Ethical approval for the ALSPAC study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees. Consent for biological samples has been collected in accordance with the Human Tissue Act (2004). Informed consent for the use of data collected via questionnaires and clinics was obtained from participants following the recommendations of the ALSPAC Ethics and Law Committee at the time. The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.

Sequencing-based approach

Whole methylome sequencing data were generated for 48 unrelated smokers and non-smokers from GS as part of this study. Two approaches were used: the TWIST methylome panel ( ~ 4 million CpG sites) and ONT sequencing ( ~ 21 million CpG sites). To ensure a robust methylation signal would be present when comparing cases against controls, only heavy smokers (12 males and 12 females chosen from a pool of 40 potential cases with tobacco use ranging from 53.9 to 87.7 pack years; see Supplementary Fig. 9) were selected as cases. Additional details on the sample selection process are provided in the Supplementary Methods. Controls were matched by age and sex using the Matchit package in R⁵², with the maximum age difference of less than 12 months (0 years) within a matched pair. Cases and controls showed a clear separation in terms of their methylation at the AHRR CpG probe cg05575921 (Supplementary Fig. 10).

Sequencing using the TWIST Human Methylome Panel was performed by the Genetics Core, Edinburgh Clinical Research Facility according to TWIST Targeted Methylation Sequencing Protocol⁵³. Sequencing using the ONT kit was performed by Edinburgh Genomics (the first 24 libraries, without basecalling) and the Genetics Core, Edinburgh Clinical Research Facility on the Oxford Nanopore PromethION 24, with R10.4.1 flow cells, running for 72 hours. Further details available in Supplementary Methods. For TWIST pre-processing of raw FASTQ files, read aligning to human reference genome (GRCh38, n = 29,401,795 total reference CpGs) with bwa-meth and quality-control of the results was performed using MethylSeq bioinformatics analysis pipeline, version 2.2.0^54,55. This analysis yielded information about the methylation level and depth of coverage (DoC) at 18,248,472 covered CpG sites.

Dorado, optimized for NVIDIA GPUs, was used for high-accuracy basecalling and modified base detection in raw ONT data. Reads were aligned to GRCh38 human reference genome. Variants were called with epi2me-labs/wf-human-variation nextflow pipeline, version 23.10.1. Methylation level and depth of coverage was measured at 28,989,402 covered CpG sites.

The bedGraph (TWIST) and bedMethyl (ONT) files were subsequently processed in R version 4.3.1⁵⁶ using Methrix package⁵⁷. As part of post-processing, loci a) of extremely low (minimal DoC = 2) and high coverage (beyond 0.99 quantile), and b) overlapping with known cytosine to thymine polymorphisms were removed from the methylation dataset. Finally, a coverage filter was applied, retaining only the loci that were covered in at least 40 samples by: a) 10 or more reads (TWIST), b) 5 or more reads (ONT). This left an analysis sample of 3,391,718 and 21,167,712 CpGs for TWIST and ONT data, respectively. CpG sites were annotated with Annotatr R package⁵⁸.

Blood-based DNAm EWAS

A blood-based EWAS of smoking was carried out in 17,865 GS individuals using BayesR⁵⁹. Before running the EWAS, each CpG was corrected for the effects of age, sex, and batch using linear regression (saving the residuals from each model as the new variable). As BayesR implicitly corrects for confounding effects without requiring a full characterization of all covariates, our models were only adjusted for measured variables i.e., estimated white blood cell proportions were not included. However, we conducted a sensitivity analysis by running an EWAS that included adjustments for estimated cell proportions. The results of this sensitivity analysis showed a strong overlap with the primary findings and are presented in Supplementary Data 17. The smoking phenotype (measured in pack years, with a pack defined as 20 cigarettes) was natural log+1 transformed and adjusted for age and sex using linear regression (again, the residuals were saved and used for the downstream analyses). Both the smoking and the CpG variables were scaled to have mean of 0 and variance of 1. The data served as inputs of a Bayesian penalised linear regression model. Four Gaussian priors were specified to model CpGs with varying effect sizes (mixture variances of 0.1%, 1%, 10% and 100%) along with a discrete spike at the origin to model CpGs with no effect. A Gibbs sampler was used to sample over the posterior distribution, conditioning on the input data. A burn-in of 5000 samples was used, after which every fifth sample was retained across 10,000 iterations. A CpG with a posterior inclusion probability of greater than 0.95 was considered as epigenome-wide significant. Previously unpublished associations were identified by searching the literature and the EWAS catalog¹².

Comparison between TWIST, ONT and EPIC850k

DNAm of 24 pairs of smokers and non-smokers from GS was profiled using the EPIC array, ONT kit, and TWIST platform (see Methods - Sequencing-based approach). During quality control, one pair was identified as mismatched due an age gap exceeding the pre-defined threshold of 12 months. To preserve the integrity of the study design and prevent potential biases in the analysis, this pair was excluded from subsequent analyses. This left an analysis sample of 46 individuals. For each DNAm profiling method, an EWAS of smoking (pack years) was conducted. The association between DNAm level at each CpG site (outcome) and binary smoking status, age and sex was modelled using linear regression. The results were displayed on a Manhattan plot generated with the CMplot R package⁶⁰. Gene names associated with CpG sites reaching a suggestive significance threshold (P < 1x10⁻⁵) were extracted and subjected gene set enrichment analysis in Functional Mapping and Annotation (FUMA) GENE2FUNC tool⁶¹, which implements a hypergeometric test. An FDR-adjusted p-value threshold of 0.05 was applied, and a minimum of 2 overlapping genes within each gene set was required.

Biomarkers of cumulative smoking

An elastic net biomarker of pack years (mCigarette) was trained in 17,865 GS individuals using the glmnet library in R⁶². As part of data pre-processing, CpG sites were filtered to 18,760 loci associated with smoking at False Discovery Rate (FDR) < 0.05 in a previous meta-analysis EWAS (n = 18,760) of tobacco use⁶. Alpha was fixed at 0.5 and the lambda value that minimised the mean prediction error was selected via 10-fold cross validation. The selected model assigned non-zero coefficients to 1255 CpGs. A single-site biomarker for smoking, based on methylation at AHRR (cg05575921), was also trained in the same subset of GS individuals (n = 17,865) using linear regression. Both mCigarette and the single-site biomarker were tested in wave 1 of LBC1936 (n = 882, mean age 70 years), while mCigarette alone was tested in ALSPAC (n = 496–1207 across four time points). These biomarkers were benchmarked against three epigenetic scores for smoking: EpiSmokEr score^7,63,64, and two scores derived from previous GS analyses on smaller subsets of the dataset - one based on BayesR weights⁵⁹, the other via lasso penalised regression by McCartney et al.¹⁰. In LBC1936, the predictive performance of mCigarette was additionally compared to that of GrimAge DNAm pack years¹⁴. Pearson’s r was calculated to estimate the degree of correlation between self-reported pack years of smoking and the studied scores. The amount of variance in pack years explained by the studied scores was assessed by comparing R² estimates of null and full models. While the null model was adjusted for age and sex, the full model also included the studied score. Incremental R² was calculated as the difference between variance explained by null and full models.

The ability of the scores to distinguish between current, former, and never smokers was assessed by AUC. Receiver operating characteristic (ROC) curves were produced using pROC R package⁶⁵. Additional prediction performance metrics such as PRAUC were obtained using MLmetrics R package⁶⁶.

Tissue specificity analyses in LBC1936

DNAm was measured in 5 brain regions (hippocampus - BA35, dorsolateral prefrontal cortex – BA46, primary visual cortex - BA17, anterior cingulate cortex - BA24, ventral/lateral inferior temporal cortex - BA20/21) from post-mortem brain samples of 14 LBC1936 individuals, with one sample missing from hippocampus. Tissue acquisition and processing details are detailed in Stevenson et al.⁴⁵. Using blood-DNAm measured at wave 1 and brain-DNAm data, exploratory EWAS analyses were performed (CpG ~ smoking category). In these analyses, smoking was treated as a continuous variable encoded as 0 = never smoker, 1 = former smoker, 2 = current smoker. Given the small sample size, nominally significant CpG-smoking associations were defined as having P < 1x10⁻⁵ and were displayed using a ggplot2 boxplot. Due to the same constraints, an enrichment analysis for significant associations was not carried out.

GWAS of smoking

Associations between genetic variants and smoking were examined using GWASs. Two phenotypes were considered: natural log(transformed pack years of smoking + 1) and an epigenetic score for smoking pack years generated by an online calculator that uses the algorithm derived for the GrimAge epigenetic clock¹⁴. After initial filtering (see Supplementary Data 3), there were 17,105 GS individuals with data available at 7,626,922 SNPs. Both GWASs were conducted using the GCTA software⁶⁷, with a Genetic Relationship Matrix fitted into a fastGWA-lmm model to account for relatedness. Each trait was adjusted for age, sex, and 20 genetic principal components. The results of these analyses, along with the findings of previous GWASs of tobacco use, were visualised using CMplot⁶⁰. Lead and methylation quantitative loci among the results of GrimAge DNAm pack years GWAS were identified using the default settings in FUMA⁶¹ and goDMC⁶⁸, respectively. GWAS catalog was accessed via FUMA website. Genetic correlations were calculated with LDSC⁶⁹.