Chronic Obstructive Pulmonary Disease Is Associated with Epigenome-Wide Differential Methylation in BAL Lung Cells.

DNA methylation patterns in chronic pulmonary obstructive disease (COPD) might offer new insights into disease pathogenesis. To assess methylation profiles in the main COPD target organ, we performed an epigenome-wide association study on BAL cells. Bronchoscopies were performed in 18 subjects with COPD and 15 control subjects (ex- and current smokers). DNA methylation was measured using the Illumina MethylationEPIC BeadChip Kit, covering more than 850,000 CpGs. Differentially methylated positions (DMPs) were examined for 1) enrichment in pathways and functional gene relationships using the Kyoto Encyclopedia of Genes and Genomes and Gene Ontology, 2) accelerated aging using Horvath's epigenetic clock, 3) correlation with gene expression, and 4) colocalization with genetic variation. We found 1,155 Bonferroni-significant (P < 6.74 × 10-8) DMPs associated with COPD, many with large effect sizes. Functional analysis identified biologically plausible pathways and gene relationships, including enrichment for transcription factor activity. Strong correlation was found between DNA methylation and chronological age but not between COPD and accelerated aging. For 79 unique DMPs, DNA methylation correlated significantly with gene expression in BAL cells. Thirty-nine percent of DMPs were colocalized with COPD-associated SNPs. To the best of our knowledge, this is the first epigenome-wide association study of COPD on BAL cells, and our analyses revealed many differential methylation sites. Integration with mRNA data showed a strong functional readout for relevant genes, identifying sites where DNA methylation might directly affect expression. Almost half of DMPs were colocated with SNPs identified in previous genome-wide association studies of COPD, suggesting joint genetic and epigenetic pathways related to disease.

Clinical Relevance

This study integrates novel findings of widespread differential methylation in BAL cells from subjects with chronic obstructive pulmonary disease (COPD) compared with control subjects with data from previous transcriptomic and genome-wide association studies. Taken together, results suggest that at some of the COPD-associated sites, DNA methylation might serve as a mediator between genetics and disease and at others be acquired alterations with functional relevance. As BAL collects cells from the parts of the lung most affected by the disease, this study might contribute to insights into COPD pathogenesis.

Chronic obstructive pulmonary disease (COPD) is a heterogeneous chronic inflammatory airway disease, characterized clinically by airflow limitation, chronic respiratory symptoms, and many systemic comorbidities. It has long been known that the key environmental risk factor for COPD is tobacco smoking, with smoke from biomass fuels and other pollutants as other major contributors. More recently, it has become clear that a failure to reach the predicted degree of peak lung function in early adulthood might also lead to chronic airflow limitation, even in the absence of environmental exposure (1, 2). According to the World Health Organization, COPD is now the third leading cause of death in the world (3).

Why some, but not all, long-time smokers develop COPD is a major yet largely unresolved question in COPD research. Large-scale genome-wide association studies (GWASs) have revealed several genetic variants related to COPD (4), and according to twin studies, ∼40–60% of COPD susceptibility is explained by genetic factors (5, 6). What additional role epigenetics in general, and DNA methylation in particular, might play in COPD pathogenesis has not been thoroughly investigated. Previous epigenome-wide association studies (EWAS) are scarce, especially those based on cells from the main target organ of the disease. In a systematic review from 2017, Machin and colleagues found no consistent differences in DNA methylation in peripheral blood (7), whereas the few existing EWASs of lung tissue and bronchial brushings have more consistently reported aberrant methylation in COPD (8–10). Different methylation patterns between target organ and surrogate tissue would be in line with what has previously been shown in transcriptomic studies of COPD (11).

The distal airways represent the main site of airflow limitation in COPD and are the parts of the lung most affected by the disease (12). To investigate DNA methylation patterns in cells from distal airways, we performed an EWAS on BAL cells from subjects with COPD and control subjects. We subsequently investigated potential functional effects of differential methylation through integration with GWAS and gene expression data and also tested the hypothesis of accelerated aging in COPD, using Horvath’s epigenetic clock (Figure 1).

Figure 1.

