A genome-wide association study of bronchodilator response in participants of European and African ancestry from six independent cohorts.

Introduction: Bronchodilator response (BDR) is a measurement of acute bronchodilation in response to short-acting β2-agonists, with a heritability between 10 and 40%. Identifying genetic variants associated with BDR may lead to a better understanding of its complex pathophysiology.
Methods: We performed a genome-wide association study (GWAS) of BDR in six adult cohorts with participants of European ancestry (EA) and African ancestry (AA) including community cohorts and cohorts ascertained on the basis of obstructive pulmonary disease. Validation analysis was carried out in two paediatric asthma cohorts.
Results: A total of 10 623 EA and 3597 AA participants were included in the analyses. No single nucleotide polymorphism (SNP) was associated with BDR at the conventional genome-wide significance threshold (p<5×10-8). Performing fine mapping and using a threshold of p<5×10-6 to identify suggestive variants of interest, we identified three SNPs with possible biological relevance: rs35870000 (within FREM1), which may be involved in IgE- and IL5-induced changes in airway smooth muscle cell responsiveness; rs10426116 (within ZNF284), a zinc finger protein, which has been implicated in asthma and BDR previously; and rs4782614 (near ATP2C2), involved in calcium transmembrane transport. Validation in paediatric cohorts yielded no significant SNPs, possibly due to age-genotype interaction effects.
Conclusion: Ancestry-stratified and ancestry-combined GWAS meta-analyses of over 14 000 participants did not identify genetic variants associated with BDR at the genome-wide significance threshold, although a less stringent threshold identified three variants showing suggestive evidence of association. A common definition and protocol for measuring BDR in research may improve future efforts to identify variants associated with BDR.

Introduction

Chronic lower respiratory diseases, which include asthma and COPD, impose a major global health burden, affecting >530 million people [1]. The first-line treatment for asthma and COPD patients includes targeting airway hyperresponsiveness by use of β2-agonists which act on β2-adrenergic receptors located on airway smooth muscle (ASM) cells in the lower respiratory tract, resulting in muscle relaxation and bronchodilation [2]. While measurement of airway hyperresponsiveness is occasionally used as a diagnostic tool for asthma, measuring acute bronchodilation in response to short-acting β2-agonists (SABAs) is more widely used to differentiate fixed from reversible airflow obstruction and determine response to treatment.

Acute bronchodilation in response to pharmacological bronchodilator, referred to as bronchodilator response (BDR), is often defined as the per cent change in forced expiratory volume in 1 s (FEV1) or forced vital capacity (FVC) occurring ∼15 min after inhalation of a SABA. In the general population, most individuals have little response to medication with 95% of asymptomatic never-smoker adults experiencing an increase in FEV1 of 9% or less after SABA use [3]. However, this response is variable and can be influenced by factors such as presence of symptoms [4], and smoking and anthropometric measurements, which are thought to explain 7–16% of the variability in BDR [5]. At least a portion of this variability appears to be influenced by genetics and has been cited to range between 10 and 40% [6, 7]. Determining the genetic variants associated with BDR may lead to a better understanding of the pathophysiology underlying this important measure of airway disease and response to treatment.

Candidate genes reported to be associated with BDR include the β2-adrenergic receptor (ADRB2) [8], adenylyl cyclase type 9 (ADCY9) [9], corticotrophin-releasing hormone receptor 2 (CRHR2) [10] and arginase 1 (ARG1) [11]. Genome-wide association studies (GWASs) and a whole-genome sequencing pharmacogenetic study of BDR in asthma and COPD participants have identified genetic variants from a number of genes, highlighting the different biological pathways likely involved in controlling BDR [12-19]. These studies suggest that BDR is a complex trait related to multiple mechanisms involving more than the β2-adrenergic receptors alone.

However, the reported associations identified thus far by GWASs have shown limited replication across studies. While this lack of replication may be related to differences in study population characteristics (such as age, race/ethnicity, lung pathology) or gene–environment interactions [6, 20], the range of reported associations indicate the complexity of this trait and the need for further studies to replicate these previous findings or identify new genetic variants associated with BDR mechanistic pathways.

By including a large number of cohorts with participants of both African and European ancestry, and participants from asthma and COPD cohorts as well as participants from a community-based study which did not recruit participants based on disease or smoking status, the aim of these ancestry-stratified and ancestry-combined GWAS meta-analyses of BDR was to identify and replicate genetic variants associated with BDR. We sought to further support our findings through validation analyses that were undertaken in two paediatric asthma cohorts of European ancestry and Hispanic/Latino ethnicity.

Methods

Study population

