ITGB5 and AGFG1 variants are associated with severity of airway responsiveness

Background Airway hyperresponsiveness (AHR), a primary characteristic of asthma, involves increased airway smooth muscle contractility in response to certain exposures. We sought to determine whether common genetic variants were associated with AHR severity. Methods A genome-wide association study (GWAS) of AHR, quantified as the natural log of the dosage of methacholine causing a 20% drop in FEV1, was performed with 994 non-Hispanic white asthmatic subjects from three drug clinical trials: CAMP, CARE, and ACRN. Genotyping was performed on Affymetrix 6.0 arrays, and imputed data based on HapMap Phase 2, was used to measure the association of SNPs with AHR using a linear regression model. Replication of primary findings was attempted in 650 white subjects from DAG, and 3,354 white subjects from LHS. Evidence that the top SNPs were eQTL of their respective genes was sought using expression data available for 419 white CAMP subjects. Results The top primary GWAS associations were in rs848788 (P-value 7.2E-07) and rs6731443 (P-value 2.5E-06), located within the ITGB5 and AGFG1 genes, respectively. The AGFG1 result replicated at a nominally significant level in one independent population (LHS P-value 0.012), and the SNP had a nominally significant unadjusted P-value (0.0067) for being an eQTL of AGFG1. Conclusions Based on current knowledge of ITGB5 and AGFG1, our results suggest that variants within these genes may be involved in modulating AHR. Future functional studies are required to confirm that our associations represent true biologically significant findings.


Background
Asthma is a common chronic respiratory disease affecting over 25 million Americans and its prevalence has risen over the past decade [1]. Airway hyperresponsiveness (AHR), in which the muscles of the airway contract in response to certain exposures, is one of the primary characteristics of asthma and one that has been correlated with current asthma severity and future lung function [2][3][4]. Bronchoprovocation challenges with methacholine or histamine have been widely used to quantify AHR in clinical and research settings [5]. These tests consist in administration of increasing dosages of a bronchoconstrictor until a specific decrease in lung function is measured (e.g., 20% drop in FEV 1 ). The mechanisms by which AHR occurs are not fully understood, but AHR has been related directly to changes in airway smooth muscle contractility as well as inflammation and airway remodeling [6,7]. Asthma is a disease with demonstrated heritability in humans, and several genes, including the IKZF3-ZPBP2-GSDMB-ORMDL3 locus, HLA-DQ, IL1RL1, IL33, TSLP, SLC22A5, SMAD3, and RORA have been consistently associated with asthma in genome-wide association studies (GWAS) [8,9]. Although AHR is often used as a quantifiable and reproducible surrogate of asthma in animal model studies of asthma genetics [10,11], the genetics of AHR have not been widely studied in humans. Nonetheless, the heritability of AHR is supported by a previous twin study [12] and positional cloning and linkage analysis studies of AHR [13,14]. Our goal was to measure the association of genetic variants with AHR severity via a GWAS.
We utilized data from 994 non-Hispanic white asthma patients to measure the association of single nucleotide polymorphisms (SNPs) with severity of AHR and found that the strongest associations were at variants in two genes: ArfGAP with FG repeats 1 (AGFG1) and integrin, beta 5 (ITGB5). After attempting to replicate primary findings in two independent populations (one composed of subjects with asthma and one composed of patients with COPD) and searching for evidence that variants modified expression levels of their respective genes, we found additional evidence to support the involvement of AGFG1 in AHR. Evidence that ITGB5 was involved in AHR was found in previous functional studies of human airway smooth muscle cells.

Subjects
The primary group of subjects consisted of 994 non-Hispanic white asthmatics from the Childhood Asthma Management Program (CAMP) [15], and subsets of clinical trials within the Childhood Asthma Research and Education (CARE) network [16], and the Asthma Clinical Research Network (ACRN) [17] participating in the NHLBI SNP Health Association Resource Asthma Resource project (SHARP). Some demographic characteristics of these cohorts are in Table 1 and further details are provided in the Additional file 1. For each trial, methacholine challenge tests were performed according to American Thoracic Society criteria [18]. Baseline AHR measures were utilized, and AHR was quantified as the natural log of the provocative concentration of methacholine that caused a 20% decrease in FEV 1 (LnPC20).

