A genome-wide analysis of the response to inhaled β2-agonists in chronic obstructive pulmonary disease.

Short-acting β2-agonist bronchodilators are the most common medications used in treating chronic obstructive pulmonary disease (COPD). Genetic variants determining bronchodilator responsiveness (BDR) in COPD have not been identified. We performed a genome-wide association study (GWAS) of BDR in 5789 current or former smokers with COPD in one African-American and four white populations. BDR was defined as the quantitative spirometric response to inhaled β2-agonists. We combined results in a meta-analysis. In the meta-analysis, single-nucleotide polymorphisms (SNPs) in the genes KCNK1 (P=2.02 × 10(-7)) and KCNJ2 (P=1.79 × 10(-7)) were the top associations with BDR. Among African Americans, SNPs in CDH13 were significantly associated with BDR (P=5.1 × 10(-9)). A nominal association with CDH13 was identified in a gene-based analysis in all subjects. We identified suggestive association with BDR among COPD subjects for variants near two potassium channel genes (KCNK1 and KCNJ2). SNPs in CDH13 were significantly associated with BDR in African Americans.The Pharmacogenomics Journal advance online publication, 27 October 2015; doi:10.1038/tpj.2015.65.

Introduction

Chronic obstructive pulmonary disease (COPD) is a disorder characterized by progressive loss of lung function. It is currently the third leading cause of death world-wide, and the global burden of disease is expected to continue to rise(1). Although cigarette smoke is the greatest risk factor for COPD, recent studies have identified several genetic risk factors for this disease(2).

Inhaled bronchodilators, including β2-agonists, play a key role in COPD management guidelines. These medications act on smooth muscle receptors in bronchial airways to produce muscle relaxation and airway dilation, resulting in improved airflow through the lungs (1), and have been shown to alleviate COPD symptoms(3). The response to inhaled bronchodilators is measured by a change in the forced expiratory volume in one second (FEV1) using standardized spirometry before and after the administration of β2-agonists. Although COPD is characterized by relatively fixed airflow limitation, up to two-thirds of COPD patients will exhibit a positive response to an inhaled bronchodilator at any one time(4).

The quantitative response to inhaled β2-agonists is a heritable trait(5), and candidate gene studies have identified several genes suggestive of association with quantitative measures of bronchodilator responsiveness (BDR)(6, 7). In addition, candidate gene(8) and genome-wide association studies (GWAS) have identified variants associated with BDR in subjects with asthma (9–11). We hypothesized that genome-wide association studies would identify associations with BDR in COPD.

Subjects and Methods

Study subjects

Details of the COPDGene, ECLIPSE, GenKOLS, and NETT studies, including study procedures, genotyping, and quality control, have been reported(12–16). COPDGene subjects were current and former smoking non-Hispanic white (NHW) or African American (AA) from the U.S. ECLIPSE subjects were Caucasian current or former smokers from Europe, North American and New Zealand. GenKOLS subjects were current and former smokers from Norway. NETT subjects were white former smokers from the U.S. All subjects had moderate to severe COPD (GOLD stage 2 or greater(17)). Subjects were excluded if they had a recent COPD exacerbation.

Spirometry

All subjects completed a respiratory questionnaire and performed standardized spirometry according to American Thoracic Society or European Respiratory Society guidelines. COPDGene, NETT, and GenKOLS subjects were tested before and approximately 20 minutes after administration of 2 puffs (180 μg) of inhaled β2-agonist (albuterol/salbutamol). ECLIPSE subjects were tested before and 15 minutes after inhalation of 400 μg β2-agonist (albuterol/salbutamol).

BDR was measured using three quantitative metrics that have been previously reported(5). BDRABS, the absolute difference in pre- versus post- bronchodilator FEV1; BDRPRED, the absolute difference in pre- versus post-bronchodilator FEV1 as a percentage of FEV1 percent predicted; and BDRBASE, the absolute difference in pre versus post bronchodilator FEV1 as a percentage of baseline FEV1.

Genotyping

All subjects were genotyped using Illumina platforms (Human Hap550 for ECLIPSE and GenKOLS, Quad610 for NETT, and Human OmniExpress for COPDGene) as previously described(13, 15). We included all variants and subjects that passed quality control, based on cluster plots (genotyped) and imputation quality (R2 ≥ 0.80) for imputed SNPs, Hardy-Weinberg equilibrium (P-value), and missingness (% threshold). Imputation was performed using MaCH and minimac with 1000 Genomes phase I v 3 European reference panels for white subjects. Cosmopolitan reference panels were used for COPDGene AA subjects. Variants with a minor allele frequency (MAF) < 1% and R2 ≤ 0.80 were excluded from analysis. Ancestry-based principal components were generated for each study using EIGENSOFT2.0(18). We performed Taqman genotyping (Applied Biosystems, Foster City, CA) for the SNPs rs114132812 and rs115067260 among 23 and 38 African American COPDGene subjects respectively, who were imputed to be carriers of the minor allele.

Statistical analysis