Six non-overlapping cohorts contributed to our primary analysis: the Framingham Heart Study (FHS) [21]; Genetic Epidemiology of COPD (COPDGene) [22]; Evaluation of COPD Longitudinally to Identify Predictive Surrogate End-points (ECLIPSE) [23]; Treatment of Emphysema with a Gamma-Selective Retinoid Agonist (TESRA) [24]; the Bergen, Norway COPD Cohort, Genetics of Chronic Obstructive Lung Disease (GenKOLS) [25]; and the Subpopulations and Intermediate Outcome Measures in COPD Study (SPIROMICS) [26]. To assess whether associations observed among adults are also present among children, we used two paediatric asthma cohorts – the Genetics of Asthma in Costa Rica Study (GACRS) [27] and the Childhood Asthma Management Program (CAMP) [28]. All cohorts have participants who self-identify as being of European ancestry (EA), and two cohorts (COPDGene and SPIROMICS) also have participants who self-identify as being of African ancestry (AA). All study protocols were approved by the respective local Institutional Review Boards, and written informed consent for genetic studies was obtained from all participants or their parents, in the case of paediatric participants.

Spirometry and bronchodilator response

Standardised spirometry was performed in all study cohorts according to the 2005 European Respiratory Society/American Thoracic Society guidelines [29]. BDR was measured as a per cent change in FEV1 and FVC following the administration of 180 μg (two inhalations) of albuterol for participants in FHS and COPDGene and 400 μg of salbutamol for participants in ECLIPSE and GenKOLS. For participants in SPIROMICS, BDR was measured following the administration of 360 μg (four inhalations) of albuterol and 68 μg (four inhalations) of ipratropium bromide.

Genotyping and quality control

Different genotyping platforms were used across these cohorts, details of which can be found in supplementary table S1. Genotyped SNPs were excluded if they had a high missingness rate (cohort-specific thresholds ranged between 1 and 5%), significant departure from Hardy–Weinberg equilibrium, high discordance rates between duplicated samples, and if they were monomorphic in the sample or had low minor allele frequency (MAF<1%). Imputation was performed using the Haplotype Reference Consortium (HRC) [30] reference panel or 1000 Genomes Project [31]. For each cohort, SNPs with poor imputation quality (imputation quality r-squared ≤0.3) were excluded from further analysis. Quantile-quantile (Q-Q) plots are provided in supplementary figure S1.

Statistical analysis

Two BDR phenotypes (Per cent change of FEV1 and Per cent change of FVC) were defined as our primary outcomes:

Per cent change of FEV1=100×(Post_FEV1 − Pre_FEV1)/Pre_FEV1

Per cent change of FVC=100×(Post_FVC − Pre_FVC)/Pre_FVC

The analyses were adjusted for age at baseline (year), sex, smoking status, mean pre-bronchodilator and post-bronchodilator FEV1 or FVC, and cohort-specific significant principal components of ancestry to account for population stratification. We included participants who have both pre- and post-bronchodilator measures for FEV1 and FVC available and had no missing covariates, excluding outliers defined as participants with an FEV1 or FVC bronchodilator response falling outside of 3 SD from the mean value for each phenotype. We also excluded participants with per cent change of FEV1 >100% or per cent change of FVC >100%.

Association analyses of SNPs were performed in each cohort using multiple linear regression models adjusting for covariates, for each outcome. The primary meta-analysis was stratified by ancestry group to minimise the effects of heterogeneity; this analysis is referred to as the stratified meta-analysis. A secondary analysis was performed by combining the EA and AA participants; the results of this secondary analysis are reported separately and referred to as the combined meta-analysis. We used a linear mixed effects model to account for familial relatedness in the FHS cohort, using expected kinship coefficient within family. Ancestry-specific inverse-variance-weighted fixed-effects meta-analysis was performed using METAL software [32]. We considered SNPs with p-value <5×10−8 as genome-wide significant [33]. We performed validation analyses on two separate paediatric asthma cohorts, of different genetic ancestry, using the same analysis strategy.

We then performed a linkage disequilibrium clumping procedure to identify distinct signals using the software PLINK. The top associated SNPs (p-value <5×10−6) located in close physical proximity (±1 Mb) were clumped by linkage disequilibrium measured by r2 >0.05, and only the strongest trait-associated SNP within each linkage disequilibrium block was kept, using CEU and YRI individuals from the 1000 Genomes (phase 1 version 3) as reference panels for EA and AA groups, respectively [34].

Validation in paediatric cohort

Validation analysis was carried out on EA participants in the CAMP study and Hispanic/Latino participants in the GACRS study using the SNPs identified from the discovery analysis (SNPs with association p<5×10−6 and MAF≥5%). Bonferroni-adjusted p-value thresholds were used: p<0.05/131=0.0004 for the stratified meta-analysis, and p<0.05/57=0.0009 for the combined meta-analysis.

Replication of previously reported BDR genetic variants