Genotyping and quality control
Genome-wide SNP genotyping for SHARP subjects was performed by Affymetrix, Inc. (Santa Clara, CA) using the Affymetrix Genome-Wide Human SNP Array 6.0. Dataset quality control (QC) included the removal of related subjects (i.e. CAMP and CARE siblings), subjects with missingness >5%, SNPs with minor allele frequency (MAF) <0.05, missingness >5%, or not passing Hardy-Weinberg equilibrium based on a threshold of 1E-03. The genomic inflation factor (λ GC ) for the remaining 442,036 SNPs genotyped in SHARP subjects was 1.007, demonstrating minimal population stratification. Further genotyping and QC measures are provided in the Additional file 1.

Statistical analysis
Imputation of all SNPs available in HapMap Phase 2 Release 22 CEU data using the Markov Chain Haplotyping software (MaCH) [19] was performed for SHARP genotype data. The primary GWAS was based on the set of 2,154,322 imputed SNPs with MAF >0.05 and a ratio of empirically observed dosage variance to the expected (binomial) dosage variance greater than 0.3, indicating good quality of imputation. The association of SNPs with LnPC20 was measured with a linear regression model using dosage data as implemented in PLINK [20].
To be sure that demographic characteristics did not have a strong influence on the association results, an additional association model was created using sex, age, height, and study as covariates. Plots of association results near specific genes were created using LocusZoom with the hg18/HapMap Phase II CEU GenomeBuild/LD Population [21]. Combined P-values with replication populations were computed using Fisher's combined probability method [22] where hypothesis tests in replication populations had one-sided alternatives, based on the direction of the association in SHARP, so that SNPs with association tests in opposite directions would not produce inappropriately small P-values. Replication of SNPs with P-values <1E-05 in the primary GWAS was attempted in two independent cohorts.

Genome-wide gene expression data
The Asthma BRIDGE project collected mRNA expression level data from 2,520 mRNA microarrays on 210 Illumina HumanHT-12 chips, among which 218 arrays were genetic control arrays and 2,302 arrays were sample arrays from 1,481 subjects. We focused on the 419 arrays (47,053 gene probes) from whole blood samples of 419 white CAMP subjects utilized in the current study. Expression data were log2 transformed and quantile normalized. The R Bioconductor GGtools package (version 3.11.27) was used to perform cis-eQTL analysis, in which chi squared tests were used to assess whether genotypes of a SNP within 50 KB from both ends of a gene were associated with the expression levels of the gene after adjusting for age and gender [23].

Replication study
Dutch Asthma GWAS (DAG). This cohort was comprised of 650 DAG subjects with doctor-diagnosed asthma and documented AHR [14,24]. All subjects had smoking history and steroid use data available at the time of the AHR test. Participants with known AHR but in remission during the test were excluded. Remission was defined as not on steroids and without 20% or greater fall in FEV 1 during the AHR test. The AHR test was conducted using histamine or methacholine as a stimulus. AHR was quantified as the difference between FEV 1 at baseline and at the dose step at which a 20% or greater FEV 1 drop was achieved, divided by the dose of stimulant used (slope). Because two protocols were used, one with a 30-second tidal breathing method and a second with a 2-minute tidal breathing phase, the AHR slopes measured with the 30-second tidal breathing method were divided by 4, in order to compensate for the 4 times greater duration of administration of stimulus. Slope values were log-transformed so that they would follow a normal distribution. To compare DAG subjects to those of the primary GWAS, PC20 was calculated in 600 of the 650 subjects that had a 20% or greater drop in FEV 1 at the highest provocation dose. Genotyping was performed using the Hapmap 317 K platform or Illumina 370 Duo Chip. Tests of AHR association were performed via linear regression, with smoking and inhaled/oral steroid use as covariates using PLINK. Lung Health Study (LHS). This cohort consisted in a subset of 3,354 individuals with COPD from the multicenter (10 centers) North American LHS [25]. The initial study population consisted of 5,887 men and women (63% male) who were smokers (aged 35-60) at enrollment with spirometric evidence of mild to moderate lung function impairment. While a history of asthma was not a criterion for study exclusion, subjects were excluded if they regularly used asthma medications, including bronchodilators and glucocorticoids. For the AHR GWAS, the subset of European American LHS participants with adequate DNA and valid methacholine measurements at baseline was utilized. Spirometry was performed following ATS criteria, and details of the methacholine challenge tests have been provided previously [26]. Samples were genotyped using the Illumina Human660W-Quad v.1_A BeadChip. Genotyped SNPs were imputed from CEU using 1000 Genomes Project data as a reference. Association of SNPs with log10 of PC20 was measured using linear regression under an additive model and while adjusting for gender, age at baseline, clinic site, log10(weight in kg), FEV1 (in L), and FEV1/FVC.
The current study was approved by the Partners Human Research Committee (Partners HealthCare, Inc., Boston, MA). Collection of data for the existing human cohort studies was approved by the Institutional Review Board of the corresponding institution(s), which ensured that all procedures followed were in accordance with the ethical standards of the responsible committee on human experimentation, as detailed in the cited works. Informed consent was obtained for all study participants.