Study overview. GO = Gene ontology; GREAT = Genomic Regions Enrichment of Annotations Tool; GWAS = genome-wide association study; KEGG = Kyoto Encyclopedia of Genes and Genomes; UCSC = The University of California, Santa Cruz, Genome Browser Database.

Methods

Study Subjects

Twenty-two subjects with COPD and 15 ever-smokers with normal lung function participated in this cross-sectional study. All subjects were recruited from the longitudinal OLIN (Obstructive Lung Disease in Northern Sweden) COPD study using predetermined criteria (13); details are provided in the data supplement.

Bronchoscopy

All subjects underwent bronchoscopy with BAL, as described in the supplement. BAL could not be performed in three subjects with COPD, because of problems tolerating bronchoscopy. In one subject with COPD, BAL recovery was too low to perform analysis.

Statistical Analyses

We performed EWAS analysis to estimate association between DNA methylation β values and COPD using robust linear regression (rlm in the MASS R package [14]). The analyses were adjusted for sex, age (years), inhaled corticosteroid use, pack-years, smoking intensity, and measured cell type count (macrophages, lymphocytes, neutrophils, and eosinophils). Covariate selection was aided by a principal-component regression analysis (see Figure E1 in the data supplement). We applied Bonferroni correction (P < 6.74 × 10−8) to account for multiple testing in EWAS results. Differentially methylated regions (DMRs) were identified using DMRcate (15). Input parameters for the DMR calling algorithms are provided in the data supplement.

Functional Gene Relationship and Pathway Enrichment Analysis

