Nearly one-half of asthmatic patients do not respond to the most commonly prescribed controller therapy, inhaled corticosteroids (ICS). We conducted an expression quantitative trait loci (eQTL) analysis using >300 expression microarrays (from 117 lymphoblastoid cell lines) in corticosteroid (dexamethasone) treated and untreated cells derived from asthmatic subjects in the Childhood Asthma Management Program (CAMP) clinical trial. We then tested the associations of eQTL with longitudinal change in airway responsiveness to methacholine (LnPC20) on ICS. We identified 2484 cis-eQTL affecting 767 genes following dexamethasone treatment. A significant over-representation of lnPC20-associated cis-eQTL [190 single-nucleotide polymorphisms (SNPs)] among differentially expressed genes (odds ratio = 1.76, 95% confidence interval: 1.35–2.29) was noted in CAMP Caucasians. Forty-six of these 190 clinical associations were replicated in CAMP African Americans, including seven SNPs near six genes meeting criteria for genome-wide significance (P < 2 × 10−7). Notably, the majority of genome-wide findings would not have been uncovered via analysis of untreated samples. These results indicate that identifying eQTL after relevant environmental perturbation enables identification of true pharmacogenetic variants.

INTRODUCTION

Asthma is the most common chronic lung disease of childhood, affecting >3 million US children and leading to >200 000 hospitalizations per year (1). For long-term control of asthma, inhaled corticosteroids (ICS) are the mainstay of asthma therapy (2). However, clinical response to ICS therapy is widely variable among individuals (3) and is associated with important side effects.

Candidate genes studies have implicated multiple pharmacogenetic variants with replicated and/or functional association with ICS response in asthma (46). Genome-wide association studies (GWAS) have recently been applied to ICS response in asthma, identifying replicable associations in the T (7) and GLCCI1 (8) genes. One method to further enrich GWAS results is to study expression quantitative trait loci (eQTL) as an intermediate phenotype to identify functional candidates. In asthma, an eQTL approach was used to identify the highly replicated ORMDL3-GSDMB locus (9). For gene-environment studies, such as pharmacogenetics, several groups have studied eQTL in immortalized cells after exposure to environmental perturbations including radiation (10) and therapeutics such as corticocorticoids (11,12). These studies have consistently reported that cis-acting eQTL in association with environmentally mediated expression change are rare. Nonetheless, we previously reported ‘drug-specific’ eQTL associated with both TNC expression change following cellular corticosteroid stimulation and clinical response to ICS in asthma (12).

Recently, we identified a functional variant in GLCCI1 associated with reduced ICS response in four asthmatic cohorts that explains ∼6% of variation in ICS lung function response (8). Notably, this promoter variant encoded for a cis-eQTL for GLCCI1 expression that was not drug-specific—i.e. the cis-signal in cellular models was noted in both the baseline and the corticosteroid-treated samples. However, GLCCI1 was strongly differentially expressed in these samples. Moreover, dexamethasone-stimulated GLCCI1 expression in CAMP lymphoblastoid cell lines (LCLs) was associated with a robust response to ICS [odds ratio, 3.22; 95% confidence interval (CI), 1.41–7.38], a finding consistent with altered expression as the functional basis for the observed genotypic associations. Thus, assuming differential expression marks the subset of genes truly responsive to treatment, post-treatment cis-eQTL within differentially expressed genes are likely to affect clinical treatment response, even if they do not fully mediate the differential expression response. We therefore hypothesized that differential expression in conjunction with cis-eQTL in corticosteroid-treated cells derived from asthmatics would harbor important clinical ICS response variants.

RESULTS

We conducted an expression eQTL analysis using 117 LCLs in paired fashion [corticosteroid (dexamethasone) treated and untreated cells], which were derived from 117 asthmatic subjects in the Childhood Asthma Management Program (CAMP) clinical trial. We then tested the associations of eQTL with longitudinal change in airway responsiveness to methacholine (LnPC20) on ICS.

Among 21 175 gene probes (tagging 17 370 unique genes), we detected 5846 gene probes (tagging 5437 unique genes) differentially expressed between dexamethasone-treated LCLs and sham-treated LCLs at false discovery rate (FDR) <5%. We performed a pathway analysis using DAVID to further evaluate the differentially expressed genes (13,14) and identified 54 biologic processes that are significantly over-represented (with FDR-adjusted P value <0.01) in response to dexamethasone, including cell apoptosis, DNA replication and repair, immune response and leukocyte and lymphocyte activation (Supplementary Material, Table S3).

Focusing on dexamethasone-treated LCLs, we identified 2484 significant cis-eQTL [2376 single-nucleotide polymorphisms (SNPs) affecting 767 genes] (Supplementary Material, Table S4). Transcript variance explained at these loci was large (median = 0.16; range: 0.07–0.86). Genes harboring cis-eQTL identified in dexamethasone-treated cell lines were more likely to be differentially expressed by dexamethasone (OR = 2.7, 95% CI: 2.3–3.1).