Results
Characteristics of all subjects are provided in Table 1. The primary GWAS measured the association of SNPs with AHR in 994 non-Hispanic white subjects. The quantile-quantile (QQ) plots for both imputed and genotyped SNPs revealed that the distributions of association P-values were similar to those expected for a null distribution [ Figure 1]. The lowest P-values (<1E-05; Table 2) among the primary GWAS are in four regions, two of which correspond to genes: ArfGAP with FG repeats 1 (AGFG1) and integrin, beta 5 (ITGB5) [ Figure 2; Additional file 1: Table S1]. Comparison of primary association results to those generated with a model that adjusted for sex, age, height, and study revealed that there were no great changes in rankings or P-values, particularly for the top SNPs [ Figure 3]. Because some of the top-ranked associations could represent biologically meaningful ones, we proceeded to attempt to replicate them in two independent populations. The SNPs in/near AGFG1 replicated at a nominally significant level (i.e. P-value <0.05) in LHS, while no replication was obtained in DAG [ Table 3].   , respectively) to each SNP in the plots is denoted by colors and was computed according to hg18/HapMap Phase II CEU data. Plots were created using LocusZoom [21].

Figure 3
Comparison of P-values generated for genotyped SNPs from the unadjusted primary GWAS vs. those from a model in which sex, age, height, and study were used as covariates. The results reveal that adjustment for these covariates did not drastically alter the ranking or P-values of SNPs.
To find further support that the associations measured represent biologically relevant findings, we searched for evidence that any are eQTL for the genes in/near them among CAMP subjects who were part of the Asthma BRIDGE genomics study. The SNPs in/near AGFG1 had nominally significant unadjusted P-values for being associated with expression levels of AGFG1 in whole blood [ Table 3, Figure 4]. None of these P-values passed Benjamini-Hochberg correction criteria as part of a genome-wide search for eQTL, nor were the expression   Table 2 are similar to those in this Figure. levels of AGFG1 significantly correlated with AHR. None of the other SNPs with P-value <1E-05 in the primary GWAS had nominally significant eQTL.