Baseline subject demographics and outcome variables were analyzed in R (v2.15.1). We excluded 20 subjects with BDR variables more than six standard deviations from the mean. We performed linear regression analysis for the three BDR variables in PLINK(19) including genotyped and imputed SNPs, adjusting for age, gender, pack-years smoking history, and ancestry-based principal components. We combined the results from all five samples in a fixed-effects meta-analysis using METAL(20). We additionally performed a gene-based test of significance among the top 20% of all SNPs using VEGAS (a Versatile Gene-based test for Genome-wide Association Studies). This method performs gene-based association testing by assigning SNPs within 50kb of a gene in accordance with the hg18 assembly and then uses simulation to account for linkage disequilibrium. All genes are tested for association with the trait of interest(21). Using the top 20% of significantly associated SNPs, we used VEGAS software to test 13,675 genes. Based on this number, we established a Bonferroni significance threshold of 3.6×10−6. We analyzed the top genes from our GWAS as well as the top genes identified through the VEGAS analysis using the functional annotation tool, DAVID(22, 23).

We tested the association of seven asthma and three COPD SNPs previously associated with BDR in asthma GWAS (9, 11) and COPD candidate gene studies (ADRB2, EPHX1, and SERPINE2) (6) with BDR in our meta-analysis results. We additionally tested the association of two SNPs from the β2-agonist receptor gene, ADRB2 (Arg16Gly, rs1042713 and Gln27Glu, rs1042714). We tested our top SNPs for their association with lung function (FEV1/FVC and FEV1) in the four COPD populations with the broadest range of lung function values: COPDGene NHW and AA, GenKOLs, ECLIPSE. We examined our top BDR variants for their association with BDR in two asthma GWAS(9, 11).

Results

The demographic data for each study population are presented in Table 1. These outcomes appear to follow a normal distribution (Supplementary Figure 1). All three outcomes are significantly correlated in COPDGene NHW, but BDRABS and BDRPRED appear more correlated (R2=0.97), compared to BDRABS and BDRBASE (R2=0.85) or BDRBASE and BDRPRED (R2=0.88).

Table 1

Study Subjects

	COPD Gene NHW	COPD Gene AA	ECLIPSE	GenKOLS	NETT
N	2792	811	1757	853	364
Population origin	US	US	Europe/US/New Zealand	Norway	US
Age	64.7 (8.2)	59.0 (8.2)	63.7 (7.1)	65.6 (10.0)	67.5 (5.8)
Males (%)	55.6	54.7	66.9	60.0	64.8
Pack-years	56.2 (27.8)	42.4 (23.1)	50.3 (27.4)	32.1 (18.5)	66.1 (30.9)
Current smoker (%)	34.7	60.8	35.4	46.9	0
BMI (kg/m2)	28.1 (6.1)	28.0 (6.8)	26.7 (5.6)	25.3 (4.9)	25.0 (3.5)
FEV1 pp (%)	49.7 (18.0)	52.3 (17.7)	47.6 (15.6)	50.7 (17.3)	28.0 (7.4)
FEV1/FVC	0.49 (0.13)	0.53 (0.12)	0.45 (0.12)	0.51 (0.13)	0.32 (0.06)
BDRABS (L)	0.10 (0.15)	0.07 (0.18)	0.12 (0.14)	0.10 (0.14)	0.09 (0.08)
BDRBASE (%)	9.2 (12.4)	7.6 (14.2)	10.7 (13.4)	8.1 (11.7)	13.4 (12.0)
BDRPRED (%)	3.6 (4.9)	2.89 (6.5)	3.9 (4.8)	3.4 (4.7)	3.3 (3.0)

Data are presented as mean (sd) or N (%). NHW: non-Hispanic white; AA: African American; FEV1pp: FEV1 percent predicted; BMI: Body mass index. BDRABS: (Post-BD FEV1 – Pre-BD FEV1); BDRBASE: (BDRABS/(Pre-BD FEV1*100)); BDRPRED: (BDRABS/(FEV1pp*100)).

The BDR outcomes were statistically different when comparing the four Caucasian populations (P<0.05 for ANOVA for all three BDR outcomes). However, these small differences are unlikely to be clinically significant. We additionally compared the NHW and AA subjects from COPDGene. The AA subjects demonstrated significantly lower bronchodilator responsiveness for all three outcomes, and these results remained significant, albeit small, after adjusting for age, gender, and pack-years (BDRABS: 0.10 L vs 0.07 L, P = 0.004; BDRPRED: 3.59% v 2.80%, P=0.02; BDRBASE: 9.24% vs 7.60%, P=0.007).

Table 2 presents the most significant SNPs from the GWAS in the non-Hispanic white subjects from COPDGene. The top SNP annotated to a gene is presented in this table. The full list of SNPs with P < 1×10−6 are presented in Supplementary Tables 1–3. SNP rs17575208, located upstream from the gene EPHA7 on chromosome 6, was significantly associated with BDRABS (β=0.11, P=8.29×10−9). This variant was also associated with BDRPRED (β=3.22, P=1.03×10−7) and BDRBASE (β=7.06, p=5.64×10−6), though the P-values were not genome-wide significant.