To further check the robustness of our cis-eQTL results, we did two additional cis-eQTL analyses based on gene expression from dex-treated cell lines by further adjusting for a RNA batch variable and for the top 10 principal components from gene expression data, respectively. These two further adjustments did not significantly change our cis-eQTL results.

We tested these 2484 cis-eQTL for association with LnPC20 in 172 CAMP Caucasians randomized to ICS therapy. As shown in Table 1, 279 (12%) cis-eQTL are significantly associated with LnPC20. There is a greater degree of enrichment in the dex-responsive genes (i.e. genes differentially expressed between dexamethasone-treated cell lines and sham-treated cell lines) in this group of 279 cis-eQTL that were associated with LnPC20, with 190 associated versus 89 not associated with a differentially expressed gene (14 versus 9%) with an OR of 1.76 (1.35–2.30). To validate the importance of these loci, we tested the 190 cis-eQTL SNPs mapped to genes that are differentially expressed between dex-treated cell lines and sham-treated cell lines (Supplementary Material, Table S5) for association with LnPC20 in the CAMP African Americans (AAs) randomized to ICS therapy. We replicated the association in 46 (24%) of the SNPs (Supplementary Material, Table S7). Notably, seven variants near six genes met criteria for genome-wide significance when we applied Bonferroni P-value correction for the Fisher combined P-values for these seven variants (affecting SLFN5, NAPRT1, GPATC2, CDKN1A, TMEM177 and SLC2A9, Table 2, with P <2 × 10−7 for the 240 000 SNP tested).

Table 1.

ORs for cis-eQTL association and SNP-treatment interaction with ICS response analyzed under dexamethasone treatment condition

cis-eQTL association
 
SNP-treatment interaction
 
 Differentially expressed* Not differentially expressed  Differentially expressed* Not differentially expressed 
Associated with airways responsivenessa,190 89 Associated with treatment-responsive differential airways responsivenessa,176 86 
Not associated with airways responsivenessa 1149 948 Not associated with treatment-responsive airways responsivenessa 1163 951 
OR [95% CI] 1.76 [1.35, 2.30]  OR [95% CI] 1.67 [1.28, 2.20]  
cis-eQTL association
 
SNP-treatment interaction
 
 Differentially expressed* Not differentially expressed  Differentially expressed* Not differentially expressed 
Associated with airways responsivenessa,190 89 Associated with treatment-responsive differential airways responsivenessa,176 86 
Not associated with airways responsivenessa 1149 948 Not associated with treatment-responsive airways responsivenessa 1163 951 
OR [95% CI] 1.76 [1.35, 2.30]  OR [95% CI] 1.67 [1.28, 2.20]  

aAirways responsiveness as measured by LnPC20 while on ICS medications over the 4-year-CAMP clinical trial.

*FDR-adjusted P-values <0.05.

Table 2.

Seven genome-wide significant dexamethasone cis-eQTL associated with LnPC20 change during ICS therapy

eQTL SNP Gene chr eQTL P-value Differential expression P-value Caucasian beta Caucasian P-value* AA beta AA P-value* Fisher's combined P-valuea 
rs2671822 SLFN5 17 1.88E−07 2.94E−07 0.31 2.57E−06 0.62 6.51E−06 4.33E−10 
rs3793371 NAPRT1 1.25E−05 2.84E−03 −0.66 1.98E−06 −2.34 4.09E−05 1.97E−09 
rs265125b GPATC2 1.26E−05 2.21E−06 0.40 2.79E−06 0.51 4.27E−04 2.56E−08 
rs4713999 CDKN1A 3.68E−08 1.51E−10 −0.40 8.22E−08 −0.34 2.31E−02 4.01E−08 
rs2272053 TMEM177 1.38E−07 7.85E−08 0.26 2.59E−03 1.43 2.03E−06 1.06E−07 
rs10939602b SLC2A9 4.42E−05 1.06E−02 −0.27 1.28E−04 −0.80 4.62E−05 1.18E−07 
rs6449090b SLC2A9 1.04E−04 1.06E−02 −0.26 2.20E−04 −0.80 4.62E−05 1.98E−07 
eQTL SNP Gene chr eQTL P-value Differential expression P-value Caucasian beta Caucasian P-value* AA beta AA P-value* Fisher's combined P-valuea 
rs2671822 SLFN5 17 1.88E−07 2.94E−07 0.31 2.57E−06 0.62 6.51E−06 4.33E−10 
rs3793371 NAPRT1 1.25E−05 2.84E−03 −0.66 1.98E−06 −2.34 4.09E−05 1.97E−09 
rs265125b GPATC2 1.26E−05 2.21E−06 0.40 2.79E−06 0.51 4.27E−04 2.56E−08 
rs4713999 CDKN1A 3.68E−08 1.51E−10 −0.40 8.22E−08 −0.34 2.31E−02 4.01E−08 
rs2272053 TMEM177 1.38E−07 7.85E−08 0.26 2.59E−03 1.43 2.03E−06 1.06E−07 
rs10939602b SLC2A9 4.42E−05 1.06E−02 −0.27 1.28E−04 −0.80 4.62E−05 1.18E−07 
rs6449090b SLC2A9 1.04E−04 1.06E−02 −0.26 2.20E−04 −0.80 4.62E−05 1.98E−07 