Functional enrichment analysis in gene ontology and pathway enrichment in Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using the ConsensusPathDB (16, 17) overrepresentation analysis tool (http://consensuspathdb.org). P values for enrichment were adjusted for multiple testing using the false discovery rate method.

Epigenetic Age

DNA methylation age was calculated using the new DNA Methylation Age Calculator by Horvath (18). Details on preprocessing and calculator options can be found in the data supplement.

Integration with GWAS

COPD-associated differentially methylated positions (DMPs) were compared with the 82 SNPs from a recent large-scale GWAS meta-analysis (4). We investigated if DMPs in the present study were located within a window 500 kb upstream and 500 kb downstream of these SNPs. Moreover, we checked whether SNPs and DMPs were located in the same gene.

DNA Methylation in Relation to Gene Expression

Correlations between DNA methylation and amounts of gene expression were tested using a previously published paired DNA methylation and gene expression BAL cell data set (19). We tested transcript concentrations of genes within a 500-kb region of the 1,155 DMPs (250 kb upstream and 250 kb downstream). The mRNA gene expression (RNA sequencing [RNA-seq]) and methylation (Illumina Infinium MethylationEPIC BeadChip Kit) were measured in BAL samples from 7 smokers and 12 nonsmokers (19). First, we created residuals for mRNA expression and residuals for DNA methylation and used Pearson correlation to evaluate correlations between expression residuals and DNA methylation residuals. These residual models were adjusted for covariates (age, sex, and smoking status). We corrected these analyses for multiple testing using Bonferroni correction. In addition, we explored the Immunological Genome Project (ImmGen) RNA-seq data set to evaluate if differentially methylated genes were likely to originate from macrophages or other cell types (20).

Results

DNA methylation data were available for 18 subjects with COPD and 15 control subjects with normal lung function (Table 1). There were no significant differences between cases and control subjects by age or body mass index. As expected, inhaled corticosteroids were used only by subjects with COPD. Pack-years were higher in the COPD group, and consistent with the results of previous studies (21), BAL recovery was lower. Both groups consisted of current and ex-smokers. Among current smokers, smoking intensity was higher in subjects with COPD. Differential cell counts revealed no significant differences between groups, and as expected, the predominant cell type in BAL was macrophages (Table 1).

Table 1.

Study Subjects

	COPD (n = 18)	Smokers with Normal LF (n = 15)	P Value
Female:male	4:14	8:7	NS
Age, y	63 (62–72)	65 (63–71.5)	NS
BMI, kg/m2	26.2 (23.7–29.0)	26.0 (24.2–27.8)	NS
Current smokers:ex-smokers	10:8	3:12	NS
Smoking intensity (among current smokers), cigarettes/d	14 (10–20)	8 (7–10)	0.049
Pack-years	36 (30–40)	15 (13.5–21.5)	<0.0001
FEV1% predicted	68.3 (44.1–72.3)	110.2 (94.4–117.9)	<0.0001
FEV1/VC	0.53 (0.43–0.62)	0.74 (0.70–0.77)	<0.0001
Use of inhaled corticosteroids, yes:no	6:12	0:15	0.02
BAL recovery, %	43.6 (32.2–55.6)	61.1 (56.1–66.7)	0.001
BAL macrophages, %	89.5 (86.4–92.8)	86.4 (81.3–92.6)	NS
BAL lymphocytes, %	8.3 (5.2–13.2)	11.8 (6.1–18.2)	NS
BAL neutrophils, %	1.0 (0.4–1.8)	0.6 (0.3–1.2)	NS
BAL eosinophils, %	0.3 (0.0–0.6)	0.2 (0.0–0.3)	NS
BAL mast cells, %	0.2 (0.0–0.5)	0.0 (0.0–0.1)	NS

Definition of abbreviations: BMI = body mass index; COPD = chronic obstructive pulmonary disease; FEV1 = forced expiratory volume in 1 second; LF = lung function; NS= not significant; VC = vital capacity (defined as the highest value of forced and slow vital capacity).

Values are given as median and interquartile range unless indicated differently. Statistical comparisons: for ratios, the Fisher exact test was used, and for other comparisons, the Mann-Whitney U test was used; a P value of less than 0.05 was considered to indicate significance. Pack-years were calculated as (number of cigarettes smoked per day/20) × number of years smoked.

Differential Methylation: Positions and Regions

DNA methylation data were analyzed for DMPs and DMRs.

We found 1,155 Bonferroni-significant DMPs associated with COPD, spanning all chromosomes (chrs) (see Figure 2 and Figure E3). DMPs were enriched for localization to CpG islands (23% of COPD-associated CpGs, 19% of all CpGs on the array; Penrichment = 0.0004) and relatively depleted in shores (15% of COPD-associated CpGs, 18% of all CpGs; Penrichment = 0.01; see Figure E3). Fifty-five percent of DMPs (638 of 1,155) displayed higher mean methylation in subjects with COPD relative to control subjects.

Figure 2.

Manhattan plot indicating the associations between COPD and DNA methylation. Black line represents the Bonferroni-corrected threshold (P < 0.05/[number of CpGs analyzed] ≈ 6.74 × 10−8). Chr = chromosome; COPD = chronic obstructive pulmonary disease.

The DMP with the lowest P value was annotated to CLTA (clathrin light chain A) on chr 9 (P = 6.00 × 10−39 for cg00202130; Table 2) and had a 1.1% higher mean methylation degree in subjects with COPD relative to control subjects. The largest effect size was observed for cg14640588, annotated to ZDHHC14 (zinc finger DHHC-type palmitoyltransferase 14) on chr 6, with 28.9% higher mean methylation in COPD (Table 3). Among DMPs with lower mean methylation in COPD, cg15969227 (nearest gene ZNF322 [zinc finger protein 322] on chr 6) had the largest effect size, with a mean methylation difference of 27.6% (Table 3). NUB1 (negative regulator of ubiquitin like proteins 1) on chr 7 was the gene with the largest proportion of significant DMPs annotated to it (3 differentially methylated CpGs/75 total CpGs = 4.0%), all with higher mean methylation in COPD.

Table 2.

The Top 10 Bonferroni-Significant Differentially Methylated Positions Associated with Chronic Obstructive Pulmonary Disease, Ranked by P Value

CpG ID	Chr	Genomic Coordinate	Annotated to Gene	Coefficient	P Value
cg00202130	9	36190707	CLTA1 *	−0.0113	6.00 × 10−39
cg27079832	15	64790751	ZNF609 *	0.0207	2.72 × 10−38
cg07088935	2	1426845	TPO *	0.1035	2.25 × 10−36
cg06958223	6	150499199	PPP1R14C *	0.0392	3.39 × 10−35
cg16928994	2	178977002	RBM45 *	−0.0150	1.93 × 10−33
cg11619181	9	136500618	DBH *	−0.0416	1.88 × 10−32
cg05738196	6	26577821	ABT1 †	0.0306	3.47 × 10−32
cg15969227	6	26745279	ZNF322 †	0.2759	9.23 × 10−28
cg14948709	16	25248396	ZKSCAN2 *	0.0468	6.77 × 10−27
cg15803210	5	31757766	PDZD2 †	−0.0445	8.31 × 10−26

Definition of abbreviations: ABT1 = activator of basal transcription 1; Chr = chromosome; CLTA1 = clathrin light chain A sequence 1; DBH = dopamine β-hydroxylase; ID = identifier; PDZD2 = PDZ domain containing 2; PPP1R14C = protein phosphatase 1 regulatory inhibitor subunit 14C; RBM34 = RNA binding motif protein 45; TPO = thyroid peroxidase; ZKSCAN2 = zinc finger with KRAB and SCAN domains 2; ZNF322 = zinc finger protein 322; ZNF609 = zinc finger protein 609.

Adjusted for sex, age, cell type, inhaled corticosteroid use, pack-years, and smoking intensity.

Illumina annotation based on The University of California, Santa Cruz (UCSC) database.

Nearest gene as reported by Genomic Regions Enrichment of Annotations Tool (GREAT) version 4.0.4.

Table 3.

The Top 10 Bonferroni-Significant Differentially Methylated Positions Associated with Chronic Obstructive Pulmonary Disease, Ranked by the Magnitude of Positive and Negative Effects (5 Differentially Methylated Positions Each)

CpG ID	Chr	Genomic Coordinate	Annotated to Gene	Coefficient	P Value
cg15969227	6	26745279	ZNF322 †	0.2759	9.23 × 10−28
cg07464716	5	157079520	SOX30 *	0.1973	4.18 × 10−8
cg02757970	2	25427350	POMC †	0.1659	7.02 × 10−10
cg02567750	1	247569605	NLRP3 †	0.1654	3.67 × 10−12
cg21834645	6	10409467	TFAP2A *	0.1625	2.03 × 10−9
cg14640588	6	158028202	ZDHHC14 *	−0.2891	1.91 × 10−8
cg08486986	16	1512134	CLCN7 *	−0.2090	6.16 × 10−8
cg02023276	7	151049248	NUB1 *	−0.1973	1.43 × 10−11
cg20095851	12	6486701	SCNN1A *	−0.1743	1.01 × 10−13
cg21627955	13	46952839	KIAA0226L *	−0.1621	1.02 × 10−8

Definition of abbreviations: CLCN7 = chloride voltage-gated channel 7; KIAA0226L = Rubicon like autophagy enhancer; NLRP3 = NLR family pyrin domain containing 3; NUB1 = negative regulator of ubiquitin like proteins 1; POMC = proopiomelanocortin; SCNN1A = sodium channel epithelial 1 subunit α; SOX30 = SRY-box transcription factor 30; TFAP2A = transcription factor AP-2 α; ZDHHC14 = zinc finger DHHC-type palmitoyltransferase 14.

Adjusted for sex, age, cell type, inhaled corticosteroid use, pack-years, and smoking intensity.

Illumina annotation based on the UCSC database.

Nearest gene as reported by GREAT version 4.0.4.

COPD-associated DMPs were annotated to 1,089 genes. Using the ImmGen RNA-seq data set (20), we identified a total of 5,328 genes expressed exclusively by nonmacrophage cell types (i.e., genes expressed by innate lymphocytes, αβ and γδ T cells, B cells, stem cells, dendritic cells, granulocytes, mast cells, basophils, eosinophils, and stromal cells but not by macrophages; see the data supplement for details). The overlap between these two gene sets was 8.4%.

Using DMRcate, we identified 7,097 DMRs. The most significant DMR (false discovery rate–adjusted P = 2.58 × 10−43) was a region on chr 15 spanning five CpG sites and overlapping the promotor of ZNF609 (zinc finger protein 609) (see Table E3).

Enrichment analyses was conducted separately for DMPs (annotated to 1,089 genes) and DMRs (annotated to 5,491 genes); results below are those overlapping in DMP and DMR analyses (for complete results, see Tables E4–E7). Gene ontology analysis showed significant enrichment for 502 terms, including transcription factor activity (molecular functions), nervous system and anatomical structure development, and regulation of metabolic processes (biological processes). KEGG pathway analyses revealed seven significant terms, many related to various cancers but also to cAMP and Rap1 signaling pathways as well as different types of cell junctions.

Epigenetic Aging

Using the DNA methylation age calculator of Horvath (18), we calculated the subjects’ DNA methylation age. As expected, there was a strong correlation between DNA methylation and chronological age (r = 0.71; see Figure E4), but no association was found between COPD status and the difference in DNA methylation age and chronological age (see the data supplement), nor were COPD-associated DMPs enriched in the CpGs included in Horvath’s epigenetic clock (data not shown).

Correlation of DNA Methylation and Gene Expression

Using a data set with paired BAL cell DNA methylation and gene expression data from Ringh and colleagues (19) (see also Table E10), we assessed the correlation between degrees of gene expression and CpG methylation. Of the 1,155 COPD-associated DMPs in the present study, 1,065 were available for analysis in Ringh and colleagues’ data set (within a ±250-kb window of a transcript). Of these, 79 unique DMPs (101 total CpG–transcript associations) correlated significantly with gene expression (Bonferroni P < 0.05; see Table E9). Fifty-four percent of all correlations were negative, with the lowest P value as well as the largest effect estimate (P = 4.81 × 10−5; ρ = −0.80) reported for cg18196647 and CPD (carboxypeptidase D) (Figure 3). Among positive correlations, the one between cg13267718 and FLI1 (friend leukemia integration 1) transcription factor had both the lowest P value and the largest effect estimate (P = 1.24 × 10−5; ρ = 0.83; Figure 3).

Figure 3.

Correlations between methylation and predicted gene expression levels for selected pairs. (A) CPD = carboxypeptidase D. (B) FLI1 = friend leukemia integration 1. We created residuals for mRNA expression and DNA methylation and then used Pearson correlation to evaluate correlations between expression residuals and methylation residuals. These residual models were adjusted for age, sex, and smoking status.

To investigate colocalization of genetic and epigenetic variation, DMPs in the present study were mapped to the 82 SNPs associated with COPD in a recent GWAS meta-analysis of 35,735 cases and 222,076 control subjects (4). We found that 447 of the 1,155 DMPs were within a ±500-kb window of one or more COPD-associated SNPs (see Table). Of these, 116 DMPs were within 100 kb of one or more SNPs. Ten colocalized SNPs and DMPs were annotated to the same gene, and of these, ADGRG6/GPR126 (adhesion G protein–coupled receptor G6) was the only gene with more than one DMP annotated to it (Table 4). Conversely, 80 of 82 COPD-associated SNPs had one or more DMPs located within a ±500-kb window. Sixty-four SNPs were within 100 kb of one or more DMPs.

Table 4.

Differentially Methylated Positions and SNPs Annotated to the Same Gene

CpG ID	Genomic Coordinate	rs ID	Genomic Coordinate	Chr	Annotated to Gene	Base Pair Difference
cg15497249	239485905	rs11579382	239737706	1	CHRM3 †	251801
cg07386720	11676054	rs2442776	11599127	3	VGLL4 *	−76927
cg13029151	29685706	rs13073544	29430921	3	RBMS3 *	−254785
cg15257376	141121677	rs7650602	141428572	3	ZBTB38 *	306895
cg24836671	148029901	rs10037493	148475407	5	HTR4 *	445506
cg14755119	19753466	rs9350191	19842430	6	ID4 †	88964
cg18128350	142706817	rs9399401	142347764	6	ADGRG6/GPR126*	−359053
cg18131185	142623588	rs9399401	142347764	6	ADGRG6/GPR126 *	−275824
cg04461802	142623433	rs9399401	142347764	6	ADGRG6/GPR126 *	−275669
cg06280150	18328109	rs9617650	18006117	22	MICAL3 *	−321992

Definition of abbreviations: ADGRG6/GPR126 = adhesion G protein–coupled receptor G6; CHRM3 = cholinergic receptor muscarinic 3; HTR4 = 5-hydroxytryptamine receptor 4; ID4 = inhibitor of DNA binding 4, HLH protein; MICAL3 = microtubule associated monooxygenase, calponin and LIM domain containing 3; RBMS3 = RNA binding motif single stranded interacting protein 3; rs = reference SNP; VGLL4 = vestigial like family member 4; ZBTB38 = zinc finger and BTB domain containing 38.

Illumina annotation based on the UCSC database.

Nearest gene as reported by GREAT version 4.0.4.

Discussion

Differential Methylation

To the best of our knowledge, this is the first EWAS of COPD on BAL cells. Measurement of more than 850,000 CpGs revealed genome-wide differential methylation in COPD. Effect sizes were large for many DMPs, comparable to those previously reported for tobacco smoking (22). As for effect direction, roughly as many DMPs had higher as lower mean methylation in subjects with COPD relative to control subjects. This is in contrast to smoking, where a large majority of DMPs are hypomethylated (19). One advantage of using BAL fluid is that it collects lung immune cells from the most distal airways, which also are the parts most affected in COPD (12). Our results may therefore provide new insights about COPD pathogenesis, although it should be noted that distal lung epithelial cells, which are also believed to play an important role, are not collected using BAL.

Among genes annotated to the 10 DMPs with the largest effect sizes, four have previously been associated with COPD and/or disease severity: NLRP3 (NLR family pyrin domain containing 3), SOX30 (SRY-box transcription factor 30), POMC (proopiomelanocortin), and ZNF322. Of these, NLRP3 (23) and SOX30 (4) have also been put forward as candidate drug targets. NLRP3 encodes for a cytoplasmic pattern recognition receptor that, as part of the NLRP3 inflammasome, is a key regulator of the inflammatory response. In mouse models, NLRP3 has been shown to be essential for the development of COPD (23), and there is increasing evidence of an important role for the NLRP3 inflammasome in COPD pathogenesis also in humans (24, 25). SOX30 is a transcription factor previously associated with progression of emphysema in patients with COPD (26). It was recently also identified as a potential causal/effector gene in an integrated analysis of COPD GWAS, gene expression, and gene regulation data (4). In the present study, the DMP annotated to SOX30 (cg07464716) correlates significantly with the expression of U7 snRNA-associated Sm-like protein LSm10, a protein previously found to be involved in chronic mucus hypersecretion in COPD (27).

The predominant BAL fluid cell type in our data set was macrophages (≈90%), and it has been shown previously that alveolar macrophage function can be modified epigenetically in COPD (28). To investigate if differential methylation in the present study was related to macrophages, we explored the ImmGen platform and available cell-specific expression profiles. First, we compiled a list of genes expressed exclusively by cell types other than macrophages (i.e., genes expressed by innate lymphocytes, αβ and γδ T cells, B cells, stem cells, dendritic cells, granulocytes, mast cells, basophils, eosinophils, and stromal cells, but not by macrophages); we then compared this list with our list of differentially methylated COPD genes. The overlap between the nonmacrophage list from ImmGen and our COPD gene list was only 8.4%, indicating that the vast majority of identified changes likely originate from macrophages (and/or unmeasured cell types that share expression profile with macrophages).

COPD, Tobacco Smoking, or Both?

Tobacco smoking has been shown to have a clear and widespread effect on the methylome, not only in peripheral blood but also in BAL cells (19, 29). As the subjects with COPD in the present study are ex- or current smokers, one challenge was to separate the signal of smoking from that of COPD. We approached this in two ways: 1) by comparing smokers with COPD to smokers without COPD (the major characteristic separating groups was thus only disease status, not smoking and disease status as in studies in which never-smokers are used as control subjects); and 2) through statistical adjustments for effects of smoking, both cumulative (pack-years) and acute (smoking intensity; i.e., number of cigarettes per day at the time of bronchoscopy). In a recent systematic review, smoking‐related changes in DNA methylation were deemed most documented for four genes: AHRR (aryl-hydrocarbon receptor repressor), F2RL3 (F2R like thrombin), DAPK (death associated protein kinase 1), and p16/CDKN2A (cyclin dependent kinase inhibitor 2A) (30). Notably, none of these was among the 1,089 genes annotated to DMPs in the present study, possibly indicating a COPD-specific signal.