Table 2

Most significant variants and genes from the GWAS of COPDGene non-Hispanic white subjects

SNP	Chr	Nearest Gene	Beta (se)	Allele	FRQ	P value
BDRABS
rs17575208	6	EPHA7	0.11 (0.02)	A	0.99	8.92×10−9
rs7778219*	7	LOC285889	−0.03 (0.01)	G	0.11	1.62×10−7
rs2367245	17	KCNJ2	0.02 (0.01)	G	0.57	1.33×10−6
rs115310518	5	CDH18	0.064 (0.014)	T	0.98	2.20×10−6
rs78008396	6	PARK2	−0.075 (0.016)	A	0.98	2.92×10−6
rs7932838*	11	SOX6	0.05 (0.01)	G	0.03	3.34×10−6
rs11775549	8	NKX2-6	−0.044 (0.0095)	G	0.94	3.39×10−6
rs34342951	8	NKX6-3	−0.021 (0.005)	G	0.72	3.42×10−6
rs11260945	1	IGSF21	−0.046 (0.01)	G	0.95	3.74×10−6
rs73671623	8	STC1	−0.044 (0.010)	G	0.94	4.01×10−6
rs56323342	11	ME3	−0.040 (0.009)	G	0.95	4.21×10−6
rs56010187	16	HS3ST4	−0.023 (0.005)	G	0.78	5.65×10−6
rs1335517	14	C14ORF37	0.027 (0.006)	T	0.85	5.85×10−6
rs181350634	3	IFT80	0.09 (0.02)	A	0.99	6.29×10−6
rs7552783*	1	KCNK1	−0.0197 (0.004)	G	0.28	7.10×10−6
BDRBASE
rs1032243	3	MIR548A3	−1.96 (0.40)	A	0.77	9.59×10−7
rs326981	5	MTRR	−2.48 (0.51)	G	0.89	1.32×10−6
rs4772755	13	LINC00460	−1.75 (0.36)	G	0.70	1.36×10−6
rs6943859*	7	KLHDC10	3.0 (0.62)	C	0.07	1.84×10−6
rs10242432	7	ZC3HC1	−2.98 (0.62)	G	0.93	1.85×10−6
rs8032265*	15	C15orf60	3.63 (0.78)	A	0.05	3.18×10−6
rs2367245	17	KCNJ2	1.57 (0.34)	G	0.57	3.34×10−6
rs8108918*	19	VAV1	1.74 (0.37)	T	0.26	3.98×10−6
rs1004790	2	LOC645949	−5.77(1.26)	A	0.98	4.66×10−6
rs17574208	6	EPHA7	7.06 (1.53)	A	0.99	5.64×10−6
BDRPRED
rs17575208	6	EPHA7	3.22 (0.60)	A	0.99	1.03×10−7
rs7778219*	7	LOC285889	−1.03 (0.20)	G	0.11	3.90×10−7
rs11775549	8	NKX2-6	−1.53 (0.31)	G	0.94	6.99×10−6
rs73671623	8	STC1	−1.53 (0.31)	G	0.94	1.27×10−6
rs2367245	17	KCNJ2	0.629 (0.13)	G	0.57	2.53×10−6
rs78008396	6	PARK2	−2.50 (0.53)	A	0.98	3.07×10−6
rs115310518	5	CDH18	2.07 (0.45)	T	0.98	4.63×10−6
rs1032243	3	MIR548A3	−0.73 (0.16)	A	0.77	4.30×10−6
rs1335517	14	C14orf37	0.89 (0.20)	T	0.85	5.65×10−6

Presenting the top SNP from each top gene with P < 5×10−6; Adjusted for age, sex, pack-years and principal components; Results filtered for MAF > 1%, R2> 0.80.

Genotyped SNP.

Chr: Chromosome; FRQ: reference allele frequency; Allele: reference allele. The full list of SNPs with P < 5×10−6 is presented in the supplementary material.

Table 3 presents the top SNPs annotated to genes having P < 5×10−6 for the COPDGene AA subjects. The full list of SNPs with P < 5×10−6 are presented in Supplementary Tables 4–6. Variants in the gene CDH13 were significantly associated with BDRABS (rs115067260; β=0.17±0.03, P=5.05×10−9) and BDRPRED (rs114132812, β=7.63±1.32, P=1.19×10−8), and showed suggestive association with BDRBASE (rs77347308 β=−17.14±3.39, P=5.35×10−7). In addition, a genotyped SNP in the gene SGCD was significantly associated with BDRBASE (rs10056066, β=7.12±1.29, P=4.86×10−8), and several rare imputed SNPs in the gene GOLGA8B were associated with the outcome BDRPRED (rs76677753, β=9.49±1.67, P=1.9×10−8). A recent GWAS using COPD subjects from the Lung Health Study population identified the variants in the gene SGCD as associated with airway responsiveness measured by methacholine challenge test in a physiologically distinct asthma cohort(24). While the response to inhaled methacholine is distinct from the response to inhaled bronchodilators, similar mechanisms of smooth muscle activation could be involved. We tested the two top SGCD SNPs from that study in our AA population. The SNP rs456290 was associated with BDRBLINE (P=0.02), and the SNP rs2642660 approached replication (P=0.08). These SNPs were not associated with BDRBLINE in the NHW population, or in the meta-analysis.