aFisher's combined P-value combines P-values for Caucasian and AA analyses.

beQTL association noted only in the dexamethasone treated cells (i.e. not significantly associated with the sham cells).

*P-values for Caucasians are two-sided, while those in AA are one-sided.

To evaluate the cis-eQTL for treatment-responsive effects, we tested the interaction between the 2484 cis-eQTL and treatment assignment (inhaled corticosteroid or placebo) on LnPC20 in 581 CAMP Caucasians. We identified 266 cis-eQTL associated with a treatment-responsive change in airways responsiveness, of which 176 mapped to differentially expressed genes (Supplementary Material, Table S6), with an OR of 1.67 (1.28–2.20), Table 1). We tested the 176 treatment-responsive eQTL for replication in the CAMP AAs and replicated the interaction in 39 loci (Supplementary Material, Table S8). Seven loci affected expression of five genes (SPATA20, ACOT4, BRWD1, ALG8 and NAPRT1) with a treatment interaction P-value that met genome-wide significance (Table 3). Figure 1 A and B shows the treatment-responsive effects of SNP rs12891009 on ACOT4 gene expression in both populations. As shown, among asthmatics taking placebo, the minor allele is associated with lower LnPC20 (risk), while among those assigned to ICS, the minor allele is associated with higher LnPC20 (protective) in both populations. For the SNPs that did not replicate in the AAs, we checked the distribution of the difference of minor allele frequencies (MAFs) between Caucasians and AAs with the distribution of the difference for SNPs that replicated in the AAs; we did not find a significant difference between the two distributions (data not shown).

Table 3.

Seven genome-wide significant dexamethasone cis-eQTL associated with treatment-responsive change in airways responsiveness

eQTL SNP Gene chr eQTL P-value Differential expression P-value Caucasian beta Caucasian P-value* AA Beta AA P-value* Fisher's combined P-valuea 
rs6504666b SPATA20 17 1.49E−04 1.34E−03 −0.45 2.05E−07 −1.48 1.78E−06 1.08E−11 
rs12891009 ACOT4 14 2.23E−05 9.19E−14 −0.28 1.23E−03 −1.05 2.30E−07 6.50E−09 
rs2037925b BRWD1 21 1.36E−05 5.32E−05 −0.42 3.02E−07 −0.48 3.91E−03 2.54E−08 
rs2836987b BRWD1 21 1.10E−05 5.32E−05 −0.42 2.37E−07 −0.48 5.25E−03 2.67E−08 
rs1380657b SPATA20 17 8.21E−05 1.34E−03 −0.42 1.30E−06 −1.25 1.63E−03 4.46E−08 
rs1144764 ALG8 11 3.95E−12 2.02E−06 −0.28 4.81E−04 −0.82 5.71E−06 5.69E−08 
rs3793371 NAPRT1 1.25E−05 2.84E−03 −0.73 1.98E−05 −2.29 2.39E−04 9.54E−08 
eQTL SNP Gene chr eQTL P-value Differential expression P-value Caucasian beta Caucasian P-value* AA Beta AA P-value* Fisher's combined P-valuea 
rs6504666b SPATA20 17 1.49E−04 1.34E−03 −0.45 2.05E−07 −1.48 1.78E−06 1.08E−11 
rs12891009 ACOT4 14 2.23E−05 9.19E−14 −0.28 1.23E−03 −1.05 2.30E−07 6.50E−09 
rs2037925b BRWD1 21 1.36E−05 5.32E−05 −0.42 3.02E−07 −0.48 3.91E−03 2.54E−08 
rs2836987b BRWD1 21 1.10E−05 5.32E−05 −0.42 2.37E−07 −0.48 5.25E−03 2.67E−08 
rs1380657b SPATA20 17 8.21E−05 1.34E−03 −0.42 1.30E−06 −1.25 1.63E−03 4.46E−08 
rs1144764 ALG8 11 3.95E−12 2.02E−06 −0.28 4.81E−04 −0.82 5.71E−06 5.69E−08 
rs3793371 NAPRT1 1.25E−05 2.84E−03 −0.73 1.98E−05 −2.29 2.39E−04 9.54E−08 

aFisher's combined P-value combines P-values for Caucasian and AA analyses.

beQTL association noted only in the dexamethasone treated cells (i.e. not significantly associated with the sham cells).

*P-values for Caucasians are two-sided, while those in AA are one-sided.

Figure 1.

Clinical association of eQTL identified via dexamethasone-treated gene expression is largely treatment specific. (A) Within CAMP Caucasians, for this variant (rs12891009) identified as a dex-eQTL, there is an interaction eQTL effect on airway responsiveness following treatment with or without inhaled corticosteroids. There are six boxplots. The left three box plots represent individuals taking placebo; the right three box plots correspond to those individuals taking ICS medications. (B) This pharmacogenetic interaction replicates for the eQTL within the CAMP AAs. WT: wild-type; HET: heterozygote; MUT: mutation.

