Chronic obstructive pulmonary disease (COPD) is one of the major causes of death worldwide. It is characterised by various degrees of airway remodelling and emphysematous lung tissue destruction. Previously, we and others have shown large transcriptional differences between COPD and non-COPD individuals in different compartments of the lungs [1]. Interestingly, even within COPD, there is large heterogeneity within these transcriptional profiles, which predict the severity of the disease [2] and the response to current therapies [3]. It is still unclear how these expression profiles are regulated, however epigenetic regulation likely plays an important role. Although many studies have found differences in DNA methylation in blood of COPD patients [4], few studies investigated the direct association between DNA methylation in the lung itself and presence or severity of COPD.

In the current study, we aimed to identify differences in airway wall DNA methylation between current smokers with and without COPD, and to assess how differences in DNA methylation are associated with the severity of airflow obstruction and hyperinflation. We also investigated the impact of differential DNA methylation on gene expression in the airway wall, to elucidate its role in COPD pathology.

We used bronchial biopsy samples from the previously described Stop-Smoking [5] (nCOPD=14, nnon-COPD=21), and Groningen and Leiden Universities Corticosteroids in Obstructive Lung Disease (GLUCOLD [6], nCOPD=76, ClinicalTrails.gov: NCT00158847) studies, where patients underwent bronchoscopy voluntarily. Patients with Alpha-1 Antitrypsin Deficiency were excluded from both studies. The local medical ethics committees approved all studies, and all subjects gave their written informed consent. DNA methylation profiling was performed using the HumanMethylation850K BeadChip (Illumina Inc., San Diego, CA, USA) and gene expression profiling using RNA-Seq [7]. The methods for DNA and RNA extraction, 850K array processing and RNA-Seq were described previously [7]. Raw DNA methylation intensity values were processed using the Minfi package [8] in R. Samples and probes failing QC were removed, and raw beta values normalized using the Dasen method as implemented in the WateRmelon package [9]. Raw RNA-Seq reads were aligned using STAR-Aligner and normalized to counts per million. Technical variation between the RNA-Seq datasets from both studies was eliminated using CombatSeq [10].

We performed differential methylation analysis comparing COPD patients and non-COPD controls using robust linear modelling, adjusting for age and sex as potential confounders. As it is well established that active smoking has significant effects on the epigenome of the airways [11] and all non-COPD controls were active smokers, we decided to only include current smoking COPD patients in our analysis, resulting in nCOPD=61 and nnon-COPD=21. We used a robust linear model comparing severity of airflow obstruction and hyperinflation as reflected by the FEV1% predicted (forced expiratory volume in 1s) and RV/TLC % predicted (residual volume/total lung capacity) in COPD patients only, adjusting for age and sex. eQTM analysis (<250,000bp) was conducted in samples with both methylation and RNA-Seq data, separately for COPD-associated CpGs and CpGs associated with disease severity. The five principal components explaining 95% of the variation in the control probes of the 850K array were used as covariates to correct for technical variation in all analyses. A Benjamini-Hochberg False Discovery Rate (FDR) <0.05 was maintained to define statistical significance. Cell proportions in bronchial biopsies were estimated using CIBERSORT as described previously [7]. M-values of DNA-methylation sites were Spearman correlated to cell proportions and a p-value <0.05 was considered significant.

The clinical characteristics of the included COPD patients (n=61) and non-COPD controls (n=21) are displayed in Table 1. We identified six differentially expressed methylation sites between COPD patients and non-COPD controls (Fig. 1A), of which three were hypermethylated (cg23800199, cg05748572, cg16560488) and three hypomethylated (cg10730398, cg01124573, cg12081267) (Fig. 1B) in COPD. Five of the six CpG sites were found in the intronic regions of protein-coding genes (DIRC3, GNG7, JMJD1C, ADGRL3 and TMEM131 respectively), which could contribute to gene splicing [12]. Of interest, cg01124573 is located within an intronic region of ADGRL3, a gene previously found to be highly expressed in asthma and associated with contraction of airway smooth muscle when activated [13]. Further, the CpG site cg16560488 was hypermethylated in COPD and located in JMJD1C, a histone H3K9 demethylase gene [14]. This finding is of interest as a previous GWAS study showed the SNP rs7899503, for JMJD1C, to be associated with the level of FEV1 in a multi-ethnic meta-analysis [15]. These findings suggest a role for histone modification, often a result of cigarette smoke exposure and high oxidant levels [16], in COPD pathogenesis. eQTM analysis showed no correlation between the degree of methylation of these 6 CpG sites and the level of gene expression of nearby genes. Even though cg01124573 was not found to be an eQTM for ADGRL3, its hypomethylation could still be important in the context of COPD as it may result in a different splicing variant of this gene.

Fig. 1.