We then searched the literature for SNPs that have reached genome-wide significance (p<5×10−8) in GWASs of BDR. Other BDR GWASs include SHARP (which combined data from CAMP, CARE and ACRN), SAGE I, SAGE II, GALA II and COPDGene. Participants in these published studies included children and adults of African, European and Hispanic/Latino ancestry with asthma or COPD. Because one of the cohorts in our meta-analysis was the COPDGene cohort, we repeated our analysis excluding all participants from this cohort for comparison. We also evaluated SNPs previously reported to be associated with asthma from a large multiancestry GWAS meta-analysis, as well as SNPs associated with COPD from several GWASs (supplementary table S4; includes references). A Bonferroni-adjusted p-value threshold of 0.0003 (0.05/143) was used to account for the total number of SNPs from the literature that were evaluated (n=143).

Fine mapping using PAINTOR

We used ancestry-specific linkage disequilibrium reference panels to conduct fine mapping analyses to help identify likely causal variants within 250-kb flanking each of the strongest trait-associated SNP identified above, using PAINTOR (Probabilistic Annotation INTegratOR) software [35]. PAINTOR is a software that integrates the strength of association signals from a test of association considering both linkage disequilibrium structure and functional annotation information to calculate a posterior probability for each SNP as the source of the association signal in a region. We used PAINTOR utilities to obtain the linkage disequilibrium correlation matrix for each locus using the 1000 Genomes phase 3 as the reference panels for EA and AA analyses separately [31]. We considered the following functional annotation categories from the annotation library provided by PAINTOR [36]: 1) lung-related annotations (E096); 2) fetal lung-related annotations (E88); 3) fetal lung fibroblast cell line-related annotations (E017); 4) lung fibroblast primary cell-related annotations (E128); and 5) smooth muscle-related annotations (E076, E078, E111). We first ran PAINTOR on each functional annotation separately and prioritised annotations based on the improvement in the fit of the model. Since none of the SNPs were labelled with any of the specified functional annotations, we ran a final PAINTOR model without any annotations to calculate the posterior probability for each SNP within our fine mapping loci.

Results

Descriptive characteristics of study cohorts

The baseline characteristics of all cohorts and participants included in the analyses are shown in table 1. A total of 10 623 participants of EA and 3597 participants of AA were included in the association analyses. TESRA differed from the remainder of the cohorts in that there were no individuals diagnosed with asthma in this cohort, and all reported smokers had a prior history of smoking but were not active smokers at the time of enrolment into the study. GenKOLS and FHS participants had the lowest reported mean pack-years of smoking; by design the FHS participants included never-smokers, while GenKOLS only enrolled current or former smokers. A description of paediatric cohorts can be found in the online supplementary material.

TABLE 1

Characteristics of the six cohorts included in the genome-wide association meta-analyses of bronchodilator response and the two paediatric cohorts used for validation analyses

	Discovery stage in adult cohorts								Validation in paediatric cohorts
	European ancestry						African ancestry		European ancestry	Hispanic/Latino
	FHS	COPDGene	ECLIPSE	GenKOLS	TESRA	SPIROMICS	COPDGene	SPIROMICS	CAMP	GACRS
Subjects n	1003	6591	1428	873	377	1779	3179	418	570	979
Age years (mean±sd)	57.9±14.0	65.1±8.8	62.7±7.5	59.7±10.3	66.6±8.0	65.2±8.2	54.7±7.2	58.5±8.8	8.92±2.13	9.19±1.88
Male n (%)	486 (48)	3442 (52)	940 (66)	507 (58)	261 (69)	993 (55.3)	1774 (56)	205 (48.8)	340 (60)	570 (58)
BMI kg·m−2 (mean±sd)	27.7±5.5	28.7±6.0	26.9±5.4	26.0±4.5	25.9±4.9	27.8±5.1	29.1±6.7	27.9±6.0	NA	NA
Current smoker n (%)	120 (12)	2582 (39)	490 (34)	398 (46)	0 (0)	598 (33.3)	2544 (80)	272 (64.8)	NA	NA
Former smoker n (%)	511 (51)	4009 (61)	938 (66)	475 (54)	377 (100)	1779 (100)§	635 (20)	418 (100)§	NA	NA
Pack-years (mean±sd)	25.3±22.7	47.3±26.0	47.4±27.1	25.3±17.1	45.9±23.7	52.2±29.0	38.3±21.4	41.2±17.4	NA	NA
Lung disease n (%)
Asthma	131 (13)	1030 (16)	244 (17)	247 (28)	NA+	300 (17.5)	681 (21)	130 (32.6)	NA	NA
COPD#	244 (24)	2789 (42)	1272 (89)	443 (51)	377 (100)	982 (54.7)	811 (26)	183 (43.6)	NA	NA
Pulmonary function at baseline
Baseline FEV1 as % predicted¶ (mean±sd)	87.4±15.3	70.0±26.0	51.2±23.5	69.5±25.3	46.30±9.3	65.2±26.0	79.0±24.1	70.0±27.1	95.0±13.2	99.7±16.5
Baseline FVC as % predicted (mean±sd)	101.1±15.6	82.9±18.8	83.5±21.8	87.0±17.3	88.8±17.3	84.2±19.2	87.9±19.2	87.0±21.1	105.4±12.6	105.4±16.0
Baseline FEV1/FVC as % predicted (mean±sd)	86.0±8.0	82.5±20.6	63.6±18.9	77.6±19.1	53.2±10.8	75.7±2.1	88.6±16.5	79.2±1.8	90.4±8.8	94.8±8.3
Baseline FEV1/FVC as % predicted (mean±sd)	67±7	63±16	48±15	64±16	42±08	56±16	70±13	62±16	80±8	84±7
Bronchodilator response
Change in FEV1 as % of baseline (mean±sd)	5.3±5.7	6.3±9.7	9.5±11.2	5.8±9.8	7.9±10.4	13.2±12.8	4.4±10.4	10.8±11.3	9.7±8.1	4.9±7.4
Change in FVC as % of baseline (mean±sd)	−0.02±4.8	4.5±11.0	8.5±12.1	2.9±8.6	6.2±8.6	8.8±11.2	2.3±11.1	6.9±10.7	1.8±4.1	1.3±5.7