Figure 1.

Clinical association of eQTL identified via dexamethasone-treated gene expression is largely treatment specific. (A) Within CAMP Caucasians, for this variant (rs12891009) identified as a dex-eQTL, there is an interaction eQTL effect on airway responsiveness following treatment with or without inhaled corticosteroids. There are six boxplots. The left three box plots represent individuals taking placebo; the right three box plots correspond to those individuals taking ICS medications. (B) This pharmacogenetic interaction replicates for the eQTL within the CAMP AAs. WT: wild-type; HET: heterozygote; MUT: mutation.

DISCUSSION

In this work, we identified multiple novel genome-wide significant pharmacogenetic loci using a drug-response eQTL model. We found that twice as many clinically significant post-treatment cis-eQTL were located within genes differentially expressed between dexamethasone-treated cells and sham-treated cells in both our main effects and gene-environment analyses. Furthermore, our results demonstrated generalizability, with nearly one-quarter of ICS response eQTL identified in Caucasians also associated in AAs. Combined, 13 unique variants met criteria for genome-wide significance for the clinical association (Tables 2 and 3). Notably, only one SNP (rs2272053) from the ICS specific treatment analysis was also associated (P-value = 0.002) in the 409 subjects randomized to non-steroid therapy. These results indicate that focusing on SNPs that impact gene expression after exposure to a pertinent therapeutic stimulus enables identification of clinically relevant pharmacogenetic loci.

While none of our genome-wide significant eQTL genes has been previously studied in asthma, several are potentially compelling pharmacogenetic loci. Acyl-CoA-thiotransferase 4 (ACOT4) (Fig. 1) is a primary regulator of lipid oxidation in human peroxisomes. Peroxisome proliferator-activated receptor gamma (PPARγ) ligands demonstrate anti-inflammatory effects over and above those of corticosteroids by inhibiting cell growth and inducing apoptosis in human airways smooth muscle cells (15). Nicotinamide phosphoribosyltransferase 1 (NAPRT1) catalyzes the conversion of nicotinic acid (NA) to NA mononucleotide (NaMN) thereby increasing cellular nicotinamide adenine dinucleotide (NAD+) levels and preventing oxidative stress of the cells. Importantly, NA also results in increased cysteinyl leukotriene production (16); cysteinyl leukotrienes augment smooth muscle constriction and airway edema in asthma and are a recognized therapeutic target (17). Therefore, variation in the NA pathway has direct relevance to asthma pathogenesis and treatment response. Figure 2 illustrates the network connecting the five genes in Table 3 and their functionally similar genes from the literature using GeneMANIA (18). By applying DAVID functional annotation, we note that the genes shown in Figure 3 are enriched for metabolic biological processes (Supplementary Material, Table S11).

Figure 2.

Networks constructed by geneMANIA based on the five genes in Table 3 show these five genes are connected based on the literature. The black circles indicate the five genes we detected, while the gray circles are genes in the literature. The edges indicate that there exist evidences from literatures linking the genes connected by the edges.

Figure 2.

Networks constructed by geneMANIA based on the five genes in Table 3 show these five genes are connected based on the literature. The black circles indicate the five genes we detected, while the gray circles are genes in the literature. The edges indicate that there exist evidences from literatures linking the genes connected by the edges.

Figure 3.

Focusing only on FDR-adjusted eQTL from untreated (basal) samples would have missed many clinically relevant eQTL. This figure illustrates this point by using the snp-gene pair (rs2037925, BRWD1) in Table 3. (A) Boxplots of expression level versus genotype based on sham-treated cell lines and dex-treated cell lines, separately. The nominal sham P-value does not meet criteria for significance following correction for multiple comparisons. (B) Boxplots of log2 difference of expression level [log2(Dex/Sham)] versus genotype. WT: wild-type; HET: heterozygote; MUT: mutation. This also demonstrates that eQTL associated with subtle changes in log-difference expression may represent clinically significant loci.

Figure 3.

Focusing only on FDR-adjusted eQTL from untreated (basal) samples would have missed many clinically relevant eQTL. This figure illustrates this point by using the snp-gene pair (rs2037925, BRWD1) in Table 3. (A) Boxplots of expression level versus genotype based on sham-treated cell lines and dex-treated cell lines, separately. The nominal sham P-value does not meet criteria for significance following correction for multiple comparisons. (B) Boxplots of log2 difference of expression level [log2(Dex/Sham)] versus genotype. WT: wild-type; HET: heterozygote; MUT: mutation. This also demonstrates that eQTL associated with subtle changes in log-difference expression may represent clinically significant loci.