(A) Volcano plot showing differential methylation in COPD. CpG sites hypermethylated in COPD are red; hypomethylated CpG sites are blue. (B) Boxplots showing methylation in beta values of the six CpG sites differentially methylated in COPD, stratified over disease status. Nominal p-values from differential methylation analysis results are shown. (C and D) Volcano plots of correlation of methylation with FEV1% predicted and RV/TLC % predicted, respectively. (E and G) Scatterplots showing the correlation of methylation with FEV1% predicted and RV/TLC % predicted, respectively, of the top eQTMs. logFC, beta values and nominal p-values from differential methylation results are shown. (F and H) Scatterplots showing the top eQTM for the previous CpG site. Beta values and nominal p-values from eQTM results are shown. (I) Cellular proportions calculated using CIBERSORT of the cell types representing each >2% of the cell population in the bronchial biopsies, stratified over disease status. Purple dots indicate individuals with COPD and green dots non-COPD controls. p-Values from a Wilcoxon rank-sum test are shown.

In patients with COPD, we found 25 CpG sites associated with the severity of airflow obstruction as reflected by the level of FEV1% predicted (Fig. 1C) and 3954 CpG sites associated with level of hyperinflation as reflected by RV/TLC % predicted (FDR <0.05) (Fig. 1D). eQTM analyses identified 2 eQTMs for CpGs associated with FEV1% predicted, including cg06445873 which was hypermethylated with a lower FEV1% predicted and associated with lower expression of HSD11B2 and TPPP3 (Fig. 1E and F). 817 eQTMs were identified for CpGs associated with RV/TLC % predicted. Two of the top CpGs, cg11734183 and cg13421038, were hypermethylated in association with higher RV/TLC % predicted and lower expression of the genes RGMB and SFSWAP, respectively (Fig. 1G and H).

To confirm that these findings are COPD specific and not related to smoking status, we additionally conducted a differential methylation analysis with ex-smokers included, correcting for age, sex and smoking status (nCOPD=90, non-COPD=21). All CpG sites that we identified as differentially methylated in current smokers in association with COPD, FEV1 and RV/TLC, were methylated in the same direction when ex-smokers were included, with 95% of sites remaining significant.

To investigate if the eQTMs were associated with altered cell proportions, we performed cellular deconvolution on bronchial biopsies. Goblet cells, cycling basal cells, ciliated cells, and fibroblasts had the highest predicted cell proportions in the airway wall biopsies, each having more than 2% of the total predicted cell population and therefore were included in subsequent analyses (Fig. 1I). Goblet cells, cycling basal cells and ciliated cells showed a significant difference in cell proportion between COPD patients and non-COPD controls. A lower FEV1% was associated with lower ciliated and goblet cell proportions, whereas a higher RV/TLC % was associated with lower goblet cell proportions. We then associated the cell types with the gene expression levels of our eQTMs. We found that 56 out of 723 eQTMs were associated with cell types, possibly indicating an association of these eQTMs with cell population shifts.

Cg06445873 hypermethylation was associated with a lower FEV1% predicted and lower HSD11B2 gene expression. HSD11B2 enzymatic activity was previously found to convert cortisol into the inactive cortisone through increasing 11-beta-dehydrogenase activity [17]. This activity is not specific to cortisol, but can have an effect on inactivation of inhaled corticosteroids used by COPD patients. Possibly, a higher efficacy of ICS due to hypermethylation of cg06445873 could contribute to less airflow obstruction in COPD patients. Additionally, SFSWAP, hypermethylated and less expressed with higher RV/TLC % predicted, regulates splicing and has been associated with regeneration and cellular differentiation [18], which is of interest given the chronic damage and repair response in COPD. The relation with splicing is also of interest as this could impact functionality of many genes and proteins involved in COPD pathology.

Cellular deconvolution revealed the presence and severity of COPD to be associated with lower predicted goblet cell proportions, which is unexpected as goblet cells have been shown to be increased in COPD. Possibly, the shift in goblet cell proportion we observed does not reflect lower absolute numbers in COPD but rather may be related to a relative increase in basal cycling cells in response to the repair of damaged epithelium in COPD.

There are some limitations in the current study. Bronchial biopsies for methylation and gene expression, though collected on the same day, were serially collected from adjacent locations within the bronchial tree. This could potentially affect the cellular composition and subsequently the concordance between methylation and gene expression, limiting our ability to identify eQTMs. Additionally, an objective test to confirm smoking status was not routinely performed. While we identified many COPD-associated epigenetic changes, it should be acknowledged this was a post-hoc analysis and the sample size was relatively small.

In conclusion, we conducted a genome-wide screening to investigate how altered DNA methylation may contribute to presence and severity of COPD. We identified CpG sites that regulate the expression of genes involved in airway smooth muscle contraction, corticosteroid metabolism and regeneration and cellular differentiation. These findings highlight that airway changes in COPD patients are associated with differential DNA methylation, specifically in the airway wall. Expanding upon prior investigations in other tissue types, our study provides direct, lung-specific insight into potential mechanisms underlying COPD pathogenesis.