FHS: Framingham Heart Study; COPDGene: Genetic Epidemiology of COPD; ECLIPSE: Evaluation of COPD Longitudinally to Identify Predictive Surrogate End-points; TESRA: Treatment of Emphysema with a Gamma-Selective Retinoid Agonist; GenKOLS: Bergen, Norway COPD Cohort; SPIROMICS: Subpopulations and Intermediate Outcome Measures in COPD Study; BMI: body mass index; FEV1: forced expiratory volume in 1 s; FVC: forced vital capacity. #: COPD definition – FEV1/FVC <0.7 and FEV1 (% predicted) <0.8; ¶: per cent predicted values are calculated using Hankinson equations; +: TESRA cohort did not include participants with asthma; §: SPIROMICS: former smoker is defined as having at least a 20 pack-year history of smoking and inclusion criteria for SPIROMICS is a minimum of 20 pack-year smoking history.

GWAS results from meta-analysis of adult cohorts

In the primary, ancestry-stratified meta-analysis of adult participants in the FHS, COPDGene, ECLIPSE, GenKOLS, TESRA and SPIROMICS samples, no BDR–SNP association achieved genome-wide significance (p<5×10−8). We only report results of variants with a MAF of 5% or greater. A total of 131 SNPs were associated with BDR at p<5×10−6, commonly considered as “suggestive” evidence of association (supplementary table S2). SNPs with an association p<5×10−6 after linkage disequilibrium clumping are shown in table 2. In the secondary, combined meta-analysis, a total of 57 SNPs were associated with BDR at p<5×10−6 (supplementary table S2), 21 SNPs after linkage disequilibrium clumping (table 2) and one SNP (rs4782614) associated with BDR as defined by per cent change in FEV1 nearing significance (p=7.5×10−8). Manhattan plots for these analyses are shown in figure 1 (stratified meta-analysis) and figure 2 (combined meta-analysis), and the Q-Q plots in supplementary figures S1 and S2. Table S5 contains cohort-specific results.

TABLE 2

Single nucleotide polymorphisms (SNPs) associated with bronchodilator response in six independent cohort meta-analyses with p<5×10−6 after linkage disequilibrium-clumping in the ancestry-stratified meta-analysis and the combined meta-analysis

Anc/Phe#	SNP	Chr:position ¶	Allele	AF+	Effect	se	p-value	Closest gene	Distance §
Ancestry-stratified meta-analysis
EA	rs13262811	8:142425325	A	0.85	0.81	0.17	1.8E-06	PTP4A3	6162
FEV1	rs140005306	2:143707656	A	0.95	1.25	0.26	2.0E-06	KYNU	within
	rs7315471	12:13333283	A	0.07	1.14	0.24	2.1E-06	EMP1	16318
	rs10220877	15:38934425	A	0.17	0.73	0.16	3.2E-06	C15orf53	54373
	rs76890972	16:84399755	A	0.09	1	0.22	4.7E-06	ATP2C2	2373
	rs72766398	2:12489567	T	0.95	1.31	0.28	3.6E-06	MIR3681HG	within
	rs8014708	14:101214590	A	0.1	−0.94	0.2	3.8E-06	DLK1	10029