Besides smoking, we also adjusted for use of inhaled corticosteroids. Epigenetic mechanisms, including DNA methylation, have long been believed to play a central role in mediating the antiinflammatory effect of glucocorticoids (31). In COPD, systemic corticosteroid use has been associated with differential methylation (32), and recently van Nijnatten and colleagues presented data indicating that this might be the case also for inhaled corticosteroids (33).

Enrichment for Biological Processes and Pathways

KEGG and gene ontology enrichment analysis identified several biologically plausible pathways and functional gene relationships. In the latter category, we found significant enrichment for known transcription factor activity, replicating the finding of one previous EWAS of lung tissue from patients with COPD (8). Of significant pathways, many were related to malignancies, which is not unexpected, as COPD is a known risk factor, independently of smoking, for several types of cancer (34). Interestingly, the two significant pathways not evidently cancer related, cAMP and Rap1 signaling, are both targeted by drugs commonly used in COPD. β-Agonists increase cAMP concentrations through activation of adenylate cyclase, which initiates a signal transduction pathway important for bronchodilation (35) and might also affect airway epithelial barrier function through the EPAC (exchange protein directly activated by cAMP)/Rap1 pathway (36). Anticholinergics enable cAMP-mediated relaxation of smooth muscle cells by inhibiting M2 muscarinic receptors (37). Corticosteroids (38) and roflumilast (39) both stimulate cAMP production, and the latter has also been suggested to exert its effect via actions on neutrophil Rap1 (39).