Discussion
Our understanding of asthma genetics has increased in recent years, partly because of the success of large-scale multi-center GWAS [8,9]. Asthma GWAS have identified a robust set of genes that are consistently associated with the disease in diverse populations (e.g. the IKZF3-ZPBP2-GSDMB-ORMDL3 locus, TSLP, and RORA), as well as genes that seem to be race/ethnicity specific (e.g. PYHIN1).
In addition to studying asthma directly, GWAS have been carried out for related traits, including bronchodilator response [27], plasma IgE concentration [28], and response to inhaled corticosteroids [29]. Compared to asthma GWAS, the study of asthma-related traits has been challenged by the smaller sized cohorts available and the lack of uniform measures of interest across cohorts. While asthma-related trait GWAS often suffer from decreased statistical power due to their smaller sample sizes, they have some advantages over the study of asthma. Asthmarelated traits can be expressed quantitatively, and this may result in increased statistical power to detect associations by capturing phenotypes that do not require arbitrary classification thresholds [30]. Moreover, the investigation of specific intermediate phenotypes of asthma may reduce phenotypic heterogeneity and facilitate the functional validation of association results as some of these traits are more easily simulated in vitro and in animal models [27]. Because AHR is a complex phenotype, our GWAS results are limited by our ability to properly define it in subjects. In a previous study of CAMP subjects, the repeatability of PC20 measures over the 4-year length of the clinical trial was high, as demonstrated by intraclass correlation coefficients, after adjusting for age, race/ethnicity, height, family income, parental education, and symptom score, of 0.67 for subjects taking budesonide or placebo and 0.73 for subjects taking nedocromil [31]. A twin study of methacholine responsiveness found that the intrapair correlation coefficient among monozygotic twins (0.67) was significantly greater than that for dizygotic twins (0.34) [12]. Positional cloning and linkage analysis studies of AHR have identified two well-validated asthma candidate genes: ADAM33 [13] and PCDH1 [14]. Although the high intraclass correlation of AHR measures and previous twin and association studies suggest that AHR is a heritable trait, AHR varies in time within individuals and is influenced by environmental factors and medication usage. Among the baseline AHR measures used in our GWAS, some were taken after placebo washout periods of varying length, while others were taken in subjects (e.g. some from CARE) who were not necessarily off of asthma medications [ Table 1]. This difference would decrease our ability to detect significant relationships, such that less significant associations would be detected than if all subjects had had a placebo washout.
We performed our GWAS using subjects with AHR measures from as many asthma trials as possible in order to maximize our statistical power to detect associations. Nonetheless, partly because AHR is not always measured in asthma clinical trials, our sample size was limited to 994. One significant difference among the trials we utilized is that the CAMP inclusion criteria included having AHR, such that 12.5 mg/dl of methacholine or less resulted in a 20% drop in FEV 1 . Nonetheless, CAMP and CARE had similarly severe AHR (mean LnPC20 0.096 (SD 1.21) and −0.25 (SD 1.35), respectively). CAMP and CARE were also similar in terms of age, sex, and baseline FEV 1 [ Table 1]. These two pediatric populations contrast with ACRN, which was composed mostly of adults, had a greater proportion of female participants, less severe AHR (mean LnPC20 0.80 (SD 0.68)), and lower baseline FEV 1 . We investigated some of the effects of trial heterogeneity on our primary association results by performing the GWAS while adjusting for sex, age, height, and study. As Figure 3 shows, these covariates did not have a strong influence on the ranks or P-values obtained. Because performing the adjustment for covariates limited our sample size further to 989 due to missing demographic variable entries for some subjects, we opted to utilize the unadjusted results.
The replication populations were also distinct from our primary populations. DAG quantified AHR differently than all other cohorts, as the slope defined by the change in FEV 1 between the stimulant dose step at which a 20% or greater FEV 1 drop was achieved vs. baseline, divided by the dose of stimulant used. This definition allowed for the inclusion of subjects who did not achieve a PC20 at the doses of stimulants administered, and hence, allowed for the inclusion of subjects with less severe AHR. Fifty of 650 DAG subjects did not have PC20 measures, and hence, the estimate of mean LnPC20 provided in Table 1 is biased in that it does not include the measures of patients with less severe AHR. The AHR severity among LHS participants was also lower, on average, than that for primary cohort subjects (mean LnPC20 1.58 (SD 0.93)). While no participants of CAMP, CARE, or ACRN were smokers at the time of baseline AHR measures, LHS was a clinical trial that specifically enrolled adult smokers with COPD, while DAG was an observational study that included both smokers and non-smokers with asthma. Because smoking is known to affect lung function, AHR measures in both of these cohorts are likely affected by the smoking status of their participants. Further, LHS was not composed of asthma patients, and hence, it is likely that some of the biological mechanisms underlying its subjects AHR do not overlap with those of asthma patients. Despite the large differences among cohorts, we attempted to replicate our primary findings in LHS and DAG to assess their generalizability.
For the primary GWAS, four independent regions had P-values <1E-05 [ Table 2]. While due to the small sample size of our primary GWAS, it is not surprising that none of these regions met traditional genome-wide significance levels, two of the top regions were composed of multiple SNPs within two genes (i.e. AGFG1 and ITGB5), which increased the likelihood that the associations were biologically significant. The top ITGB5 association was at SNP rs848788 (P-value 7.2E-07), while the top AGFG1 association was at SNP rs6731443 (P-value 2.5E-06), both intronic SNPs of the corresponding genes. We utilized the AceView tool [32] to gather current information for these genes. The ITGB5 gene maps to chromosome 3 at 3q21.2, and covers 139.46 kb, from 124620254 to 124480791 (NCBI 37, August 2010). ITGB5 is a highly expressed gene in many tissues, including lung, with 15 different mRNAs and 10 probable alternative promoters. The gene has been studied widely and its product, integrin β5, is involved in cell adhesion and integrin-mediated signaling. Of most relevance to AHR, one study by Tatler, et al., investigated whether contraction agonists could promote TGF-β activation in human airway smooth muscle (HASM) cells and found that integrin ανβ5 was a mediator of this activation [33]. Specifically, lysophosphatidic acid (LPA) and methacholine activation of TGF-β, a cytokine that has been implicated in airway remodeling in asthma, was shown to be abrogated by an ανβ5 blocking antibody. Further, the β5 cytoplasmic domain of integrin ανβ5 was found to be involved in LPA activation because a polymorphism in this subunit reversed the integrin ανβ5 activation of TGF-β. Comparison of normal and asthmatic HASM cells found that asthmatic HASM had increased LPA-induced, integrin ανβ5-mediated TGF-β activity, and that this increase was not due to increased cell surface expression of integrin ανβ5. While we were unable to replicate the primary GWAS ITGB5 associations in two independent populations, and there was no evidence that the corresponding SNPs were eQTL in whole blood of CAMP subjects, the study by Tatler, et al., suggests that polymorphisms in ITGB5 could play a role in modulating HASM response to contraction agonists via mechanisms that do not involve changes to the expression levels of the genes whose products form integrin ανβ5. Further, it is possible that our top ITGB5 SNPs are eQTLs of this gene HASM cells and not in whole blood, which is the tissue for which we currently have eQTL measures. Thus, it is possible that the ITGB5 association we observed truly reflects a genetic change that modulates AHR.
Our association results for AGFG1 variants were stronger than for ITGB5 ones in that we were able to replicate them in one of two independent populations (lowest combined P-value across all cohorts was rs6731443 2.14E-06) and the corresponding SNPs were nominally significant eQTL of AGFG1. The AGFG1 gene maps to chromosome 2, at 2q36.3, covering 89.07 kb, from 228336873 to 228425938 (NCBI 37, August 2010). AGFG1 is also a highly expressed gene in many tissues, including lung, with 17 different mRNAs and 4 probable alternative promoters. The protein encoded by this gene (a.k.a. RIP, HRB) binds the activation domain of the human immunodeficiency virus (HIV) Rev protein when Rev is assembled onto its RNA target, and is required for the nuclear export of Revdirected RNAs [34]. While many studies of AGFG1 focus on its relationship to HIV, AGFG1 has also been found to be involved in influenza A genome trafficking from the host nucleus to plasma membrane [35]. Because asthma development and severity are known to be influenced by respiratory pathogens, including influenza [36], current knowledge of AGFG1 suggests that if our association and eQTL data for this gene represent true biologically significant findings, their relationship with AHR may involve changes in immune response or susceptibility to external factors (e.g., influenza). Further functional validation of the eQTL results, including via quantitative RT-PCR of AGFG1 in subjects with various genotypes of the SNPs listed in Table 3 and for various tissues/cell types, would help clarify whether any relationship between the identified AGFG1 variants and its expression levels truly exists. Because the nominal replication of AGFG1 associations was observed in LHS, a clinical trial that measured baseline AHR in smokers who did not have current asthma, it is possible that the AGFG1 association reflects biological processes modulating AHR that are not unique to asthma.

Conclusions
An AHR GWAS among 994 asthma patients found that the most strongly associated SNPs, rs848788 (P-value 7.2E-07) and rs6731443 (P-value 2.5E-06), were in the ITGB5 and AGFG1 genes, respectively. The ITGB5 association did not replicate in two independent populations, nor was there any evidence that the corresponding SNP was an eQTL of ITGB5. The AGFG1 result replicated at a nominally significant level in one independent population of COPD subjects (LHS P-value 0.012), and the SNP had a nominally significant unadjusted P-value (0.0067) for being eQTL of AGFG1. While future functional studies are required to validate the potential involvement of these SNPs in modulating AHR, current knowledge of both genes suggests that our associations may represent biologically significant findings.