Table 3

Most significant variants and genes from the GWAS of COPDGene African American subjects

SNP	Chr	Nearest Gene	Beta (se)	Allele	FRQ	P value
BDRABS
rs115067260	16	CDH13	0.17 (0.03)	A	0.97	5.05×10−9
rs140948272	17	PITPNA	0.33 (0.06)	C	0.99	5.24×10−9
rs76677753	15	GOLGA8B	0.25 (0.05)	A	0.99	7.96×10−8
rs78060357	7	PLXNA4	0.08 (0.02)	A	0.90	8.06×10−8
rs145442019	16	HSBP1	0.203 (0.038)	A	0.99	8.58×10−8
rs13345720*	19	RNASEH2A	−0.05 (0.010)	C	0.34	1.28×10−7
rs114871691	12	BTBD11	0.114 (0.02)	A	0.95	1.32×10−7
rs60085550	1	OR6F1	0.157 (0.030)	G	0.97	1.76×10−7
rs2295599*	6	SYCP2L	0.07 (0.014)	A	0.10	5.31×10−7
rs149447163	2	UBR3	0.191 (0.038)	G	0.98	5.62×10−7
rs145986148	16	KCNG4	0.169 (0.034)	A	0.98	6.59×10−7
rs12529809	6	ELOVL2	−0.07 (0.01)	G	0.88	8.09×10−7
rs387092	16	MLYCD	0.20 (0.04)	A	0.99	9.90×10−7
rs4742936	9	ABCA1	0.088 (0.018)	C	0.94	9.95×10−7
BDRBASE
rs10056066*	5	SGCD*	7.12 (1.29)	A	0.07	4.86×10−8
rs143376495	6	TFAP2B	−19.47 (3.63)	G	0.99	1.06×10−7
rs2295599*	6	SYCP2L*	6.0 (1.13)	A	0.10	1.38×10−7
rs11651753*	17	MIR4315-2	4.22 (0815)	T	0.22	2.81×10−7
rs75762663	13	PRR20E	−19.22 (3.80)	A	0.99	5.34×10−7
rs77347308	16	CDH13	−17.14 (3.39)	A	0.99	5.35×10−7
rs77094751	18	SALL3	−10.23 (2.03)	G	0.96	5.65×10−7
rs62385074	5	LCP2	−10.17 (2.03)	G	0.96	6.32×10−7
rs141697850	4	TRAM1L1	−15.05 (3.01)	C	0.99	7.02×10−7
rs12529809	6	ELOVL2	−5.674 (1.14)	G	0.88	7.85×10−7
rs2575515	4	GRXCR1	3.525 (0.71)	T	0.47	8.33×10−7
rs7848679	9	ABCA1	6.96 (1.41)	G	0.93	9.31×10−7
rs76999017	2	INPP5D	−8.263 (1.67)	G	0.95	9.70×10−7
BDRPRED
rs114132812	16	CDH13	7.63 (1.32)	A	0.98	1.19×10−8
rs76677753	15	GOLGA8B	9.49 (1.67)	A	0.99	1.90×10−8
rs13345720*	19	RNASEH2A	−1.83 (0.35)	C	0.34	1.60×10−7
rs114871691	12	BTBD11	4.116 (0.79)	A	0.95	2.06×10−7
rs72653721	6	SYCP2L	−2.746 (0.53)	C	0.88	2.91×10−7
rs147388556	1	LPHN2	4.17 (0.81)	T	0.96	3.32×10−7
rs11245387	10	FAM175B	−1.76 (0.34)	T	0.40	3.49×10−7
rs145442019	16	HSBP1	7.10 (1.39)	A	0.99	3.55×10−7
rs143859190	2	GALNT13	8.10 (1.58)	T	0.99	3.67×10−7
rs188637550	8	MTUS	6.89 (1.36)	T	0.98	4.60×10−7
rs78060357	7	PLXNA4	2.88 (0.57)	A	0.90	5.08×10−7
rs12529809	6	ELOVL2	−2.64 (0.52)	G	0.88	5.37×10−7
rs188524320	8	PCM1	7.01 (1.40)	C	0.99	6.39×10−7
rs4742936	9	ABCA1	3.28 (0.66)	C	0.94	7.17×10−7
rs17116498	8	FGL1	7.07 (1.42)	C	0.99	7.32×10−7
rs139218005	8	ATP6V1H	4.935 (0.89)	A	0.97	9.74×10−7

Presenting the top SNP from each top gene with P < 1×10−6; Adjusted for age, sex, pack-years and principal components; Results filtered for MAF > 1%, R2> 0.80.