Epigenetic Clock

Many features of normal aging are prominent in COPD lungs (e.g., increased cellular senescence, oxidative stress, stem cell exhaustion [40]), and it has been proposed that COPD is a condition of accelerated aging (41). DNA methylation age, calculated using Horvath’s epigenetic clock, has been validated to correlate with chronological age, not only in peripheral blood but also in BAL cells from healthy subjects (42). As epigenetic alteration is one of the hallmarks of aging (43), we hypothesized that subjects with COPD would have an increased difference between DNA methylation age and chronological age compared with control subjects. However, no significant difference was found, suggesting that if epigenetic alteration is indeed involved in accelerated aging of BAL cells, it is likely through mechanisms other than DNA methylation. Histone deacetylase enzyme activity, for example, regulates histone acetylation and has previously been associated with disease severity in COPD. Interestingly, one of the transcription factors identified in the ontology analysis in the current study, sirtuin 1, is a histone deacetylase associated not only with COPD but also with aging (44).

Integration with Gene Expression

To identify potential functional effects of differential methylation, we investigated correlation with gene expression in a separate BAL cell data set. We found many DMPs for which methylation was strongly correlated with gene expression, indicating that these CpGs may have a direct functional effect on disease manifestation.

Among the 10 genes with the highest correlations between DMP methylation and gene expression, 7 have previously been associated with COPD and/or lung function (FLI1, C11orf45 [chromosome 11 open reading frame 45], ASF1A [anti-silencing function 1A histone chaperone], GOSR1 [Golgi SNAP receptor complex member 1], NSRP1 [nuclear speckle splicing regulatory protein 1], SHH [sonic hedgehog signaling molecule], and PCDHAC1 [protocadherin α subfamily C, 1]) and 4 with lung cancer (FLI1, ASF1A, CPD, and SHH), whereas PCDHB3 (protocadherin β 3) and PGAM2 (phosphoglycerate mutase 2) have not previously been linked to lung disease or lung function.

