Search tips
Search criteria 


Logo of plosgenPLoS GeneticsSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)View this Article
PLoS Genet. 2012 July; 8(7): e1002824.
Published online 2012 July 5. doi:  10.1371/journal.pgen.1002824
PMCID: PMC3390407

Genome-Wide Association Analysis in Asthma Subjects Identifies SPATS2L as a Novel Bronchodilator Response Gene

Carole Ober, Editor


Bronchodilator response (BDR) is an important asthma phenotype that measures reversibility of airway obstruction by comparing lung function (i.e. FEV1) before and after the administration of a short-acting β2-agonist, the most common rescue medications used for the treatment of asthma. BDR also serves as a test of β2-agonist efficacy. BDR is a complex trait that is partly under genetic control. A genome-wide association study (GWAS) of BDR, quantified as percent change in baseline FEV1 after administration of a β2-agonist, was performed with 1,644 non-Hispanic white asthmatic subjects from six drug clinical trials: CAMP, LOCCS, LODO, a medication trial conducted by Sepracor, CARE, and ACRN. Data for 469,884 single-nucleotide polymorphisms (SNPs) were used to measure the association of SNPs with BDR using a linear regression model, while adjusting for age, sex, and height. Replication of primary P-values was attempted in 501 white subjects from SARP and 550 white subjects from DAG. Experimental evidence supporting the top gene was obtained via siRNA knockdown and Western blotting analyses. The lowest overall combined P-value was 9.7E-07 for SNP rs295137, near the SPATS2L gene. Among subjects in the primary analysis, those with rs295137 TT genotype had a median BDR of 16.0 (IQR = [6.2, 32.4]), while those with CC or TC genotypes had a median BDR of 10.9 (IQR = [5.0, 22.2]). SPATS2L mRNA knockdown resulted in increased β2-adrenergic receptor levels. Our results suggest that SPATS2L may be an important regulator of β2-adrenergic receptor down-regulation and that there is promise in gaining a better understanding of the biological mechanisms of differential response to β2-agonists through GWAS.

Author Summary

Bronchodilator response (BDR) is an important asthma phenotype that measures reversibility of airway obstruction by comparing lung function before and after the administration of short-acting β2-agonists, common medications used for asthma treatment. We performed a genome-wide association study of BDR with 1,644 white asthmatic subjects from six drug clinical trials and attempted to replicate these findings in 1,051 white subjects from two independent cohorts. The most significant associated variant was near the SPATS2L gene. We knocked down SPATS2L mRNA in human airway smooth muscle cells and found that β2-adrenergic receptor levels increased, suggesting that SPATS2L may be a regulator of BDR. Our results highlight the promise of pursuing GWAS results that do not necessarily reach genome-wide significance and are an example of how results from pharmacogenetic GWAS can be studied functionally.


Asthma is a chronic respiratory disease that affects over 20 million Americans and 300 million people worldwide [1], [2]. A hallmark characteristic of asthma is reversible airway obstruction, which is commonly measured via a bronchodilator response (BDR) test, in which the reduction of bronchoconstriction after administration of a short-acting reliever drug is quantified [3]. β2-agonists, the most common short-acting reliever drugs used during BDR tests and for asthma therapy, act in part by stimulating β2-adrenergic receptors (β2ARs) on airway smooth muscle cells to reduce bronchoconstriction via subsequent increases in cyclic adenosine monophosphate (cAMP) and protein kinase A (PKA) [3]. Although a comprehensive pathophysiologic understanding of BDR has not been obtained, it is a complex trait involving interactions among various tissues and cells, including inflammatory [4], airway epithelium [5], smooth muscle [6], and the autonomic nervous system [7]. In addition to being used for the diagnosis of asthma, BDR tests can be used to measure whether inhaled β2-agonists are effective in patients. Although short-acting β2-agonists are widely used clinically as asthma rescue medications, they are variably efficacious among patients [8]. Studying BDR may thus provide information regarding both the pathophysiology and pharmacogenetics of asthma.

The search for genetic variants that modify asthma susceptibility has resulted in the most recent multi-center asthma genome-wide association studies (GWAS) providing strong statistical evidence for the association of many genes, including the IKZF3-ZPBP2-GSDMB-ORMDL3 locus, HLA-DQ, IL1RL1, IL18RL1, IL33, TSLP, SLC22A5, SMAD3, and RORA, with asthma [9], [10]. Functional experiments to identify the role that these genes play in asthma pathophysiology are hindered by the complexity of the asthma phenotype. Familial aggregation [11] and genetic association studies [12] have provided suggestive evidence for a genetic contribution to interindividual differences in BDR. Candidate genes reported to be associated with BDR include β2-adrenergic receptor (ADRB2) [13], [14], adenylyl cyclase type 9 (ADCY9) [15], corticotrophin-releasing hormone receptor 2 (CRHR2) [16], and arginase 1 (ARG1) [17], [18]. While BDR is a complex phenotype, functional studies of BDR candidate genes are simpler than those for a general asthma phenotype because this pharmacogenetic phenotype can be readily simulated in vitro via stimulation of cells with β2-agonists.

In this study, we performed a GWAS of BDR in 1,644 non-Hispanic white asthmatics and found that the strongest evidence of association with BDR was at variants near the Homo sapiens spermatogenesis associated, serine-rich 2-like (SPATS2L) gene. We attempted to replicate the primary findings in two independent populations and investigated the function of SPATS2L via mRNA knockdown experiments and found evidence to support its involvement in BDR.


Primary BDR GWAS