We identified few log-difference cis-eQTL (Supplementary Material, Table S9), that is, cis-eQTL associated with dexamethasone-induced changes in gene expression versus sham. Although two of the four (rs10505755, affecting C12orf59; rs962817 affecting ANKRD15) were associated with LnPC20 response in the 172 CAMP Caucasians, this finding was not replicated in CAMP AAs. The rarity of cis-eQTL meeting stringent criteria in differential response to treatment is consistent with the prior studies of environmental perturbation eQTL (1012). For instance, a recent study by Mangravite et al. (19), only six such differential cis-eQTL in LCLs treated with statins were detected. Similarly, Maranville et al. (11), noted just nine cis-eQTL in a study of dexamethasone-induced gene expression changes in HapMap LCLs. Notably, despite differences in passage number and treatment duration, there was significant overlap of their eQTL in our dataset, as previously reported (11). In contrast, Smirnov et al. (10) identified extensive trans-acting variants in cell lines exposed to radiation. Exploring the role of trans-acting log-difference cis-eQTL is an important analysis that is beyond the scope of this paper.

While we focus on the SNPs associated (i.e. FDR-adjusted P-value < 0.05) to gene expressions from dexamethasone-treated cell lines (denoted as dex-eQTL) rather than SNPs associated (i.e. FDR-adjusted P-value < 0.05) to gene expressions from sham-treated samples (denoted as sham-eQTL) in this manuscript, we emphasize that among the 190 dex-eQTL that were associated with LnPC20 and were associated to genes that were differentially expressed between dexamethasone-treated and sham-treated cell lines, only 116 (60%) SNPs were also sham-eQTL. (The P-values and FDR adjusted P-values of eQTL tests for the other 40% (74) SNPs were shown in Supplementary Material, Table S10). Therefore, focusing only on eQTL from untreated (basal) samples would have missed many clinically relevant eQTL, including the majority of the genome-wide significant associations (Tables 2, 3 and Fig. 3). Of the 74 dex-eQTL, 26 (35%) were also nominal eQTL for the difference in expression (log-difference cis-eQTL with P < 0.05), but only at P-values less stringent than the FDR used in the log-difference analysis (Fig. 3 and Supplementary Material, Table S10). Thus, clinically relevant log-difference (i.e. treatment response) cis-eQTL are likely to be missed by focusing on genome-wide thresholds for identification. Despite this, many sham-eQTL were indistinguishable from dex-eQTL in terms of their genotypic effect on gene expression. We reiterate that the sham-eQTL overlapping with our clinically relevant dex-eQTL all mapped to differentially expressed genes. Therefore, it is likely that in these genes the treatment effect on expression presents an additive interaction between the SNP and corticosteroid treatment that influences ICS response (Supplementary Material, Fig. S1).

To assess consistency of our cis-eQTL results with other datasets, we applied the web tool GTEx (Genotype-Tissue Expression) eQTL browser (20). We found that 326 (13.1%) of the 2484 dex-eQTL we detected were previously published in the literature based on Lymphoblastoid tissues and that 269 (10.8%) of the 2484 cis-eQTL have also been detected by other researchers in the brain tissue. Unfortunately, lung tissue data are not yet available in GTEx.

Several potential limitations to our analysis exist. First, our clinical sample sizes were relatively small. However, the longitudinal nature of our trial increased statistical power. Additionally, we were able to identify novel replicable airways responsiveness loci in both our main effects and interaction analyses despite the sample size. Secondly, we tested immortalized B-lymphocytes rather than primary B-cells from trial participants. However, with regard to LCLs, Ding et al. (21) recently reported that 70% of cis-eQTL in LCLs is shared with skin, suggesting that our LCL analysis may also have large overlap with primary cells. Moreover, B-lymphocytes are crucial inflammatory mediators in asthma that produce IgE in response to cytokines and T-cell interactions to propagate allergic asthma (2225). In turn, IgE levels are very highly correlated with PC20, our primary outcome phenotype (26). The fact that the results from differential expression and cis-eQTL analysis in LCL showed enrichment of LnPC20-association SNP further highlights the biological relevance of our findings.

In summary, cis-eQTL derived from steroid-treated cell lines, in conjunction with differential expression, identified multiple novel, genome-wide loci in asthma that would not have been noted via the analysis of basal expression alone. These replicated eQTL represent true pharmacogenetic loci in asthma. It is likely that our method of utilizing both eQTL data and gene expression response to environmental perturbation to identify SNPs associated with clinically important outcomes can be used more broadly to identify novel pharmacogenetic and other gene-environment loci.

MATERIALS AND METHODS

Population

The CAMP Genetics Ancillary Study was approved by each study center’s internal review board, and informed consent/assent was obtained from all participants and their parents. The CAMP clinical trial enrolled 1041 children aged 5–12 years (27,28) with mild-to-moderate asthma, who were randomly assigned to treatment with budesonide (a corticosteroid), nedocromil, or placebo and followed for a mean of 4.8 years. Because there was no difference in treatment effects of nedocromil or placebo, those groups are analyzed jointly in this work. We focus in this analysis on a subset of 581 Caucasian (172 randomized to ICS) and 97 AA (32 randomized to ICS) CAMP subjects for whom genome-wide genotyping (GWAS) is available. The descriptive statistics of these subjects are summarized in Supplementary Material, Table S1.