The most significant positive correlation was found for cg13267718 and FLI1. This protooncogene encodes for a transcription factor believed to play a role in both small-cell lung cancer (45) and non–small-cell lung cancer (46), both comorbidities of COPD. FLI1 has also been associated with lung function (47) and, specifically in COPD, has been found to be upregulated in activated neutrophils (48). Dysregulated neutrophilic inflammation has long been considered a prominent feature of COPD, contributing to tissue destruction and the emblematic progressive decline in lung function (49).

One major question in COPD research is why some, but not all, long-time smokers develop COPD. Genetic determinants have been proposed to partially explain this variation in disease susceptibility, and in a large-scale GWAS meta-analysis, Sakornsakolpat and colleagues recently identified 82 loci associated with the disease (4). In the present study, 38.7% of DMPs were colocalized with one or more of these COPD-associated SNPs, indicating that these might be sites where cis-acting DNA motifs and/or DNA-binding factors shape DNA methylation and, conversely, that the remaining 61.3% of DMPs are more likely influenced by environmental determinants. Correlation with gene expression was found among DMPs in both these groups, suggesting that in COPD, DNA methylation might at some sites serve as a mediator between genetics and disease and at others be acquired alterations with functional relevance.

Among genes annotated to sites of colocalized epigenetic and genetic variation, ADGRG6/GPR126 stood out as the only one with more than one DMP annotated to it. This gene was differentially methylated in COPD also in one previous EWAS on epithelial cells from bronchial brushings (10). In GWASs and gene expression studies, ADGRG6/GPR126 has been associated not only with COPD (50) but also independently with DlCO/Va (51). Technically, DlCO/Va is an index of how efficiently carbon monoxide is transferred over the alveolar–capillary membrane. In practice, DlCO/Va is often low in subjects with emphysema, and it is used together with other measurements to distinguish among different pulmonary pathologies (52). ADGRG6/GPR126 codes for a G protein–coupled receptor that modulates the VEGF (vascular endothelial growth factor) cascade in hypoxia-induced angiogenesis (53), which in turn might affect DlCO/Va through decreased perfusion. In summary, findings from multiple tissues and several omics layers indicate a possible role for ADGRG6/GPR126 in COPD pathogenesis that might be mediated through epigenetic mechanisms.