Figure 1 is an overview of our study design. Characteristics of the subjects used in the primary GWAS are provided in Table 1. We utilized 1,644 non-Hispanic white subjects from six clinical trials to measure the association of SNPs to BDR. After QC filters, 469,884 SNPs genotyped in CAMP/LOCCS/LODO/Sepracor and either genotyped or imputed in CARE and ACRN were used to test for the association of SNPs to BDR. We utilized genotyped SNPs for CAMP/LOCCS/LODO/Sepracor because these cohorts, who were all genotyped using Illumina platforms, had the largest sample size. Due to the poor overlap of Illumina and Affymetrix platform SNPs, we utilized HapMap Phase 2 imputed SNPs for CARE and ACRN, so that the maximal number of SNPs in all cohorts could be analyzed. The quantile-quantile (QQ) and Manhattan plots revealed that the distribution of association P-values was similar to that expected for a null distribution and that no P-values met genome-wide statistically significant levels (Figures S1 and S2). To expand the primary association results further, all SNPs available in the June 2010 release of the 1000 Genome Project (1000GP) data were imputed using MaCH in each of the three primary groups of genotype data and overall BDR GWAS results were re-computed. Among SNPs contained in the primary GWAS, imputed and genotyped P-values were similar, particularly for those with low P-values (Figure S3). Some imputed regions had P-values lower than those of the primary GWAS, but the results in most of these regions were not supported by primary GWAS data (Figure S2). Thus, we proceeded to attempt to validate the primary GWAS findings based on the combined genotyped CAMP/LOCCS/LODO/Sepracor SNP results and HapMap Phase 2 imputed CARE and ACRN SNP results. The top five primary GWAS SNPs with P-value<1E-05 are in Table 2. Further details on these regions and all primary GWAS SNPs with P-value<1E-04 are in Tables S1 and S2. Further details on all 1000GP imputed GWAS SNPs with P-value<1E-05 are in Tables S3

Figure 1
Study overview.
Table 1
Characteristics of GWAS subjects.
Table 2
Primary GWAS Top Results (SNPs with Combined P-values <1e-05).

Replication of Primary GWAS Findings

We attempted to replicate in SARP all of the SNPs with primary GWAS P-values<1E-04 (Table S5). Three had nominally significant P-values (i.e. <0.05), and two of these associations supported the top 5 primary GWAS associations (Table 3). The lowest combined P-value for all primary GWAS plus SARP data was 7.7E-07 for rs295137. The region of BDR association spanning this SNP was in the 5′UTR region of SPATS2L, a gene of unknown in vivo function and paralog of SPATS2 (Figure 2). The effect of the rs295137 genotype on BDR is shown in Figure 3, and a plot of the residuals of the linear regression fit of BDR while adjusting for age, sex, and height is shown in Figure S4. We sought further evidence of association for the two SPATS2L SNPs with lowest P-values in our primary GWAS in a second independent population: DAG. There was no evidence for association in this cohort (rs295137 P-value = 0.21; rs295114 P-value = 0.21), and combined P-values for these two SNPs across all cohorts were 9.7E-07 and 1.6E-06 (Table 3). To investigate whether our top combined association represented a biologically significant finding, we sought experimental evidence that SPATS2L was involved in bronchodilator response.

Figure 2
Region of association near SPATS2L SNPs to BDR.
Figure 3
Summary of BDR by genotype of the SNP (rs295137) near SPATS2L with lowest P-value among all subjects in the primary populations.
Table 3
Replication study of top SNPs in two independent populations (SARP and DAG).

SPATS2L mRNA Expression Changes when the β2-Agonist Pathway Is Modified in Human Airway Smooth Muscle Cells

We found one public gene expression array experiment (GSE13168) that would help to address the question of whether SPATS2L is differentially expressed in response to changes in the BDR pathway. We compared the levels of expression of two SPATS2L and one SPATS2 probes in human airway smooth muscle (HASM) cells that stably expressed a PKA inhibitor vs. a GFP control at baseline and when stimulated with the pro-asthmatic cytokines interleukin-1β (IL1β), epidermal growth factor (EGF), or both. A trend of differential expression was observed for the SPATS2 and one SPATS2L probes, but not a second SPATS2L probe (Figure S5). None of the comparisons met a Benjamini-Hochberg adjusted significance threshold, but nominally significant P-values were obtained for the SPATS2 probe under all conditions and for one of the SPATS2L probes under the condition of EGF and IL1β stimulation (Table S6). According to the Gene Enrichment Profiler, the two SPATS2L probes are highly expressed in lung, and all three probes are highly expressed in smooth muscle, especially the SPATS2 one. Overall, there are strikingly different tissue-specific expression patterns for each probe (Figure S6).

Knockdown of SPATS2L mRNA Leads to Increased β2AR Levels

We further investigated the involvement of SPATS2L in the β2-adrenergic response pathway by knocking down SPATS2L mRNA using two different small interfering RNAs (siRNA) and measuring subsequent changes in β2AR protein levels. The knockdown efficiency of the siRNAs was >80% reduction of SPATS2L mRNA as measured by qRT-PCR, and the corresponding increases in β2AR (normalized against the control β-actin protein) levels were 1.88- (SD 0.41) and 1.86- (SD 0.30) fold for the two SPATS2L siRNAs (Figure 4).

Figure 4
Effect of siRNA-mediated SPATS2L knockdown on β2-adrenergic receptor levels in HASM cells.

Primary GWAS Results in Previously Identified BDR Candidate Genes

The association of SNPs with BDR at SNPs in/near genes (i.e. ADRB2, ADCY9, CRHR2, ARG1) previously reported as being associated with BDR was measured in our primary GWAS imputed data (Figure S7). Nominally significant (P-value<0.05) SNPs were found in ADCY9, CRHR2, and ARG1, but not in ADRB2 (Tables S7 and S8). The SNPs with lowest P-values within 50 kb of these genes were: rs2531988 for ADCY9 (3.2E-03), rs12533248 for CRHR2 (0.029), and rs6929820 for ARG1 (0.012).


In recent years, many GWAS that have successfully identified risk-modifying loci for a wide range of complex diseases have been published, but progress toward understanding how the loci and genes identified are functionally related to diseases has been slow [19]. The relationship of genes and gene variants to pharmacogenetic traits is often easier to test functionally than that for complex diseases because pharmacogenetic traits are more amenable to in vitro testing. However, compared to GWAS of complex diseases, GWAS of pharmacogenetic traits have been challenged by the relatively small size of drug clinical trials, which has caused many studies to be underpowered for obtaining genome-wide significant associations [20]. Nonetheless, successful pharmacogenetic GWAS have led to the identification of loci involved in modulating response to inhaled corticosteroids among asthma patients [21], warfarin dose [22], and lipid-lowering response to statins [23].