EA FVC	rs7096465	10:4250500	T	0.09	0.98	0.20	1.1E-06	LINC00702	within
	rs10894263	11:130699110	A	0.56	0.58	0.12	8.1E-07	LOC100507431	14977
	rs8046017	16:84378488	A	0.57	−0.56	0.12	4.3E-06	WFDC1	15031
	rs10020466	4:181306647	A	0.71	0.59	0.13	3.8E-06	LINC00290	678595
	rs989808	4:47166781	T	0.16	0.73	0.16	2.6E-06	GABRB1	within
AA	rs116556600	21:15154158	A	0.43	−1.26	0.25	5.8E-07	MIR8069-1	57563
FEV1	rs11130868	3:6249417	C	0.6	1.2	0.24	9.8E-07	GRM7-AS3	424627
	rs35870000	9:14801710	A	0.12	−1.72	0.37	3.4E-06	FREM1	within
	rs2532841	10:119069842	A	0.66	1.16	0.25	4.8E-06	PDZD8	within
	rs4852100	2:240571827	A	0.76	−1.34	0.29	4.9E-06	LOC150935	112726
	rs58795512	2:40312398	A	0.94	−2.54	0.53	1.6E-06	SLC8A1-AS1	within
	rs491175	2:45151690	C	0.6	1.26	0.24	2.4E-07	SIX3-AS1	15602
	rs186201422	10:133458567	T	0.063	−2.85	0.62	3.7E-06	N/A	N/A
	rs144465006	11:18273453	A	0.077	−2.35	0.51	3.5E-06	N/A	N/A
	rs10426116	19:44591553	T	0.059	2.55	0.52	1.1E-06	ZNF284	within
AA FVC	rs2400730	3:135796524	A	0.29	1.53	0.29	1.6E-07	PPP2R3A	within
	rs10495859	2:36241146	T	0.82	−1.77	0.34	2.0E-07	LOC100288911	340745
	rs12406488	1:88071848	A	0.87	−1.8	0.39	3.7E-06	LINC01364	234510
	rs11680558	2:11522391	A	0.086	2.68	0.56	1.4E-06	LINC00570	11715
	rs41265547	2:189872891	C	0.21	1.51	0.32	1.9E-06	COL3A1	within
	rs6864533	5:20534743	T	0.093	2.43	0.48	3.3E-07	CDH18	within
	rs1736601	6:167683260	T	0.64	1.24	0.27	3.5E-06	UNC93A	21542
	rs7336381	13:114963436	A	0.77	−1.42	0.31	3.4E-06	CDC16	36925
	rs80202087	15:78653983	A	0.095	−2.12	0.46	4.7E-06	CRABP1	13411
	rs1974858	16:26743681	A	0.42	−1.23	0.26	2.8E-06	C16orf82	334537
	rs3935851	8:31133902	T	0.1	2.41	0.52	3.1E-06	WRN	102625
	rs13437873	7:48189850	A	0.74	1.36	0.29	2.8E-06	ABCA13	21206
	rs2217516	8:94401023	A	0.57	1.24	0.26	1.6E-06	LINC00535	within
	rs12713765	2:72023225	A	0.59	1.27	0.25	6.4E-07	DYSF	109332
	rs113520112	2:110362336	C	0.093	2.24	0.44	3.5E-07	SEPT10	within
	rs74491400	2:224088691	T	0.062	2.49	0.53	2.4E-06	KCNE4	168334
	rs28672350	4:24166929	A	0.88	1.89	0.4	2.6E-06	PPARGC1A	275229
	rs55699229	8:19114322	T	0.9	−2.25	0.47	1.7E-06	LOC100128993	11290
	rs144212455	15:53363818	C	0.94	−2.59	0.54	1.9E-06	ONECUT1	281609
	rs199964754	18:63446314	T	0.51	1.42	0.3	2.4E-06	CDH7	within
Combined meta-analysis
FEV1	rs1953833	1:239941671	A	0.71	−0.54	0.12	3.6E-06	CHRM3	within
	rs140005306	2:143707656	A	0.95	1.20	0.26	3.0E-06	KYNU	within
	rs6780151	3:6249820	A	0.60	1.26	0.26	1.2E-06	N/A	N/A
	rs9868510	3:105700120	A	0.51	−0.52	0.11	2.1E-06	CBLB	112233
	rs9385215	6:122074667	A	0.10	−0.88	0.19	2.7E-06	GJA1	303777
	rs13262811	8:142425325	A	0.86	0.76	0.16	2.2E-06	PTP4A3	6162
	rs7079679	10:102055546	A	0.53	0.50	0.11	3.0E-06	PKD2L1	within
	rs186201422	10:133458567	T	0.06	−2.85	0.62	3.7E-06	N/A	N/A
	rs144465006	11:18273453	A	0.08	−2.35	0.51	3.5E-06	N/A	N/A
	rs11055267	12:13329855	A	0.93	−1.14	0.23	1.1E-06	EMP1	19746
	rs4782614	16:84396105	A	0.32	0.64	0.12	7.5E-08	ATP2C2	6023
	rs76890972	16:84399755	A	0.08	1.02	0.21	1.6E-06	ATP2C2	2373
	rs10426116	19:44591553	T	0.06	2.56	0.52	1.1E-06	ZNF284	within