Strengths and Limitations

One strength of the present study is its well-characterized study population. All subjects were recruited from the longitudinal and population-based OLIN COPD study, which provided granular and reliable data on smoking habits, use of inhaled corticosteroids, and other possible confounders (13). In the comprehensive statistical analysis, these data were then used as covariates to sift out a signal as COPD specific as possible. Another strength is that we investigated cells from the lungs, not from peripheral blood. Even though COPD is a systemic disease, when trying to uncover its underlying epigenetic mechanisms, investigating the main target organ is likely to produce more relevant data than using surrogate tissue. In addition, to account for potential cell type effects, we adjusted our models using measured cell count.

One limitation of the present study is that it included relatively few participants. Although this is true for most studies dependent on bronchoscopies, a larger study population would have been desirable. Also, although we did detect widespread differential methylation associated with COPD, previous studies of lung diseases other than COPD indicate that a larger sample size might be needed to identify a stable and replicable signature (54). As this is the first published EWAS on COPD BAL cells, 1) one major limitation was that replication in an independent COPD BAL data set was not possible; and 2) correlation between BAL cell DNA methylation and gene expression had to be assessed in a non-COPD data set, and results are thus dependent on the assumption that methylation–expression patterns are fairly stable among individuals and groups. Although replication in peripheral blood data sets would have been possible, transcriptomic studies of COPD have shown that gene expression patterns differ between target organ and surrogate tissue such as peripheral blood (11). And because DNA methylation is one of the main mechanisms regulating tissue-specific gene expression (55), it seems reasonable to assume that DNA methylation patterns would also differ among compartments. Thus, results of a replication analysis in a peripheral blood data set, whether positive or negative, would be hard to interpret with any certainty. Finally, a cross-sectional study such as the present one cannot determine the causal direction of an association. In other words, it cannot conclude whether the identified epigenetic variation contributed to COPD susceptibility, or if it was a result of the disease, postdisease processes, or even disease-associated drug interventions.

Conclusions

We found epigenome-wide differential methylation in COPD lung cells from BAL. For many DMPs, effect sizes were large, and some also displayed a strong correlation with gene expression. Almost half of DMPs were located in close proximity to sites where genetic variation (SNPs) is known to be associated with COPD, sites where DNA sequence and/or DNA-binding factors might dictate methylation status. Conversely, remaining DMPs are more likely influenced by external factors such as tobacco smoke and air pollutants. Our interpretation is that genetic and acquired epigenetic factors both might affect COPD susceptibility. As this is the first ever EWAS of COPD BAL cells, the results need to be validated in future studies.