Genotyped SNP.

Chr: Chromosome; FRQ: reference allele frequency; Allele: reference allele. The full list of SNPs with P < 5×10−6 is presented in the supplementary material.

The CDH13 SNPs were in LD (rs115067260 and rs114132812, R2=0.60, D’=1.0). These variants were monomorphic in the Caucasian populations. We tested additional variants within CDH13 in the COPDGene NHW subjects for association with BDR. Rs4783331 was nominally associated with BDRABS (β=0.11±0.02, P=9.39×10−5), and five additional SNPs in this gene had P<0.001.

We tested our model assumptions of normal distributions of the BDR traits, focusing on BDRABS in the African Americans. The BDR traits appeared to fit a normal distribution (Supplementary figure 1). We additionally examined the residuals from linear regression for the BDR outcomes in the African American subjects, adjusted for age, pack-years, and gender, which appeared normally distributed and were consistent with our model assumption (Supplementary figure 2). We performed an inverse normal transformation of the BDRABS phenotype, and tested this trait in a GWAS. The CDH13 SNPs, rs115067260 (P=4.46×10−6) and rs114132812 (P=5.10×10−6), remained the top associations. Taqman genotyping verified the imputation accuracy of the CDH13 SNPs. 22/23 imputed carriers of the SNP rs114132812 were heterozygous for the minor allele, and for rs115067260, 35/37 subjects were verified as heterozygous and one subject homozygous for the minor allele.

We performed a meta-analysis of the results of all five study populations (Table 4, Figure 1). A SNP in the potassium channel, subfamily K, Member 1 gene (KCNK1) demonstrated suggestive association with the outcome BDRABS (β =−0.0142, P=2.02×10−7). A SNP upstream from the gene KCNJ2 (rs9898686) showed suggestive association with all three traits: BDRABS, (β =0.014, P=2.05×10−7), BDRBASE (β=1.26, P=1.83×10−7) and BDRPRED (β =0.44, P=1.22×10−6). Variants in the gene MC5R (melanocortin 5 receptor) were suggestively associated with BDRPRED (rs12956045, β= −0.45, P=4.69×10−7). Several other variants in the KCNJ2 region also demonstrated nominal association, and these variants were in linkage disequilibrium (R2>0.80). SNPs upstream from KCNJ2 were recently identified as associated with lung function (measured by FEV1) in a joint meta-analysis of SNP and SNP-by-smoking effects in a population-based sample(25). We examined the top SNP from that analysis, rs11654749. This SNP was nominally associated with BDRABS in COPDGene NHW subjects, with the same effect direction (β=0.011, P=0.007). We performed the meta-analysis without the AA population (Supplementary Table1), with very similar results as that including all five study populations together.

Table 4

Meta-analysis of top variants associated with BDR from 5 COPD populations

SNP	Chr	Gene	Allele	Effect	Std Err	P Value
BDRABS
rs61824320	1	KCNK1	A	−0.0142	0.0027	2.02×10−7
rs9898686	17	KCNJ2	T	0.0138	0.0027	2.05×10−7
rs68120025	17	KCNJ2	C	0.0138	0.0027	2.05×10−7
rs68193035	17	KCNJ2	A	0.0138	0.0027	2.05×10−7
rs9906150	17	KCNJ2	A	0.0137	0.0026	2.06×10−7
rs9899756	17	KCNJ2	C	0.0138	0.0027	2.08×10−7
rs67158616	17	KCNJ2	C	0.014	0.0028	5.35×10−7
rs12956045	18	MC5R	A	−0.013	0.0026	8.30×10−7
BDRPRED
rs68193035	17	KCNJ2	A	1.26	0.24	1.79×10−7
rs68120025	17	KCNJ2	C	1.26	0.24	1.79×10−7
rs9906150	17	KCNJ2	A	1.25	0.24	1.82×10−7
rs9898686	17	KCNJ2	T	1.26	0.24	1.83×10−7
rs9899756	17	KCNJ2	C	1.25	0.24	1.85×10−7
rs11871999	17	KCNJ2	T	1.26	0.26	9.36×10−7
BDRBASE
rs12956045	18	MC5R	A	−0.45	0.09	4.69×10−7
rs2278843	10	PKD2L1	G	0.42	0.08	7.11×10−7
rs7079679	10	PKD2L1	A	0.42	0.08	7.26×10−7
rs28650403	19	NWD1	G	−1.26	0.26	9.82×10−7
rs28428860	19	NWD1	C	−1.26	0.26	9.94×10−7
rs28668077	19	NWD1	C	−1.26	0.26	9.95×10−7

Most significant variants associated with BDR among the meta-analysis including COPDGene NHW and AA, ECLIPSE, GenKOLS, and NETT. Allele: reference allele.

Figure 1

Manhattan plots for meta-analysis results for each bronchodilator responsiveness outcome