One of the difficulties specific to BDR GWAS is the complexity of the BDR phenotype. Regardless of how it is quantified, BDR is highly variable among asthma patients due to the time-dependent variation in baseline FEV1 and the influence of external environmental factors [24]. BDR can be quantified in various ways, with slightly varying resulting classification of patients as responders or non-responders. For our study, we selected the definition most widely utilized in clinical and human asthma research settings: percent change in baseline FEV1 following administration of a standard dose of short-acting inhaled β2-agonist [25]. We have attempted to control for the known relationship between baseline lung function and BDR [24], [26] by (1) selecting a definition of BDR that standardizes the change in FEV1 by dividing by baseline FEV1 and (2) by using age, sex, and height, which together account for a large portion of the variability in baseline lung function, as covariates in our statistical models. Because BDR tests are routinely performed during asthma clinical trials to use as inclusion criteria and to monitor outcomes among patients, we were able to utilize subjects from several diverse asthma clinical trials that were not specifically designed to study the pharmacogenetics of BDR. Most of these trials included a wash-out period that reduces modification of BDR due to concomitant medication administration, but LOCCS and some CARE and ACRN subjects were administered BDR tests at a time when they were not necessarily off of medications (Table 1). Subjects without a placebo washout, and especially those who were on ICS (e.g. LOCCS subjects), may be expected to have less BDR than those on placebo. The relationship of the magnitude of BDR to the various gene loci could therefore be blunted and show a less significant relationship than would be expected if all studies had incorporated a placebo washout.

In addition to variable washout periods, the cohorts had other significant differences in their design. Two trials consisted of children with asthma (i.e. CAMP, CARE), while the others consisted mostly of adults. The gender composition varied from 25% to 62% male. We attempted to control for age, gender and height, all of which are known to influence BDR, by including them as covariates in the association analysis. The mean and range of BDR also varied among trials. Of most significance, because the Sepracor trial used BDR greater than 15% as a criterion for inclusion, its subjects had markedly greater BDR than those of other trials. We attempted to control for this difference and any other trial specific differences among the cohorts that were pooled in the primary analysis by including trial as a covariate in the association model. There were additional differences among trials that were not taken into account. For example, SARP and Sepracor were composed of subjects with more severe asthma than those of other cohorts. Some ACRN, CARE, and SARP subjects were administered a different amount of albuterol during their BDR tests than those of CAMP, Sepracor, LOCCS, and LODO. DAG subjects were administered a different beta-agonist (i.e. Salbutamol) at a different concentration than that used with subjects of all other trials. DAG and SARP subjects were not participants of clinical trials, so there was greater heterogeneity of subjects within those cohorts. Despite the heterogeneity among trials, we utilized as many subjects as possible in an attempt to increase our statistical power to detect associations of SNPs with BDR. We reasoned that any associations detected despite the heterogeneity of the trial populations would be those most likely to generalize to all asthma patients. Another expected consequence of the trial heterogeneity is that our association results do not replicate in all cohorts. While having the largest number of subjects provides the greatest statistical power to detect statistically significant associations that are most generalizable across the clinical trials, we may be missing associations that are specific to the individual trials. For example, the subjects within clinical trials representing different ranges of asthma severity, age, and baseline characteristics may have genetic associations that are unique to subjects with their specific trial characteristics. The small sample size of each individual clinical trial makes detection of trial-specific associations more challenging. Despite the cohort heterogeneity, our meta-analysis identified a strong association that suggests a novel gene is involved in BDR.

Our top association was at SNP rs295137, with a combined P-value across all cohorts of 9.7E-07. This P-value does not meet conventional genome-wide significance thresholds (e.g. Bonferroni corrected minimally significant P-value would be 0.05/469,884 = 1.1E-07), but performing searches through public data sources and the fact that other pharmacogenetic GWAS have discovered biologically important results without genome-wide significant associations led us to pursue our top association further. The region of association surrounding rs295137 is in the 5′UTR of SPATS2L (Figure 2). This gene maps to chromosome 2 at 2q33.1, covering 176.78 kb from 201170592 to 201347368 (NCBI 37, August 2010). According to data gathered via the AceView [27] tool, SPATS2L is a complex locus that may have at least 30 spliced variants, its in vivo function is unknown, and it is a highly expressed gene in many tissues, with the greatest number of GenBank accessions belonging to lung. In gene-trap experiments in myoblasts, SPATS2L (a.k.a. SGNP) was found to be involved in ribosomal biogenesis and translational control in response to oxidative stress [28].

The availability of one public expression array experiment that utilized HASM cells expressing a PKA inhibitor (PKI) to modify the β2-adrenergic pathway allowed us to perform a preliminary search for evidence that SPATS2L may be involved in BDR. We found that a probe for SPATS2, the paralog of SPATS2L, was significantly differentially expressed in PKI vs. control cells at baseline and when stimulated with pro-inflammatory cytokines (EGF, IL1β, or both). One SPAST2L probe followed this trend but had a nominally significant P-value only under the condition of stimulation with both EGF and IL1β, while the other SPATS2L probe did not exhibit any changes. As illustrated in Figure S6, the tissue-specific expression patterns of the three probes varied widely. While all were expressed in smooth muscle, the SPATS2 probe's relative expression in this tissue was markedly greater than that of the SPATS2L probes. Taken together, the expression patterns are consistent with tissue and isoform dependent changes in SPATS2L gene products. While the public dataset SPATS2L results were inconclusive based on the differences among probes, they suggested that SPATS2L expression may change when PKA is inhibited in HASM cells.