Airways responsiveness, as measured by the provocative concentration of methacholine causing a 20% decrement in FEV1 (PC20) values were obtained at 0, 8, 20, 32 and 44 months. Airway responsiveness was determined by methacholine testing with the Wright nebulizer technique using methods as previously described (29,30). PC20 was log transformed (LnPC20) for all analyses.

Gene expression arrays

LCLs derived from 151 CAMP Caucasian participants were utilized. Each cell line culture was split into two equal parts: one part was treated with 10−6m dexamethasone (a corticosteroid); the other part was sham treated. After 6h, expression levels were measured using the Illumina HumanRef8 v2 BeadChip (Illumina, San Diego, CA). Each pair of dexamethasone-treated and sham-treated arrays was in the same batch. The resulting gene expression profile data include expression levels for 22 184 gene probes and 2 arrays for each individual. There were 464 arrays in total. In data pre-process stage, we removed 35 technically failed arrays, 65 arrays from subjects with missing clinical variables, 104 arrays without genotype data, 2 arrays without information about the principal components of genotype data, 1 array with the ratio of the 95th percentile of expression level to the 5th percentile <10, and 31 replicated arrays. For each of the two types (dexamethasone-treated and sham treated) of arrays, we kept only the array with the largest median expression level for subjects with replicated arrays. We also removed 49 probes with bad chromosome annotation and 960 probes in X or Y chromosome. We then did vst transformation and quantile normalization for all the remaining 234 arrays together (117 dexamethasone-treated arrays and 117 sham-treated arrays) to reduce the effects of technical noises and to make the distribution of expression level for each array closer to normal distribution. Data pre-processing steps are summarized in the Supplementary Material, Table S2. After filtering, we evaluated 21 175 gene probes (17 370 unique genes). We limited cis-eQTL analysis to 117 Caucasian individuals who were randomized to ICS treatment in the clinical trial and for whom GWAS genotyping data was available.

SNP genotyping

SNP genotyping was performed using Illumina’s HumanHap550v3 Genotyping BeadChip (Illumina, Inc., San Diego, CA). Detailed QC methods have been published previously (31). Briefly, we removed markers for low clustering scores, for flanking sequences that failed to map to a unique position on the HG17 reference genome sequence, monomorphic markers, those with five or more Mendel errors, missing in >5% of subjects, minor allele frequency (MAF) <5%, and mapping to the X or Y chromosome. Among the remaining 431 484 SNPs, the average completion rate was >99%. We assessed genotype reproducibility by plating four subjects on each of 14 genotyping plates. All of these replicates had at least 99.8% concordance. We used Eigenstrat software (32) to adjust for four principal components of population substructure in all analyses.

Statistical analysis

As noted above, each of the 117 CAMP cell lines was split into two equal parts; one exposed to dexamethasone, the other to sham treatment (placebo). We first assessed the effect of dex treatment on gene expression levels by testing if β1 = 0 for linear regression log2(Dex/Sham) = β1 + β2age + β3gender + β4-8PCs + error term.

We then performed cis-eQTL analyses to detect for SNP effects on gene expression, by testing a general linear model in an additive model (i.e. testing if β1 = 0 for log2(Dex) = β0 + β1SNP + β2age + β3gender + β4–8PCs + error term) focusing on log2 expression following dexamethasone treatment. We used this model to test for cis-eQTL (SNPs within 50 kb of a gene tested with the 21 175 gene probes).

We next evaluated the associations of the cis-eQTL with LnPC20 over time in 172 CAMP Caucasian subjects treated with ICS using generalized estimating equations (GEE). In the GEE model, we adjusted for potential confounding factors: gender and longitudinal measurements of age and height (i.e. treating age and height as time-dependant covariates given known change in PC20 with age).

In addition to our main effects model, we performed a SNP-treatment group interaction analysis on the entire cohort of 581 subjects for the eQTL, including the 172 subjects randomized to ICS and 409 subjects randomized to non-steroid therapy in the GEE model.

Given our main hypothesis, we specifically assessed whether cis-eQTL associated with LnPC20 are also more likely to affect genes that are differentially expressed after dexamethasone treatment. Finally, we tested the generalizability of the impact of the cis-eQTL on PC20 by testing whether they also affect methacholine responsiveness in the CAMP AA subjects.

For all analyses except replication tests, we used two-sided tests and adjusted the P-values using Benjamini and Hochberg's method to control multiple testing (33). For clinical association tests, if multiple cis-eQTL SNP were in perfect LD, we tested only one of them for association with lnPC20. A test is claimed as significant if the FDR-adjusted P-value < 0.05. For replication, we used one-sided test and un-adjusted P-value with a P-value threshold of <0.05. The statistical software R was used to perform all analyses. Bioconductor (www.bioconductor.org, last accessed date for www.bioconductor.org is 28 April 2014) packages limma, GGtools, and gee were used to perform gene expression differential analysis, eQTL analyses and GEE analyses, respectively (34,35).