We examined SNPs previously associated with BDR in asthma GWAS (9, 11, 26, 27) (Supplementary Table 9). In asthma, the SNP rs4452682 in the gene SLC22A23 was associated with BDR(11). This SNP was nominally associated with BDRBASE in the COPD meta-analysis (β =0.63, P=2.5×10−3), although the effect size was in the opposite direction. A rare variant in the gene SLC24A4 was also previously associated with BDR in an asthma study(26). Although this SNP, rs77441273, was not present in our GWAS meta-analysis, several SNPs in this gene demonstrated nominal association with BDRABS (rs60243508, β=0.017, P=0.008). The gene SPATA13 was previously associated with BDR in a gene-based GWAS among African Americans with asthma(27). Although this gene did not replicate in our gene-based VEGAS analysis and the top reported SNP did not replicate in our meta-analysis, we tested all 464 of the genotyped and imputed SNPs available in the gene SPATA13 for their association with BDRABS. Twenty-eight of these SNPs were nominally associated with BDRABS in the COPD meta-analysis, including rs9511156 (β=0.02, P=0.007); however these associations were not significant after correction for multiple testing. In addition, among the COPDGene AA population, the top SNP rs9507294 from the asthma study showed nominal association with BDRBASE (β=1.94 , P=0.05). The remainder of the asthma GWAS SNPs were not significantly associated with BDR in the COPD analyses, and none of the asthma BDR genes were significant in the gene-based VEGAS analysis. The ADRB2 codon 16 and 27 SNPs (rs1042713 and rs1042714) that have been previously identified in asthmatic populations were not significantly associated with BDR in the COPD GWAS. None of the top COPD variants were associated with BDR when examined in two prior asthma GWAS(9, 11). The SNPs from candidate genes previously associated with BDR in COPD populations (6), including EPHX1 (rs3753661), SERPINE2 (rs6712954), and ADRB2 (rs1042717), were not associated with BDR in our analysis.

We tested the top SNPs from the BDR meta-analysis for association with lung function in four of the five COPD populations (COPDGene NHW and AA, GenKOLs, ECLIPSE). SNP rs61824320 in the gene KCNK1 was significantly associated with FEV1/FVC (β= −0.0042 , P=0.03) and FEV1 (β= −0.02 , P=0.04) (M. Cho and S. Lutz, personal communication). The remainder of top SNPs from the BDR meta-analysis were not significantly associated with lung function.

In the VEGAS gene-based analysis (Supplementary Tables 10–12), using the P value results from the meta-analysis of all five studies, KCNK1 approached genome-wide significance for the outcome BDRABS (P=8.4×10−5) and BDRPRED (P=3.8×10−4). The gene SH2B adaptor protein 3, SH2B3, was a top gene for all three traits, approaching genome wide significance for BDRABS (P=3.20×10−5), BDRPRED (7.0 × 10−5) and BDRBASE (1.50 × 10−5). We examined the top genes from the GWAS in the VEGAS analysis. KCNJ2 and SGCD were not significant in the gene-based analysis; however, CDH13 showed nominal significance for all three traits. We additionally performed Gene Ontology analysis of the top 50 genes from the GWAS and from the VEGAS analysis using the functional annotation software, DAVID. Among the top genes in the GWAS, there was enrichment for ion channel and ion transport genes, including KCNK1, FXYD1, and PKD2L1. Among the top genes from the VEGAS analysis, there was enrichment for chemotaxis and lipid biosynthesis.

Discussion

This is the first genome-wide association study of the response to inhaled β2-agonists among COPD subjects. In a meta-analysis including five COPD populations and over 5000 subjects, we identified several genetic variants associated with the response to inhaled β2-agonist bronchodilators. In the African-Americans from COPDGene, several variants in the genes CDH13 (Cadherin-13), SGCD (Sarcoglycan delta), and GOLGA8B (golgin A8 family, member B) demonstrated genome-wide significance for their association with the response to β2-agonists. In the non-Hispanic white COPDGene population, SNPs upstream and within the potassium channel genes, KCNJ2 and KCNK1 respectively, approached genome-wide significant association with the response to inhaled bronchodilators. This association remained in a meta-analysis including all five COPD case populations, although not at the genome-wide significance level.

In the primary analysis among the COPDGene NHW population, the SNP rs17575208 on chromosome 6 upstream from the gene Ephrin-type A receptor 7, EPHA7, was genome-wide significant with a P value of 7.23×10−9 for the outcome BDRABS, 8.3×10−8 for BDRPRED and 4.15×10−6 for BDRBASE. Although mutations in this gene have been found in resected non-small cell lung cancer human specimens, little is known about a potential role for EPHA7 in COPD, asthma, or BDR. In our meta-analysis, this SNP had a P value of 2.32 × 10−5 for the outcome BDRABS and 1.30×10−4 for the outcome BDRPRED, however the association was not significant among the other populations tested (P>0.05). This is an imputed SNP with minor allele frequency of 1.4% but good imputation quality (R2 = 0.82), and therefore may suggest a promising gene for future studies.