FVC	rs113520112	2:110362336	C	0.09	2.24	0.44	3.6E-07	SEPT10	within
	rs73803631	4:24159527	A	0.91	2.08	0.45	4.5E-06	PPARGC1A	267827
	rs55699229	8:19114322	T	0.90	−2.25	0.47	1.7E-06	N/A	N/A
	rs75638212	10:22950466	T	0.86	−0.74	0.15	1.0E-06	PIP4K2A	within
	rs144212455	15:53363818	C	0.94	−2.59	0.54	1.9E-06	ONECUT1	281609
	rs4782614	16:84396105	A	0.32	0.57	0.12	2.5E-06	ATP2C2	6023
	rs71355331	18:45645902	T	0.94	1.12	0.24	3.6E-06	ZBTB7C	within
	rs199964754	18:63446314	T	0.51	1.42	0.30	2.4E-06	CDH7	within

#: ancestry/phenotype; ¶: Position: genome reference used hg19; +: AF: allele frequency; §: distance from closest gene measured in base pairs; “within” denotes SNPs that are located within the gene.

FIGURE 1

Plot of p-values by chromosome for change in forced expiratory volume in 1 s (FEV1) and forced vital capacity (FVC) in response to bronchodilator in the ancestry-stratified meta-analysis. Ancestry and phenotype: a) change in FEV1 as per cent predicted in European ancestry (EA) participants; b) change in FVC as per cent predicted in EA participants; c) change in FEV1 as per cent predicted in African ancestry (AA) participants; d) change in FVC as per cent predicted in AA participants.

FIGURE 2

Plot of p-values by chromosome for change in forced expiratory volume in 1 s (FEV1) and forced vital capacity (FVC) in response to bronchodilator in the combined meta-analysis. a) Change in FEV1 as per cent predicted in European ancestry (EA) and African ancestry participants. b) Change in FVC as per cent predicted in EA and AA participants.

Validation analyses yielded no significant associations after adjusting for multiple comparisons (p<0.05/342=0.00015 for the stratified meta-analysis, and p<0.05/57=0.0009 for the combined meta-analysis) (supplementary table S3). A test of proportions to compare the direction of effect of those SNPs identified in the validation analysis yielded p=0.005 for CAMP and p=0.57 for GACRS. Therefore, the direction of effect in the CAMP study, in which participants are EA, was consistent with the direction in our primary analysis, while that in GACRS, which is made up of Hispanic/Latino participants, was not.

Replication of BDR, COPD and asthma previously reported SNPs

None of the SNPs that have previously attained genome-wide statistical significance in published GWASs of BDR [12-19], COPD or asthma reached statistical significance in our analysis after correcting for multiple comparisons (p<0.0003) (supplementary table S4).

Fine mapping and regional plots for biologically relevant SNPs

We performed fine mapping and examined regional plots for the SNPs showing “suggestive” evidence of an association (i.e. p<5×10−6) in the ancestry-stratified meta-analysis to identify the most relevant signals. We examined SNPs with a posterior probability above 70% after fine mapping to narrow our focus on variants of potential interest (table 3). A total of seven SNPs with MAF ≥5% and a BDR-association p<5×10−6 had a posterior probability >70%. For these seven SNPs, the number of SNPs included in the 95% credible sets ranged from 2 to 8. We examined the genes that are most closely associated with these SNPs for possible biological relevance in BDR. Regional plots for the two genes encoding proteins in biological pathways plausibly related to BDR based on fine mapping results among participants of AA (FREM1 and ZNF284) are shown in figure 3. There were no significant recombination peaks between FREM1 and rs35870000 or between ZNF284 and rs10426116, indicating the signal for association obtained from these two SNPs in our study may be related to these two genes. The regional plot for the SNP showing the strongest association in the combined meta-analysis, rs4782614 (p=7.5×10−8), is shown in figure 4.

TABLE 3

Single nucleotide polymorphisms (SNPs) with p<5×10−6 and posterior probability >70% after fine mapping using PAINTOR (Probabilistic Annotation INTegratOR) software

Ancestry/phenotype	SNP	Chr:position#	Allele	p-value	Posterior probability	95% credible set	Closest gene	Distance ¶
EA FEV1	rs13262811	8:142425325	A	1.8×10−6	0.78	6	PTP4A3	6162
AA FEV1	rs35870000	9:14801710	A	3.4×10−6	0.75	3	FREM1	within
	rs144465006	11:18273453	A	3.5×10−6	0.79	8	N/A	N/A
	rs10426116	19:44591553	T	1.1×10−6	0.73	6	ZNF284	within