Knockdown of SPATS2L in HASM cells resulted in significantly increased β2AR protein levels, suggesting that SPATS2L may affect BDR by directly modulating β2AR protein expression. In HASM, β2-agonists exert their effects exclusively via the β2AR [6]. The relaxation of HASM occurs after the binding of β2-agonists to β2ARs via increased levels of cAMP followed by PKA activation. PKA activation leads to changes in gene transcription via activation of cAMP response element binding protein (CREB). Because β2ARs are the gateway to the effects of β2-agonists in HASM cells, modulations, such as SPATS2L inactivation, that increase the levels of β2ARs in HASM cells may lead to both greater relaxation in response to β2-agonists in the short term and greater differences in gene transcription in the longer term. Further study is needed to elucidate the precise mechanism by which SPATS2L regulates β2AR and consequently modifies BDR. Among our primary GWAS subjects, those whose SPATS2L SNP rs295137 has the TT genotype have greater BDR than those with CT or TT genotypes (median BDR 16.0 vs. 10.9). In one of the simplest scenarios, it is possible that the increased BDR among subjects with the TT genotype results from this genotype playing a direct role in decreasing transcription of SPATS2L, which in turn results in increased β2AR levels. Further work is required to understand how specific SNP associations in/near SPATS2L affect SPATS2L function and/or expression and how such effects impact β2AR signaling and BDR. Because the observed influence of our most strongly associated SNP genotype on BDR is relatively small (Figure 3), our current data do not support the development of any personalized therapeutics based solely on variants in/near SPATS2L.

In addition to studying top primary GWAS SNPs, we attempted to replicate findings from previous BDR candidate gene association studies. Specifically, we measured association between BDR and ADRB2 [13], [14], ADCY9 [15], CRHR2 [16], and ARG1 [17] variants. Notably, these previous findings are not entirely independent of those from the current GWAS: CAMP was a primary population utilized to identify associations in ADRB2, ADCY9, CRHR2, and ARG1 in previous reports; LODO and Sepracor were replication populations in the CRHR2, and ARG1 reports; and LOCCS was a replication population in the ARG1 report. At a nominal significance level, we replicated gene-level associations for all of the candidate genes other than ADRB2. This gene, which encodes the β2AR, is the most studied gene related to BDR and SNPs and haplotypes in this gene have been related to decreased pulmonary function [29], response to β2-agonist treatment [30], an increased frequency of asthma exacerbations [31], and BDR [13], [14]. Initial reports of ADRB2 associations were very promising and suggested that variants of this candidate gene would be reliable markers of BDR pharmacogenetics. However, a meta-analysis of 21 studies of ADRB2 polymorphisms found that most of the positive associations identified in individual studies cannot be compared to findings in other studies due to different study designs, phenotypes considered and selective reporting, making the number of variants with true replications very small and questioning the utility of ADRB2 polymorphisms for generalizable pharmacogenetic tests [32]. Our inability to find associations with ADRB2 variants is consistent with the view that BDR genetics are complex: no individual SNPs or genes are responsible for a large proportion of BDR variability observed among all asthmatics. Our results suggest that genes other than the previously identified candidate genes are more strongly associated with BDR and that functional studies of these regions may yield important insights into BDR biology despite not having strong effects or generalizing to all populations.

In summary, a BDR GWAS among asthma patients from eight cohorts found that the most strongly associated SNP, rs295137, had a combined P-value of 9.7E-07. This association led us to SPATS2L, a gene of unknown in vivo function that we showed may be involved in BDR via the down-regulation of β2AR levels. Our results support the notion that there is promise in pursuing GWAS results that do not necessarily reach genome-wide significance and are an example of the way in which results from pharmacogenetic GWAS can be studied functionally.

Materials and Methods

Ethics Statement

Each study was approved by the Institutional Review Board of the corresponding institution, which ensured that all procedures followed were in accordance with the ethical standards of the responsible committee on human experimentation. Informed consent was obtained for all study participants.


The primary group of subjects consisted of 1,644 non-Hispanic white asthmatics from the following drug clinical trials: Childhood Asthma Management Program (CAMP) [33], Leukotriene Modifier or Corticosteroid Salmeterol study (LOCCS) [34], Effectiveness of Low Dose Theophylline as an Add-on Treatment in Asthma trial (LODO) [35], a medication trial conducted by Sepracor, Inc. [36], [37], and subsets of clinical trials within the Childhood Asthma Research and Education (CARE) network [38], and the Asthma Clinical Research Network (ACRN) [39] participating in the NHLBI SNP Health Association Resource (SHARe) Asthma Resource project (SHARP). Some basic characteristics of these cohorts are in Table 1 and further details are provided in Text S1. BDR tests were performed according to American Thoracic Society criteria with Albuterol as the β2-agonist [25], unless otherwise noted. Baseline BDR measures were utilized, and BDR was quantified as the percent change in FEV1 in response to administration of a β2-agonist [i.e. (post-BD FEV1 – pre-BD FEV1)/pre-BD FEV1].

Genotyping and Quality Control

Genome-wide genotyping for CAMP subjects (n = 546) was performed on the HumanHap550 Genotyping BeadChip or Infinium HD Human610-Quad BeadChip by Illumina, Inc (San Diego, CA) at the Channing Laboratory. LOCCS (n = 135), LODO (n = 114), and Sepracor (n = 401) subjects were genotyped at the Riken Center for Genomic Medicine using the Infinium HD Human610-Quad BeadChip. CARE (n = 207) and ACRN (n = 241) subjects were genotyped on Affymetrix 6.0 genotyping chip by Affymetrix, Inc. (Santa Clara, CA). Data from those subjects genotyped using Illumina technologies was combined into a primary dataset with 469,884 overlapping SNPs having missingness <1%, passing HWE (P-value threshold of 1E-03), and having minor allele frequency (MAF)>0.05. EIGENSTRAT was used to identify 23 outliers (not included in counts above) based on being outside of six standard deviations of the top four principal components during five iterations [40]. The genomic inflation factor (λGC) of the remaining 1,196 subjects was 1.002, demonstrating minimal population stratification. CARE and ACRN dataset quality control (QC) also included the removal of four related subjects (i.e. CARE siblings), SNPs with MAF<0.05, missingness >5%, or not passing HWE based on a threshold of 1E-03. The λGC for CARE and ACRN genotype data were 1.02 and 0.98, demonstrating minimal population stratification among subjects within each group. Comprehensive genotyping and QC measures are provided in Text S1.

Statistical Analysis

