Tobacco smoking induces changes in true DNA methylation, hydroxymethylation and gene expression in bronchoalveolar lavage cells

Abstract

Background: While smoking is known to associate with development of multiple diseases, the underlying mechanisms are still poorly understood. Tobacco smoking can modify the chemical integrity of DNA leading to changes in transcriptional activity, partly through an altered epigenetic state. We aimed to investigate the impact of smoking on lung cells collected from bronchoalveolar lavage (BAL).

Methods: We profiled changes in DNA methylation (5mC) and its oxidised form hydroxymethylation (5hmC) using conventional bisulphite (BS) treatment and oxidative bisulphite treatment with Illumina Infinium MethylationEPIC BeadChip, and examined gene expression by RNA-seq in healthy smokers.

Findings: We identified 1667 total 5mC + 5hmC, 1756 5mC and 67 5hmC differentially methylated positions (DMPs) between smokers and non-smokers (FDR-adjusted P <.05, absolute Δβ >0.15). Both 5mC DMPs and to a lesser extent 5mC + 5hmC were predominantly hypomethylated. In contrast, almost all 5hmC DMPs were hypermethylated, supporting the hypothesis that smoking-associated oxidative stress can lead to DNA demethylation, via the established sequential oxidation of which 5hmC is the first step. While we confirmed differential methylation of previously reported smoking-associated 5mC + 5hmC CpGs using former generations of BeadChips in alveolar macrophages, the large majority of identified DMPs, 5mC + 5hmC (1639/1667), 5mC (1738/1756), and 5hmC (67/67), have not been previously reported. Most of these novel smoking-associating sites are specific to the EPIC BeadChip and, interestingly, many of them are associated to FANTOM5 enhancers. Transcriptional changes affecting 633 transcripts were consistent with DNA methylation profiles and converged to alteration of genes involved in migration, signalling and inflammatory response of immune cells.

Interpretation: Collectively, these findings suggest that tobacco smoke exposure epigenetically modifies BAL cells, possibly involving a continuous active demethylation and subsequent increased activity of inflammatory processes in the lungs. FUND: The study was supported by the Swedish Research Council, the Swedish Heart-Lung Foundation, the Stockholm County Council (ALF), the King Gustav's and Queen Victoria's Freemasons' Foundation, Knut and Alice Wallenberg Foundation, Neuro Sweden, and the Swedish MS foundation.

Keywords: Alveolar macrophages; DNA hydroxymethylation; DNA methylation; EPIC; Enhancers; Epigenetics; Oxidative stress; Smoking.

Fig. 1

Feature-specific differences in β-value distribution between total 5mC + 5hmC and true 5mC. Violin plots representing the distribution of β-values from 735,794 probes after SQN normalisation, plotted as densities in relation to CpG islands (a) and genomic features (b). Yellow plots represent standard bisulphite-treated samples with a combination of both 5mC and 5hmC (total BS methyl). Green plots represent oxidative bisulphite-treated samples with 5mC (true 5mC, oxBS). Plots representing smoker densities are shown in upper rows (a-b) and non-smoker plots in lower rows (a-b). Differences in distributions between 5mC and 5mC + 5hmC was tested using Wilcoxon signed-rank test was considered significant with a P value <.05. * P < .05; ** P < .01; *** P < .001; n.s = not significant. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Smoking-associated DMPs are predominantly hypomethylated in 5mC and total 5mC + 5hmC with enrichment in gene bodies, non-CGI context and enhancer sites. Horizontal bar plots (a-d) illustrating relative frequencies of DMPs associated with smoking. (a) Relative frequencies of hypermethylated and hypomethylated 5mC + 5hmC, 5mC, and 5hmC DMPs. Percentage of hypermethylated and hypomethylated 5mC + 5hmC, 5mC and 5hmC across (b) CGI-related features (CpG islands, shores, shelves, open sea), (c) gene features (TSS1500, TSS200, 1stExon, 5’UTR, Body, 3’UTR, ExonBnd, IGR), and (d) enhancers. The distribution of all EPIC array probes included in our analysis are shown as EPIC bkg (background) for comparison (b-d). Enrichment/depletion analysis was performed using Chi-square test on frequencies, adjusting P values for multiple testing (Bonferroni). * P < .05 compared to non-smokers, ** P < .01 compared to non-smokers, *** P < .001 compared to non-smokers.