SUPPLEMENTARY MATERIAL

Supplementary Material is available at HMG online.

FUNDING

This work was supported by National Institutes of Health (R01 HL092197, K23 HG003983, R01 NR13391) and U01 HL065899).

ACKNOWLEDGEMENTS

We thank all 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.

Conflict of Interest statement. All authors have no conflict of interests to report. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

REFERENCES

1
NHIS
National health interview survey
2005
Hyattsville, MD
National Center for Health Statistic (NCHS), Centers for Disease Control and Prevention
2
Expert Panel Report 3 (EPR-3): Guidelines for the Diagnosis and Management of Asthma-Summary Report 2007
J. Allergy Clin. Immunol.
 
120
S94
138
National Asthma Education and Prevention Program (2007)
3
Drazen
J.M.
Silverman
E.K.
Lee
T.H.
Heterogeneity of therapeutic responses in asthma
Br. Med. Bull.
 
2000
56
1054
1070
4
Tantisira
K.G.
Hwang
E.S.
Raby
B.A.
Silverman
E.S.
Lake
S.L.
Richter
B.G.
Peng
S.L.
Drazen
J.M.
Glimcher
L.H.
Weiss
S.T.
TBX21: a functional variant predicts improvement in asthma with the use of inhaled corticosteroids
Proc. Natl. Acad. Sci. U. S. A.
 
2004
101
18099
18104
5
Tantisira
K.G.
Lake
S.
Silverman
E.S.
Palmer
L.J.
Lazarus
R.
Silverman
E.K.
Liggett
S.B.
Gelfand
E.W.
Rosenwasser
L.J.
Richter
B.
et al.  
Corticosteroid pharmacogenetics: association of sequence variants in CRHR1 with improved lung function in asthmatics treated with inhaled corticosteroids
Hum. Mol. Genet.
 
2004
13
1353
1359
6
Tantisira
K.G.
Silverman
E.S.
Mariani
T.J.
Xu
J.
Richter
B.G.
Klanderman
B.J.
Litonjua
A.A.
Lazarus
R.
Rosenwasser
L.J.
Fuhlbrigge
A.L.
et al.  
FCER2: a pharmacogenetic basis for severe exacerbations in children with asthma
J. Allergy Clin. Immunol.
 
2007
120
1285
1291
7
Tantisira
K.G.
Damask
A.
Szefler
S.J.
Schuemann
B.
Markezich
A.
Su
J.
Klanderman
B.
Sylvia
J.
Wu
R.
Martinez
F.
et al.  
Genome-wide association identifies the T gene as a novel asthma pharmacogenetic locus
Am. J. Respir. Crit. Care Med.
 
2012
185
1286
1291
8
Tantisira
K.G.
Lasky-Su
J.
Harada
M.
Murphy
A.
Litonjua
A.A.
Himes
B.E.
Lange
C.
Lazarus
R.
Sylvia
J.
Klanderman
B.
et al.  
Genomewide association between GLCCI1 and response to glucocorticoid therapy in asthma
N. Engl. J. Med.
 
2011
365
1173
1183
9
Moffatt
M.F.
Kabesch
M.
Liang
L.
Dixon
A.L.
Strachan
D.
Heath
S.
Depner
M.
von Berg
A.
Bufe
A.
Rietschel
E.
et al.  
Genetic variants regulating ORMDL3 expression contribute to the risk of childhood asthma
Nature
 
2007
448
470
473
10
Smirnov
D.A.
Morley
M.
Shin
E.
Spielman
R.S.
Cheung
V.G.
Genetic analysis of radiation-induced changes in human gene expression
Nature
 
2009
459
587
591
11
Maranville
J.C.
Luca
F.
Richards
A.L.
Wen
X.
Witonsky
D.B.
Baxter
S.
Stephens
M.
Di Rienzo
A.
Interactions between glucocorticoid treatment and cis-regulatory polymorphisms contribute to cellular response phenotypes
PLoS Genet.
 
2011
7
e1002162
12
Grundberg
E.
Adoue
V.
Kwan
T.
Ge
B.
Duan
Q.L.
Lam
K.C.
Koka
V.
Kindmark
A.
Weiss
S.T.
Tantisira
K.
et al.  
Global analysis of the impact of environmental perturbation on cis-regulation of gene expression
PLoS Genet.
 
2011
7
e1001279
13
Huang da
W.
Sherman
B.T.
Lempicki
R.A.
Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources
Nat. Protoc.
 
2009
4
44
57
14
Huang da
W.
Sherman
B.T.
Lempicki
R.A.
Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists
Nucleic Acids Res.
 
2009
37
1
13
15
Patel
H.J.
Belvisi
M.G.
Bishop-Bailey
D.
Yacoub
M.H.
Mitchell
J.A.
Activation of peroxisome proliferator-activated receptors in human airway smooth muscle cells has a superior anti-inflammatory profile to corticosteroids: relevance for chronic obstructive pulmonary disease therapy
J. Immunol.
 