AA FVC	rs11680558	2:11522391	A	1.4×10−6	0.79	4	LINC00570	11715
	rs55699229	8:19114322	T	1.7×10−6	0.93	2	LOC100128993	11290
	rs3935851	8:31133902	T	3.1×10−6	0.88	4	WRN	102625

#: position: genome reference used hg19; ¶: distance from closest gene measured in base pairs; “within” denotes SNPs that are located within the gene.

FIGURE 3

Regional and recombination plots near a) FREM1 and b) ZNF284.

FIGURE 4

Regional and recombination plots near ATP2C2.

Expression quantitative trait loci in lung and whole blood

Next, we examined whether “suggestive” SNPs from our study with p<5×10−6 had previously been reported as significant expression quantitative trait loci (eQTLs) in lung tissue and whole blood using data from the Genotype Tissue Expression (GTEx) Project. A total of 10 SNPs from the stratified meta-analysis and four SNPs from the combined meta-analysis matched previously reported eQTLs from GTEx samples of lung and whole blood (table 4). The rs35870000 SNP is an eQTL for RP11-408A13.4 in whole blood, although there is no known clinical significance for this gene. The rs4782614 SNP is located near the ATP2C2 gene and is an eQTL for ATP2C2 in whole blood.

TABLE 4

Single nucleotide polymorphisms (SNPs) from ancestry-stratified and combined meta-analyses with p<5×10–6 and corresponding expression quantitative trait loci (eQTLs) in lung tissue and whole blood from the Genotype Tissue Expression (GTEx) Project

Ancestry/phenotype	SNP	Nearest gene	GWAS p-value	Lung eQTL	WB eQTL
Ancestry-stratified meta-analysis
EA FEV1	rs13262811	PTP4A3	1.8E-06	N/A	SLC45A4
EA FVC	rs7096465	LINC00702	1.1E-06	LINC00702	N/A
	rs7079903	LINC00702	1.1E-06	LINC00702	N/A
	rs11252396	LINC00702	1.4E-06	LINC00702	N/A
	rs7087175	LINC00702	2.9E-06	LINC00702	N/A
AA FEV1	rs35870000	FREM1	3.4E-06	N/A	RP11-408A13.4
	rs2532841	PDZD8	4.8E-06	PDZD8	PDZD8
AA FVC	rs7336381	CDC16	3.4E-06	CDC16, UPF3A	CDC16
	rs1736600	UNC93A	3.6E-06	RP11-568A7.3	N/A
	rs1757120	UNC93A	4.2E-06	RP11-568A7.3	N/A
Combined meta-analysis
FEV1	rs7079679	PKD2L1	3.0E-06	CWF19L1, BLOC1S2	CWF19L1
	rs13262811	PTP4A3	2.2E-06	N/A	SLC45A4
	rs4782614	ATP2C2	7.5E-08	N/A	ATP2C2
FVC	rs4782614	ATP2C2	2.5E-06	N/A	ATP2C2

GWAS: genome-wide association study; WB: whole blood.

Discussion

In this GWAS meta-analysis of BDR no SNP reached genome-wide significance in association with either per cent change in FEV1 or per cent change in FVC. This study has a number of strengths that represent an important contribution to the literature on the genetics of BDR which include using a large number of independent cohorts, having participants of both AA and EA in the primary analysis, and of Hispanic/Latino ancestry in the validation analysis, and finally, the inclusion of participants with and without obstructive lung disease.

With a total of 10 623 participants of EA and 3597 participants of AA, this study had 80% power to detect SNPs explaining 0.4% and 1.1% of the variance in the outcome for our stratified EA and AA analyses, respectively, at a genome-wide significance level of 5×10−8. Therefore, we had adequate power to detect SNPs with small effects among EA participants and SNPs with larger effects among AA participants. However, despite having adequate power, we were unable to replicate findings of previous GWASs, even at a higher p-value cut-off (p<5×10−6). Although this is a negative finding, it is nonetheless an important one given the size of this meta-GWAS.

Unlike GWASs for disease phenotypes, such as asthma, or spirometric measurements of pulmonary function, such as the FEV1/FVC ratio, genetic variants reported to be associated with BDR have not been well replicated across studies. This lack of replication is probably due to a number of factors including the physiological complexity of the BDR phenotype and variation in the definition of BDR across studies, as well as variations in medication protocols used in different studies, between-study differences in race/ethnicity and pulmonary pathology, and finally, environmental factors affecting BDR. For example, there is evidence that at least one environmental factor, inhaled corticosteroid treatment, can modulate the effect of a gene on BDR response [20]. Environmental factors that increase airway inflammation or epithelial injury may also modify the effects of genetic variants on observed BDR, but further gene–environment interaction studies of BDR are needed to identify such interactions. Unfortunately, gene–environment interaction studies require a very large number of participants before adequate statistical power is achieved.