Due to the poor overlap among SNPs genotyped on the Illumina and Affymetrix platforms, imputation of all SNPs available in HapMap Phase 2 Release 22 CEU data using the Markov Chain Haplotyping software (MaCH) [41] was performed for ACRN and CARE genotyped data. The primary GWAS consisted in the set of 469,884 SNPs that were successfully genotyped in those cohorts using Illumina arrays (i.e., CAMP/Sepracor/LOCCS/LODO) and imputed with HapMap Phase 2 data in those cohorts genotyped with Affymetrix arrays (i.e., ACRN and CARE) with a ratio of empirically observed dosage variance to the expected (binomial) dosage variance greater than 0.3, indicating good quality of imputation.

To expand the association results, imputation of all SNPs available in the June 2010 release of the 1000 Genome Project (1000GP) data using MaCH was performed for each of the three primary groups of genotype data. An overlapping set of 4,571,615 imputed SNPs had a MAF>0.05 and ratio of empirically observed dosage variance to the expected (binomial) dosage variance greater than 0.5, indicating good quality of imputation.

The association of SNPs with BDR was measured with a linear regression model as implemented in PLINK [42] in the three sets of data: 1) CAMP/Sepracor/LOCCS/LODO, 2) ACRN, 3) CARE. Association of imputed SNPs was carried out using dosage data. Covariates for the CAMP/Sepracor/LOCCS/LODO group included age, gender, height, and study. Covariates for the CARE and ACRN groups included age, gender, and height. To get the primary GWAS results, CARE and ACRN P-values were combined with those of the CAMP/Sepracor/LOCCS/LODO group by using Liptak's combined probability method [43] where hypothesis tests in CARE and ACRN had one-sided alternatives, based on the direction of the association in CAMP/Sepracor/LOCCS/LODO, so that SNPs with association tests in opposite directions would not produce inappropriately small P-values. The overall λGC was 1.002 in the primary set of GWAS results and 1.000 in the 1000GP imputed data GWAS.

Plots of association results near specific genes were created using LocusZoom with the hg18/1000 Genomes June 2010 CEU GenomeBuild/LD Population [44].

Replication Studies

1) Severe Asthma Research Program (SARP)

Replication of primary GWAS P-values was attempted in 501 European American subjects from SARP who were recruited to meet the American Thoracic Society (ATS) definition of mild to severe asthmatics, with enrichment for severe asthma [45]. Genotyping was performed on the Illumina Infinium HD Human 610-Quad Bead Chip at Wake Forest University. Tests of BDR association were performed using linear regression in R [46] with age, sex, and height as covariates.

2) Dutch Asthma GWAS (DAG)

Replication of two SNPs near SPATS2L with lowest P-values in the primary GWAS was attempted in 550 DAG subjects with asthma and BDR data that were phenotyped in a single center (i.e, Beatrixoord, Haren, the Netherlands) [18], [47]. All patients refrained from bronchodilator use for at least 8 hours, and stopped inhaled steroids for 2 weeks before pulmonary function testing if clinically possible. Reversibility was assessed measuring spirometry before and 15 minutes after inhaling 800 µg of salbutamol. Genotyping was performed using the Hapmap 317K platform or Illumina 370 Duo Chip. Tests of BDR association were performed using linear regression, with age, gender and height as covariates using Predictive Analysis Software (PASW) version 18.0.

Genome-Wide Gene Expression Data

The publicly available Gene Expression Omnibus (GEO) dataset, GSE13168, corresponding to an experiment in which human airway smooth muscle (HASM) cell cultures were generated from four donor trachea to test for the effects of glucocorticoids and PKA inhibition on the HASM transcriptome using the Affymetrix Human Genome U133A platform was used [48]. We tested for the involvement of our top primary GWAS gene in the β2-adrenergic pathway by comparing the differential expression of genes in cells stably expressing a PKA inhibitor (PKI) vs. control at baseline and in the presence of pro-inflammatory cytokines interleukin-1β (IL1β), epidermal growth factor (EGF), or both. The expression array contained two SPATS2L probes (i.e., 215617_at, 222154_s_at) and one SPATS2 probe (i.e., 218324_s_at). The probe for the paralog of SPATS2L was included to account for the possibility of non-specific binding of SPATS2L mRNA to the SPATS2 probe. Analyses were conducted in R [46]. Pre-processing of raw signal intensities was performed with RMA [49] and differential expression was quantified using the limma package [50]. Tissue-specific expression of these probes was assessed using 557 microarrays from 126 human normal primary tissues in the Gene Enrichment Profiler [51].

SPATS2L siRNA Knockdown and β2AR Western Blotting Analysis

Primary HASM cells were isolated from aborted lung transplant donors with no chronic illness. The tissue was obtained from the National Disease Resource Interchange (NDRI) and their use approved by the University of Pennsylvania IRB. HASM cell cultivation and characterization were described previously [52], [53]. Passages 4 to 7 HASM cells maintained in Ham's F12 medium supplemented with 10% FBS were used in all experiments. 2×105 HASM cells were grown overnight and then transfected with 50 nM siRNA by using DharmaFECT 1 reagent (Thermo Scientific, Lafayette, CO, USA). About 72 h post transfection, cells were washed with PBS and lysed with NP-40 lysing buffer (50 mM Tris-HCl pH7.5, 150 mM NaCl, 0.5% Nonidet P-40) containing protease inhibitor cocktail (Roche Diagnostics Corporation, Indianapolis, IN, USA). Protein samples were denatured 10 min at 50°C, separated on NuPAGE 4–12% Bis-Tris gels (Invitrogen, Grand Island, NY, USA) and transferred to PVDF membranes (Bio-Rad Laboratories, Hercules, CA, USA). Immunoblot signals were developed using SuperSignal West Pico (Pierce Protein Research Products, Thermo Fisher Scientific, Rockford, IL, USA) and quantified by ImageJ software. Non-targeting control siRNA, SPATS2L-specific siRNA 1 (sense 5′- guc agu cca uug auu guc u(dT)(dT)-3′, antisense 5′- aga caa uca aug gac uga c(dT)(dT) -3′) and SPATS2L-specific siRNA 2 (sense 5′-caa ccu gug uug uag cag u(dT)(dT)-3′, antisense 5′- acu gcu aca aca cag guu g(dT)(dT) -3′) were obtained from Sigma-Aldrich (Mission siRNA; St. Louis, MO, USA). Antibodies for β2AR (H20) and β-actin were from Santa Cruz Biotechnology, Inc. (Santa Cruz, CA, USA). Experiments were performed in triplicate.