2003
170
2663
2669
16
Saareks
V.
Ylitalo
P.
Mucha
I.
Riutta
A.
Opposite effects of nicotinic acid and pyridoxine on systemic prostacyclin, thromboxane and leukotriene production in man
Pharmacol. Toxicol.
 
2002
90
338
342
17
Drazen
J.M.
Israel
E.
O'Byrne
P.M.
Treatment of asthma with drugs modifying the leukotriene pathway
N. Engl. J. Med.
 
1999
340
197
206
18
Warde-Farley
D.
Donaldson
S.L.
Comes
O.
Zuberi
K.
Badrawi
R.
Chao
P.
Franz
M.
Grouios
C.
Kazi
F.
Lopes
C.T.
et al.  
The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function
Nucleic Acids Res.
 
2010
38
W214
W220
19
Mangravite
L.M.
Engelhardt
B.E.
Medina
M.W.
Smith
J.D.
Brown
C.D.
Chasman
D.I.
Mecham
B.H.
Howie
B.
Shim
H.
Naidoo
D.
et al.  
A statin-dependent QTL for GATM expression is associated with statin-induced myopathy
Nature
 
2013
502
377
380
20
The Genotype-Tissue Expression (GTEx) project
Nat. Genet.
 
2013
45
580
585
21
Ding
J.
Gudjonsson
J.E.
Liang
L.
Stuart
P.E.
Li
Y.
Chen
W.
Weichenthal
M.
Ellinghaus
E.
Franke
A.
Cookson
W.
et al.  
Gene expression in skin and lymphoblastoid cells: refined statistical method reveals extensive overlap in cis-eQTL signals
Am. J. Hum. Genet.
 
2010
87
779
789
22
Oettgen
H.C.
Geha
R.S.
IgE in asthma and atopy: cellular and molecular connections
J. Clin. Invest.
 
1999
104
829
835
23
Tsitoura
D.C.
Yeung
V.P.
DeKruyff
R.H.
Umetsu
D.T.
Critical role of B cells in the development of T cell tolerance to aeroallergens
Int. Immunol.
 
2002
14
659
667
24
Geha
R.S.
Jabara
H.H.
Brodeur
S.R.
The regulation of immunoglobulin E class-switch recombination
Nat. Rev. Immunol.
 
2003
3
721
732
25
Pearlman
D.S.
Pathophysiology of the inflammatory response
J. Allergy Clin. Immunol.
 
1999
104
S132
S137
26
Sears
M.R.
Burrows
B.
Flannery
E.M.
Herbison
G.P.
Hewitt
C.J.
Holdaway
M.D.
Relation between airway responsiveness and serum IgE in children with asthma and in apparently normal children
N. Engl. J. Med.
 
1991
325
1067
1071
27
The Childhood Asthma Management Program (CAMP): design, rationale, and methods. Childhood Asthma Management Program Research Group
Control Clin. Trials
 
1999
20
91
120
28
Long-term effects of budesonide or nedocromil in children with asthma. The Childhood Asthma Management Program Research Group
N. Engl. J. Med.
 
2000
343
1054
1063
29
Cockcroft
D.W.
Killian
D.N.
Mellon
J.J.
Hargreave
F.E.
Bronchial reactivity to inhaled histamine: a method and clinical survey
Clin. Allergy
 
1977
7
235
243
30
Juniper
E.F.
Frith
P.A.
Dunnett
C.
Cockcroft
D.W.
Hargreave
F.E.
Reproducibility and comparison of responses to inhaled histamine and methacholine
Thorax
 
1978
33
705
710
31
Himes
B.E.
Hunninghake
G.M.
Baurley
J.
Rafaels
N.
Sleiman
P.
Strachan
D.
Wilk
J.
Willis-Owen
S.
Klanderman
B.
Lasky-Su
J.A.
et al.  
Genome-wide association analysis identifies PDE4D as an asthma-susceptibility gene
Am. J. Hum. Genet.
 
2009
84
581
593
32
Price
A.L.
Patterson
N.J.
Plenge
R.M.
Weinblatt
M.E.
Shadick
N.A.
Reich
D.
Principal components analysis corrects for stratification in genome-wide association studies
Nat. Genet.
 
2006
38
904
909
33
Benjamini
Y.
Hochberg
Y.
Controlling the false discovery rate: a practical and powerful approach to multiple testing
J. R. Stat. Soc. B
 
1995
57
289
300
34
Carey
V.J.
Morgan
M.
Falcon
S.
Lazarus
R.
Gentleman
R.
GGtools: analysis of genetics of gene expression in bioconductor
Bioinformatics
 
2007
23
522
523
35
Gentleman
R.C.
Carey
V.J.
Bates
D.M.
Bolstad
B.
Dettling
M.
Dudoit
S.
Ellis
B.
Gautier
L.
Ge
Y.
Gentry
J.
et al.  
Bioconductor: open software development for computational biology and bioinformatics
Genome Biol.
 
2004
5
R80