We identified several SNPs on chromosome 17 upstream from KCNJ2, also known as KIR2.1 or Potassium inwardly-rectifying channel, subfamily J, member 2, that were suggestively associated with BDR. Inwardly-rectifying potassium channels were initially described in cardiac, skeletal, and brain tissue, and especially in smooth muscle of small arterioles. These channels play a role in potassium-mediated constriction in response to hypoxemia or ischemia(28). KCNJ2 encodes Kir2.1, a strong inwardly rectifying potassium channel, which has been identified in small and large bronchial smooth muscle cells and plays a role in membrane depolarization (29) (30). Although the exact role for these channels in smooth muscle relaxation remains to be determined, it is possible that these channels play a role in the response to increased extracellular potassium, such as that induced by hypoxemia or acidic environments, leading to membrane hyperpolarization and smooth muscle relaxation(28, 30).

Genotyped and imputed SNPs within the potassium channel gene, KCNK1, located on chromosome 1, were also associated with BDRABS. This gene is also known as TWIK-1, and encodes the 2-pore protein potassium channel subfamily K member 1, or inward rectifying potassium channel protein TWIK-1. Two-pore potassium channels have been identified in a lung epithelial cell line(31). TWIK-1 transcript has been identified on the apical membrane of human bronchial epithelial cells, and has been suggested to play a role in hyperpolarization of membrane action potential(32). The top variant in this gene was additionally associated with measures of lung function, suggesting a potential link to COPD severity. Further studies will be necessary to confirm the roles of these potassium channel genes in BDR and COPD severity.

The identification of multiple potassium channel genes in the single SNP and gene-based analyses suggests a potential role for these channels in moderating the response to inhaled bronchodilators. Further, both the GWAS and gene-based analysis were enriched for ion channel genes. Cellular potassium levels play a key role in maintaining membrane potential, and potassium channels have been demonstrated to mediate the effects of β-agonists(2). Other potassium channel genes, such as the KCNQ voltage activated channels, have been found to ameliorate methacholine bronchoconstriction in rat lung models(33), and these inwardly-rectifying potassium channels may play a similar role in bronchial smooth muscle relaxation. The potassium channel opening medication cromakalin has been tested in animal and human asthma subjects(34). However, limited knowledge is available about the role of potassium channels in mediating smooth muscle relaxation. The identification of potassium-channel genes suggests the importance of revisiting this class of medications for COPD and asthma therapeutics.

We noted a statistically significant, but clinically small, difference in the response to inhaled β2-agonists between the non-Hispanic white and African American subjects in COPDGene. To our knowledge, this is the first demonstration that NHW and AA subjects with COPD may respond differently to inhaled β2-agonists. Because of this difference, we examined the AA population alone for variants associated with BDR. In AA subjects, SNPs in the genes CDH13, SGCD, and GOLGA8B were associated with bronchodilator responsiveness. CDH13 encodes the protein T-cadherin, which functions as an adiponectin receptor(35) and is expressed in mouse lungs in response to allergen stimulation with ovalbumin (36). Adiponectin is protective against allergen-induced inflammatory cell response in mouse lungs and airway hyperresponsiveness(37) and T-cadherin knock-out mice demonstrate reduced immune response, airway hyperresponsiveness, and mucus hyperplasia compared to wild-type mice (36). Elevated adiponectin levels have been associated with increased radiographic measures of percent emphysema and lower response to inhaled bronchodilators among subjects with COPD(38).

SGCD encodes the dystrophin-sarcoglycan complex protein subunit sarcoglycan-δ. This protein complex is expressed in skeletal and cardiac muscle and is thought to play a role in limb girdle muscular dystrophy(39). The delta-sarcoglycan complex has been identified in airway smooth muscle cells, and plays a role in mediating the transition of airway smooth muscle cells from contractile to proliferative phenotypes(40), suggesting a possible role in COPD pathogenesis. In addition, variants in the SGCD gene were recently associated with airway hyperresponsiveness in a GWAS among COPD subjects from the large multicenter Lung Health Study. These variants were also nominally associated with BDRBLINE in our population.

Variants in the gene GOLGA8B (Golgin A8 Family, Member B, GOLGA8B) were associated with BDRPRED and approached significant association with BDRABS. Although these were imputed SNPs with low minor allele frequency, the imputation was of good quality. GOLGA8B encodes a golgi system autoantigen, and this region has been associated with myopia in a large GWAS(41), although a potential role in bronchodilator responsiveness is unknown.

We examined the response to an inhaled bronchodilator as a quantitative variable using three measures that have been used in prior epidemiologic and genetic studies(5, 6). In a family-based study, Palmer and colleagues demonstrated that both BDRABS and BDRPRED have greater than 30% heritability, while BDRBASE is less heritable(5), suggesting that all three outcomes were suitable phenotypes to test for genetic association. The absolute change in FEV1 after β2-agonist administration is the most straightforward measure, but it does not account for baseline lung function, which is reduced in COPD. In contrast, BDRBASE has been shown to correlate with baseline lung function(42). These variables were all highly correlated. As no single measure appears to be the most comprehensive, we analyzed all three traits.