Supporting Information

Figure S1

Quantile–quantile plot. Comparison of primary GWAS P-values to those expected for a null distribution. There is little evidence of deviation of measures at the tail, obscuring the distinction among SNPs having low p-values representing true associations vs. those SNPs having low p-values by chance.


Figure S2

Manhattan plot. The x-axis denotes position along each chromosome. The y-axis denotes −Log10(P) corresponding to association P-values. The 1000GP imputed results are in light blue and green. The primary GWAS results are in dark blue and green. Some of the lowest imputed P-values did not have corresponding primary GWAS p-values, while other regions with low P-values contained both primary and imputed GWAS results. We prioritized regions based on primary GWAS results.


Figure S3

P-value comparison between Primary GWAS vs. 1000GP imputed results, for those SNPs contained in both datasets. The results are highly correlated (r2 = 0.99).


Figure S4

Residuals of linear regression fit of BDR ~ age + sex + height vs. rs295137 genotypes.


Figure S5

Distribution of normalized signal intensities across different experimental conditions for two SPATS2L and one SPATS2 probes corresponding to a subset of GSE13168 arrays where the effects of protein kinase A inhibition (PKI) in human airway smooth muscle cells was assessed at baseline and following stimulation with epidermal growth factor (EGF), interleukin 1 beta (IL1b), or both. Corresponding P-values and log-Fold Change differences are in Table S4. A) SPATS2 probe 218324_s_at, B) SPATS2L probe 215617_at, C) SPATS2L probe 222154_s_at.


Figure S6

Tissue-specific gene expression patterns of Affy U133A probes for SPATS2 and SPATS2L (a.k.a. LOC26010).


Figure S7

Plots of imputed association results near previously identified BDR candidate genes. A) ADCY9, B) ADRB2 C) ARG1, D) CRHR2.


Table S1

Primary GWAS SNP details for SNPs with Primary GWAS P<1E-04.


Table S2

Primary GWAS P-values and Beta coefficients for SNPs with Primary GWAS P<1E-04. CARE and ACRN P-values are 1-sided based on the direction in CAMP/LOCCS/LODO/Sepracor. CAMP/LOCCS/LODO/Sepracor reference allele was minor allele as shown in Table 1. SE = Standard Error for corresponding Beta coefficient. Rsq = MACH R-squared value for imputed SNP. Combined P-values were obtained using the Liptak method, with weights proportional to population size.


Table S3

Primary GWAS 1000GP Imputed SNP details for SNPs with P<1E-05. MAF = Minor Allele Frequency. Rsq = MACH R-squared value for imputed SNP.


Table S4

Primary GWAS 1000GP Imputed SNP P-values and Beta coefficients for SNPs with P<1E-05. CARE and ACRN P-values are 1-sided based on the direction in CAMP/LOCCS/LODO/Sepracor. SE = Standard Error for corresponding Beta coefficient. Combined P-values were obtained using the Liptak method, with weights proportional to population size.


Table S5

SARP Replication Results for SNPs with Primary GWAS P<1E-04. P-values are 1-sided based on the direction in CAMP/LOCCS/LODO/Sepracor. MAF = minor allele frequency.


Table S6

Unadjusted P-values and log fold-change quantifying differential expression of SPATS2L and SPATS2 probes for GEO dataset GSE13168, in which human ASM cell lines expressing a PKA inhibitor vs. a GFP control were compared at baseline and when stimulated with IL1b, EGF, or both.


Table S7

Primary GWAS 1000GP Imputed SNP details for SNPs that have nominally significant p-values (<0.05) in, or within 50,000 KB of, genes (i.e. ADRB2, ADCY9, CRHR2, ARG1) previously identified as being associated with BDR. There were no such SNPs near ADRB2. MAF = Minor Allele Frequency. Rsq = MACH R-squared value for imputed SNP.


Table S8

Primary GWAS 1000GP Imputed SNP P-values and Beta coefficients for SNPs that have nominally significant p-values (<0.05), and are in, or within 50,000 KB of, genes (i.e. ADRB2, ADCY9, CRHR2, ARG1) previously identified as being associated with BDR. There were no such SNPs near ADRB2. CARE and ACRN P-values are 1-sided based on the direction in CAMP/LOCCS/LODO/Sepracor. SE = Standard Error for corresponding Beta coefficient. Combined P-values were obtained using the Liptak method, with weights proportional to population size.


Text S1

Detailed subject characteristics and genotype data quality control.



We thank all CAMP subjects for their ongoing participation in this study. We acknowledge the CAMP investigators and research team, supported by the National Heart, Lung, and Blood Institute (NHLBI), for collection of CAMP Genetic Ancillary Study data. All work on data collected from the CAMP Genetic Ancillary Study was conducted at the Channing Laboratory of the Brigham and Women's Hospital under appropriate CAMP policies and human subject's protections. We thank the American Lung Association Asthma Clinical Research Centers (ALA-ACRC) for use of LOCCS and LODO study samples.


The authors have declared that no competing interests exist.