Fig. 3

Hypomethylation and hyperhydroxymethylation in smokers. Overlapping total 5mC + 5hmC (yellow), 5mC (green) and 5hmC (purple) with adjusted P value <.05 (absolute Δβ threshold: 5mC + 5hmC and 5mC >0.15; 5hmC >0.05). (a) Venn diagram illustrating number of DMPs and overlaps between 5mC + 5hmC, 5mC, and 5hmC. (b) Boxplot showing 5mC + 5hmC, 5mC and 5hmC Δβ for the 10 overlapping DMPs between 5mC and 5hmC. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

Correlation of promoter, gene body, and enhancer methylation with gene expression. Plots showing genes with differences in both DNA methylation and expression in promoters (a-b), gene body (c-d), and enhancers (e-f) of smokers-associated DMPs and genes. Scatterplots of mean gene expression values (log2 fold change) and total methyl (5mC + 5hmC) mean Δβ in promoter (a), gene body (c), and enhancer (e). Orange dots represent genes with significant smoking-associated changes in both total 5mC + 5hmC (>0.15 absolute Δβ) and gene expression (log2 fold change >1), green dots represent significant expression but not methylation, and purple dots represent changes in 5mC + 5hmC but not gene expression (b,d,f). Correlation plots of selected genes showing both methylation (5mC + 5hmC, β-values) and gene expression (normalised gene count), with Pearson correlation coefficient and P value (b,d,f). Smokers (n = 7) are represented by yellow dots, and (n = 12) non-smokers by blue dots. Promoters are represented by TSS200 and TSS1500 CpG sites. Enhancer sites are present at various genomic features and can overlap with both promoter and gene body sites. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

Functional annotations of differentially expressed genes. All (a) and top (b) biological processes and diseases associated with differentially expressed genes obtained with Ingenuity Pathway Analysis (IPA). Top canonical pathways (c) associated with differentially expressed genes obtained with IPA. (d) Schematic representation of the top gene interaction network obtained with IPA, with downregulated and upregulated genes depicted in blue and red, respectively. (a-c) Significance is represented as –log₁₀P value and colours indicate predicted activation z-score, with decreased, no effect, and increased activation in blue, white and red colours respectively. n.a. prediction not available. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

Functional annotations of differentially methylated and expressed genes. (a) Multidimensional scaling of GO Biological processes terms associated with differentially methylated (BS-DMP: total 5mC + 5hmC) (green), differentially expressed (DE) (orange) genes according to semantic similarities. Findings from annotated DE genes that harbour BS-DMPs are depicted in purple. GO terms were obtained using over-representation analysis and visualization was generated by REVIGO. (b) Top common biological processes and diseases associated with differentially methylated (BS-DMP) (green) and DE (orange) genes obtained with Ingenuity Pathway Analysis (IPA), with IPA findings from genes displaying both BS-DMP (FDR < 0.05, Δβ > 0.15) and transcriptional change in purple. (c) Top canonical pathways associated to both DE- and DMP-genes. (d) Representation of the genes network obtained with IPA on both DE- and DMP-genes and visualised using STRING. Grey line gradient and thickness indicates the strength of data support (darker thick grey representing stronger confidence) and colours represent different cluster (Markov Clustering set at 6). (e) Top upstream regulators for both differentially methylated and expressed genes, with colours depicting different classes of regulators. Significance is represented as –log₁₀P value. DE, differentially expressed, BS-DMPs, bisulphite-generated differentially methylated positions (5mC + 5hmC). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)