Quantitative measures of the response to inhaled bronchodilators differ from the binary definition used by the American Thoracic Society and European Respiratory Society(1). Prior studies have suggested that this binary outcome does not identify a phenotypically distinct subset of COPD patients, since the presence or absence of a bronchodilator response does not predict clinical outcomes and demonstrates intra-individual variability (42, 43). In contrast, linkage and candidate gene studies have previously identified genetic risk factors for quantitative measures of bronchodilator responsiveness(5, 6). In addition to our new GWAS results, these genetic associations suggest that there are distinct genetic risk factors that play a role in determining the quantitative response to inhaled β2-agonists. Identifying these markers may help to identify COPD patients who demonstrate a greater response to β2-agonists, or who may be unlikely to benefit and therefore should be prescribed alternative medications.

Several GWAS have identified SNPs and genes associated with BDR in asthma (8, 11, 26, 27, 44). We examined the top variants from the asthma studies for association with BDR in the COPD populations, as well as SNPs in the ADRB2 gene that have previously been associated with clinical response to long-acting β2-agonist administration in COPD(44–46). A variant in the SLC22A23 gene was nominally associated with BDR in COPD, although the effect size was in the opposite direction as originally reported. A previously identified intergenic SNP, rs11252394, on chromosome 5 also suggested significant association with BDRABS. However, none of the other asthma BDR variants demonstrated significant association. Although ADRB2 may be a candidate gene for bronchodilator responsiveness in asthma, variants from this gene have not demonstrated consistent association with BDR in COPD (6, 47, 48).

Our study has several limitations. It is common for investigators to replicate top GWAS findings in a replication population. However, in order to improve our power to detect an association, we used all available COPD cohorts in our meta-analysis to perform the largest GWAS of BDR to date. Our most significant findings are in biologically plausible genes, and the effect sizes are similar across all included cohorts. It is encouraging that we did find some cross-over between our top hits and those in asthma populations for both bronchodilator responsiveness and airway hyperresponsiveness despite the fact that these are different study populations. Although we identified several variants upstream from KCNJ2 as associated with BDRABS and BDRPRED, these results failed to meet genome-wide significance in the meta-analysis. An examination of each population demonstrates that these variants all had a similar effect size (Supplementary Figure 1, Supplementary Table 8), suggesting that the lack of significance may be related to sample size. In contrast, studies that have identified variants associated with lung function in the general population have had sample sizes up to ten times larger than the current study(49). This is the only genome-wide study of the response to inhaled bronchodilators in COPD performed to date. We were specifically interested in identifying genes associated with BDR in COPD populations, and therefore we are limited to available COPD cohorts for this analysis. Although the GWAS meta-analysis did not demonstrate genome-wide significance, the top SNP is upstream from the gene KCNJ2, and the protein-product potassium channel is relevant to the phenotype being studied. In addition, variants in this same region have previously been associated with lung function. We are additionally limited by the use of a one-time measurement of bronchodilator responsiveness. Although BDR as a binary trait is not necessarily stable in an individual COPD patient over serial measurements (43), we used quantitative outcomes in this analysis. Post-bronchodilator FEV1 is a stable phenotype over time(5), diminishing the noise in the quantitative measures. However, the presence of intra-individual variability may have diluted our ability to identify a significant genome-wide association. Future longitudinal studies in these populations that can account for intra-individual variability may better identify genetic risk factors for these outcomes.

Within the COPDGene African American subjects, there were several SNPs that demonstrated genome-wide significance, but were of low minor allele frequency, including imputed SNPs in CDH13. However, all SNPs had a minor allele frequency greater than 1%, with excellent imputation quality (R2 ≥ 90%). Although the BDR outcomes were normally distributed, minimal skewing could result in false positive associations especially among variants with low minor allele frequency. In order to test our assumption of normality, we performed inverse normal transformation of the BDR outcomes and tested for variants for association with this transformed outcome. The order of SNPs was preserved with this transformed outcome, suggesting that our assumption of normality was correct. We confirmed imputation accuracy through direct genotyping. Furthermore, animal studies provide good evidence that CDH13 is biologically plausible as a gene potentially involved in the bronchodilator pathway. Although the top variants associated with BDR in the AA population were monomorphic in the Caucasian populations, several other CDH13 variants were nominally associated with BDR. In addition, the gene-based test identified this gene as nominally associated with bronchodilator responsiveness even though the top SNPs in the AA analysis were not included in the gene-based test.

In summary, in the largest COPD pharmacogenetics GWAS to date, we demonstrated that variants upstream from the gene KCNJ2 are associated with response to an inhaled short acting β2 -agonist bronchodilator. In addition, several SNPs in the genes CDH13 and SGCD were significantly associated with BDR in African Americans. These results may point to novel assessments or potential novel therapeutic pathways for COPD. Future studies will require larger COPD populations to identify genome-wide significant variants, and functional studies will help to identify a role for the SNPs and genes highlighted in the GWAS.