The CAMP Genetics Ancillary Study is supported by U01 HL075419, U01 HL65899, P01 HL083069, R01 HL086601, and T32 HL07427 from the NHLBI, National Institutes of Health (NIH). Additional support was provided by NIH U10 HL064287, U10 HL064288, U10 HL064295, U10 HL064305, U10 HL064307, U01 HL064313, RC2 HL101487, and U01 HL65899, an NIH Pharmacogenomics Research Network (PGRN)–RIKEN Center for Genomic Medicine (CGM) Global Alliance and by funding from the BioBank Japan project that was supported by the Ministry of Education, Culture, Sports, Sciences, and Technology of the Japanese government. BEH was supported by NIH K99 HL105663. Dutch Asthma GWAS (DAG) was supported by the Netherlands Asthma Foundation grant AF (AF 95.09, AF 98.48, AF, and AF and a grant from the University Medical Center Groningen. GHK was supported by a Ter Meulen Fund grant from the Royal Netherlands Academy of Arts and Sciences. RAP was supported by NIH HL097796 and ES013508. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


1. American Lung Association. Trends in Asthma Morbidity and Mortality. New York, NY: Epidemiology and Statistics Unit, Research and Program Services, American Lung Association; 2006.
2. Global Initiative for Asthma Management and Prevention. NHLBI/WHO Workshop Report. National Institutes of Health, Bethesda: US Department of Health and Human Services; 1995.
3. Nelson HS. Beta-adrenergic bronchodilators. N Engl J Med. 1995;333:499–506. [PubMed]
4. Loza MJ, Penn RB. Regulation of T cells in airway disease by beta-agonist. Front Biosci (Schol Ed) 2010;2:969–979. [PMC free article] [PubMed]
5. Salathe M. Effects of beta-agonists on airway epithelial cells. J Allergy Clin Immunol. 2002;110:S275–281. [PubMed]
6. Shore SA, Moore PE. Regulation of beta-adrenergic responses in airway smooth muscle. Respir Physiol Neurobiol. 2003;137:179–195. [PubMed]
7. Jartti T. Asthma, asthma medication and autonomic nervous system dysfunction. Clin Physiol. 2001;21:260–269. [PubMed]
8. Drazen JM, Silverman EK, Lee TH. Heterogeneity of therapeutic responses in asthma. Br Med Bull. 2000;56:1054–1070. [PubMed]
9. Moffatt MF, Gut IG, Demenais F, Strachan DP, Bouzigon E, et al. A large-scale, consortium-based genomewide association study of asthma. N Engl J Med. 2010;363:1211–1221. [PMC free article] [PubMed]
10. Torgerson DG, Ampleford EJ, Chiu GY, Gauderman WJ, Gignoux CR, et al. Meta-analysis of genome-wide association studies of asthma in ethnically diverse North American populations. Nat Genet. 2011;43:887–892. [PMC free article] [PubMed]
11. Niu T, Rogus JJ, Chen C, Wang B, Yang J, et al. Familial aggregation of bronchodilator response: a community-based study. Am J Respir Crit Care Med. 2000;162:1833–1837. [PubMed]
12. Hawkins GA, Weiss ST, Bleecker ER. Clinical consequences of ADRbeta2 polymorphisms. Pharmacogenomics. 2008;9:349–358. [PubMed]
13. Drysdale CM, McGraw DW, Stack CB, Stephens JC, Judson RS, et al. Complex promoter and coding region beta 2-adrenergic receptor haplotypes alter receptor expression and predict in vivo responsiveness. Proc Natl Acad Sci U S A. 2000;97:10483–10488. [PubMed]
14. Silverman EK, Kwiatkowski DJ, Sylvia JS, Lazarus R, Drazen JM, et al. Family-based association analysis of beta2-adrenergic receptor polymorphisms in the childhood asthma management program. J Allergy Clin Immunol. 2003;112:870–876. [PubMed]
15. Tantisira KG, Small KM, Litonjua AA, Weiss ST, Liggett SB. Molecular properties and pharmacogenetics of a polymorphism of adenylyl cyclase type 9 in asthma: interaction between beta-agonist and corticosteroid pathways. Hum Mol Genet. 2005;14:1671–1677. [PubMed]
16. Poon AH, Tantisira KG, Litonjua AA, Lazarus R, Xu J, et al. Association of corticotropin-releasing hormone receptor-2 genetic variants with acute bronchodilator response in asthma. Pharmacogenet Genomics. 2008;18:373–382. [PMC free article] [PubMed]
17. Litonjua AA, Lasky-Su J, Schneiter K, Tantisira KG, Lazarus R, et al. ARG1 is a novel bronchodilator response gene: screening and replication in four asthma cohorts. Am J Respir Crit Care Med. 2008;178:688–694. [PMC free article] [PubMed]
18. Vonk JM, Postma DS, Maarsingh H, Bruinenberg M, Koppelman GH, et al. Arginase 1 and arginase 2 variations associate with asthma, asthma severity and beta2 agonist and steroid response. Pharmacogenet Genomics. 2010;20:179–186. [PubMed]
19. Manolio TA. Genomewide association studies and assessment of the risk of disease. N Engl J Med. 2010;363:166–176. [PubMed]
20. Daly AK. Genome-wide association studies in pharmacogenomics. Nat Rev Genet. 2010;11:241–246. [PubMed]
21. Tantisira KG, Lasky-Su J, Harada M, Murphy A, Litonjua AA, et al. Genomewide association between GLCCI1 and response to glucocorticoid therapy in asthma. N Engl J Med. 2011;365:1173–1183. [PMC free article] [PubMed]
22. Takeuchi F, McGinnis R, Bourgeois S, Barnes C, Eriksson N, et al. A genome-wide association study confirms VKORC1, CYP2C9, and CYP4F2 as principal genetic determinants of warfarin dose. PLoS Genet. 2009;5:e1000433. doi: 10.1371/journal.pgen.1000433. [PMC free article] [PubMed]
23. Barber MJ, Mangravite LM, Hyde CL, Chasman DI, Smith JD, et al. Genome-wide association of lipid-lowering response to statins in combined study populations. PLoS ONE. 2010;5:e9763. doi: 10.1371/journal.pone.0009763. [PMC free article] [PubMed]
24. Sharma S, Litonjua AA, Tantisira KG, Fuhlbrigge AL, Szefler SJ, et al. Clinical predictors and outcomes of consistent bronchodilator response in the childhood asthma management program. J Allergy Clin Immunol. 2008;122:921–928 e924. [PMC free article] [PubMed]
25. American Thoracic Society. Lung function testing: selection of reference values and interpretative strategies. Am Rev Respir Dis. 1991;144:1202–1218. [PubMed]
26. Waalkens HJ, Merkus PJ, van Essen-Zandvliet EE, Brand PL, Gerritsen J, et al. Assessment of bronchodilator response in children with asthma. Dutch CNSLD Study Group. Eur Respir J. 1993;6:645–651. [PubMed]
27. Thierry-Mieg D, Thierry-Mieg J. AceView: a comprehensive cDNA-supported gene and transcripts annotation. Genome Biol. 2006;7(Suppl 1):S12 11–14. [PMC free article] [PubMed]
28. Zhu CH, Kim J, Shay JW, Wright WE. SGNP: an essential Stress Granule/Nucleolar Protein potentially involved in 5.8s rRNA processing/transport. PLoS ONE. 2008;3:e3716. doi: 10.1371/journal.pone.0003716. [PMC free article] [PubMed]
29. Israel E, Drazen JM, Liggett SB, Boushey HA, Cherniack RM, et al. The effect of polymorphisms of the beta(2)-adrenergic receptor on the response to regular use of albuterol in asthma. Am J Respir Crit Care Med. 2000;162:75–80. [PubMed]
30. Israel E, Chinchilli VM, Ford JG, Boushey HA, Cherniack R, et al. Use of regularly scheduled albuterol treatment in asthma: genotype-stratified, randomised, placebo-controlled cross-over trial. Lancet. 2004;364:1505–1512. [PubMed]
31. Taylor DR, Drazen JM, Herbison GP, Yandava CN, Hancox RJ, et al. Asthma exacerbations during long term beta agonist use: influence of beta(2) adrenoceptor polymorphism. Thorax. 2000;55:762–767. [PMC free article] [PubMed]
32. Contopoulos-Ioannidis DG, Alexiou GA, Gouvias TC, Ioannidis JP. An empirical evaluation of multifarious outcomes in pharmacogenetics: beta-2 adrenoceptor gene polymorphisms in asthma treatment. Pharmacogenet Genomics. 2006;16:705–711. [PubMed]
33. Childhood Asthma Management Program Research Group. The Childhood Asthma Management Program (CAMP): design, rationale, and methods. Control Clin Trials. 1999;20:91–120. [PubMed]
34. Peters SP, Anthonisen N, Castro M, Holbrook JT, Irvin CG, et al. Randomized comparison of strategies for reducing treatment in mild persistent asthma. N Engl J Med. 2007;356:2027–2039. [PubMed]
35. American Lung Association Asthma Clinical Research Centers. Clinical trial of low-dose theophylline and montelukast in patients with poorly controlled asthma. Am J Respir Crit Care Med. 2007;175:235–242. [PubMed]
36. Silverman ES, Palmer LJ, Subramaniam V, Hallock A, Mathew S, et al. Transforming growth factor-beta1 promoter polymorphism C-509T is associated with asthma. Am J Respir Crit Care Med. 2004;169:214–219. [PubMed]
37. Silverman ES, Baron RM, Palmer LJ, Le L, Hallock A, et al. Constitutive and cytokine-induced expression of the ETS transcription factor ESE-3 in the lung. Am J Respir Cell Mol Biol. 2002;27:697–704. [PubMed]
38. Guilbert TW, Morgan WJ, Krawiec M, Lemanske RF, Jr, Sorkness C, et al. The Prevention of Early Asthma in Kids study: design, rationale and methods for the Childhood Asthma Research and Education network. Control Clin Trials. 2004;25:286–310. [PubMed]
39. Denlinger LC, Sorkness CA, Chinchilli VM, Lemanske RF., Jr Guideline-defining asthma clinical trials of the National Heart, Lung, and Blood Institute's Asthma Clinical Research Network and Childhood Asthma Research and Education Network. J Allergy Clin Immunol. 2007;119:3–11; quiz 12–13. [PMC free article] [PubMed]
40. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38:904–909. [PubMed]
41. Willer CJ, Sanna S, Jackson AU, Scuteri A, Bonnycastle LL, et al. Newly identified loci that influence lipid concentrations and risk of coronary artery disease. Nat Genet. 2008;40:161–169. [PubMed]
42. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–575. [PubMed]
43. Liptak T. On the combination of independent tests. Magyar Tud Akad Mat Kutato Int Kozl. 1958;3:171–197.
44. Pruim RJ, Welch RP, Sanna S, Teslovich TM, Chines PS, et al. LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics. 2010;26:2336–2337. [PMC free article] [PubMed]
45. Moore WC, Bleecker ER, Curran-Everett D, Erzurum SC, Ameredes BT, et al. Characterization of the severe asthma phenotype by the National Heart, Lung, and Blood Institute's Severe Asthma Research Program. J Allergy Clin Immunol. 2007;119:405–413. [PMC free article] [PubMed]
46. R Development Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2008.
47. Koppelman GH, Meyers DA, Howard TD, Zheng SL, Hawkins GA, et al. Identification of PCDH1 as a novel susceptibility gene for bronchial hyperresponsiveness. Am J Respir Crit Care Med. 2009;180:929–935. [PMC free article] [PubMed]
48. Misior AM, Deshpande DA, Loza MJ, Pascual RM, Hipp JD, et al. Glucocorticoid- and protein kinase A-dependent transcriptome regulation in airway smooth muscle. Am J Respir Cell Mol Biol. 2009;41:24–39. [PMC free article] [PubMed]
49. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy–analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20:307–315. [PubMed]
50. Smyth GK. Limma: linear models for microarray data. In: Gentleman R, Carey V, Dudoit S, Irizarry R, Huber W, editors. Bioinformatics and Computational Biology Solutions using R and Bioconductor. New York: Springer; 2005. pp. 397–420.
51. Benita Y, Cao Z, Giallourakis C, Li C, Gardet A, et al. Gene enrichment profiles reveal T-cell development, differentiation, and lineage-specific transcription factors including ZBTB25 as a novel NF-AT repressor. Blood. 2010;115:5376–5384. [PubMed]
52. Panettieri RA, Murray RK, DePalo LR, Yadvish PA, Kotlikoff MI. A human airway smooth muscle cell line that retains physiological responsiveness. Am J Physiol. 1989;256:C329–335. [PubMed]
53. Cooper PR, Mesaros AC, Zhang J, Christmas P, Stark CM, et al. 20-HETE mediates ozone-induced, neutrophil-independent airway hyper-responsiveness in mice. PLoS ONE. 2010;5:e10235. doi: 10.1371/journal.pone.0010235. [PMC free article] [PubMed]

Articles from PLoS Genetics are provided here courtesy of Public Library of Science