The findings of our study and previous GWAS studies of BDR reinforce the notion that the molecular pathways influencing BDR are complex and involve multiple genes. While no SNP reached genome-wide significance in our analyses, several genes lie in close proximity to SNPs with p-values that suggested association (p<5×10−6) and may have biological functions that could influence BDR. Using fine mapping and regional plots to hone in on signals of potential biological importance, we highlight three SNPs that, to our knowledge, have not been previously reported as being associated with BDR: rs35870000 and rs10426116, both of which had suggestive associations with BDR as defined by per cent change in FEV1 (p=3.4×10−6 and p=1.1×10−6 respectively) only among participants of AA in the stratified meta-analysis, and rs4782614, which had a suggestive association with BDR as defined by per cent change in FEV1 (p=7.53×10−8) and FVC (p=2.53×10−6) in the combined meta-analysis.

The rs35870000 SNP is a missense variant located within the gene FREM1 (FRAS1 related ECM1), with no other known genes within 60 kbp of this SNP. The Toll-like/interleukin (IL)-1 receptor regulator (TILRR) which regulates the IL-1 receptor type I (IL1R1) is encoded by a spliced variant of this FREM1gene [37]. Inhibition of IL1R1 has been shown to decrease IL-5- and IgE-induced changes in ASM responsiveness [38]. It is therefore conceivable that ASM contraction, and by extension BDR, can be modulated by FREM1 through the action of cytokines widely implicated in asthma pathophysiology. Notably, the association between rs35870000 and BDR was seen in the AA subgroup in our study, a population in which anti-IL-5 therapy has been reported to be most effective [39].

The rs10426116 SNP lies within the ZNF284 (zinc finger protein 284) gene on chromosome 19 and is a 3′ UTR variant (the 3′UTR region often contains regulatory regions that influence gene expression). In a gene–environment study of BDR, a region on chromosome 19 containing multiple zinc finger proteins was identified as showing an association with BDR that was modified by inhaled corticosteroids [20]. As pointed out by the authors of that study, and corroborated by the findings of the current study, further exploration of the mechanisms by which zinc finger proteins affect BDR may be worth pursuing.

The rs4782614 SNP lies within 6023 bp from the ATP2C2 (ATPase secretory pathway Ca2+ transportin 2) gene and is a significant eQTL for ATP2C2 in whole blood (as reported in GTEx). The ATP2C2 gene encodes the calcium-ATPase SPCA2 protein that is highly expressed in the Golgi apparatus and plays an important role in calcium homeostasis [40], an essential component of muscle contraction. Furthermore, intronic SNPs near ATP2C2 have been associated with frequent COPD exacerbations in a COPD GWAS [41]. Notably, and perhaps not surprisingly, other GWASs have reported enrichment of ion channels in association with BDR [16].

A major limitation of our findings is the lack of genome-wide significant associations despite our large sample size attained by pooling several studies. Although SNPs with nominally significant (p<5×10−6) associations with BDR are highlighted above, none reached conventional genome-wide significance, despite combined samples sizes of 10 623 EA and 3597 AA participants. In addition, our validation cohorts were fundamentally different from the discovery cohorts in that they were paediatric studies of asthma and did not include participants with AA. Associations observed in adults might not be present in children, especially given the possible presence of age-by-genotype interaction that has been previously reported [42]. Heterogeneity across studies included in the meta-analysis and in the paediatric cohorts may have limited our ability to detect genome-wide associations. One potential source of such heterogeneity is the difference in the medication protocols that were used to obtain the BDR in the various studies. The ascertainment of race and ethnicity by self-report leads to uncertainty regarding admixture of genetic ancestries and may be confounded by aspects of social identity, but it is unlikely to be a limitation in this study as we controlled for ancestry by including principal components as covariates in our model. Finally, another limitation relevant to all GWASs, including this one, is that the studied SNPs are often in linkage disequilibrium with the causal variants, and further experimental validation will be required to move from statistical association to mechanistic inferences about causation.

Future studies should consider including gene–environment interactions when studying BDR, such as the effect of first-hand and second-hand smoking or the use of inhaled corticosteroids, as this may help tease out the most important association signals related to BDR. Consideration of alternative, and unified, methods of defining BDR for research purposes should also be evaluated.

Please note: supplementary material is not edited by the Editorial Office, and is uploaded as it has been supplied by the author.

Supplementary material 00484-2021.SUPPLEMENT