The phosphatidylinositide 3-kinase (PI3K) signaling pathway is a determinant of zileuton response in adults with asthma.

Variable responsiveness to zileuton, a leukotriene antagonist used to treat asthma, may be due in part to genetic variation. While individual SNPs were previously associated with zileuton-related lung function changes, specific quantitative trait loci (QTLs) and biological pathways that may contribute have not been identified. In this study, we investigated the hypothesis that genetic variation within biological pathways is associated with zileuton response. We performed an integrative QTL mapping and pathway enrichment study to investigate data from a GWAS of zileuton response, in addition to mRNA expression profiles and leukotriene production data from lymphoblastoid cell lines (LCLs) (derived from asthmatics) that were treated with zileuton or ethanol (control). We identified 1060 QTLs jointly associated with zileuton-related differential LTB4 production in LCLs and lung function change in patients taking zileuton, of which eight QTLs were also significantly associated with persistent LTB4 production in LCLs following zileuton treatment (i.e., 'poor' responders). Four nominally significant trans-eQTLs were predicted to regulate three candidate genes (SELL, MTF2, and GAL), the expression of which was significantly reduced in LCLs following zileuton treatment. Gene and pathway enrichment analyses of QTL associations identified multiple genes and pathways, predominantly related to phosphatidyl inositol signaling via PI3K. We validated the PI3K pathway activation status in a subset of LCLs demonstrating variable zileuton-related LTB4 production, and show that in contrast to LCLs that responded to zileuton, the PI3K pathway was activated in poor responder LCLs. Collectively, these findings demonstrate a role for the PIK3 pathway and its targets as important determinants of differential responsiveness to zileuton.

Introduction

The development of innovative approaches that integrate genomic, proteomic, and metabolomic data has expanded the number of high-priority candidate genes and biomarkers identified for complex diseases, including asthma [1-5]. These strategies, which combine RNA and/or protein and metabolite expression profiling with genetic association studies, are a breakthrough in the search for the biological mechanisms that underlie complex diseases. SNPs that are tightly linked to or capable of directly influencing gene expression represent QTLs that regulate expression of mRNAs (i.e., eQTLs), or production of metabolites and small molecules (i.e., pQTLs) [6-8], respectively. Through QTL mapping, variation in molecular, transcriptomic, and phenotypic traits can be correlated with variation in DNA sequences [7, 9–11]. Correlation of molecular and phenotypic traits can also be informative for inferring molecular networks that contribute to asthma pathogenesis or drug treatment responses.

In this study, we applied integrative methods to investigate the relationship between patient genotype and treatment response to zileuton, a leukotriene modifying agent [12-14]. Leukotrienes are synthesized from arachidonic acid (AA) by 5-lipoxygenase (5-LO) [15, 16]. After export from the cell, leukotrienes modulate important immunological activities including neutrophil activation, chemotaxis, eosinophil migration and smooth muscle contraction [16, 17]. Zileuton exerts its effects through inhibition of 5-LO-mediated production of leukotrienes, in particular leukotriene B4 (LTB4) [18]. Clinical response to leukotriene modifiers, including zileuton, is highly variable and repeatable within individuals, suggesting that polymorphic variation in genes in the leukotriene pathway may play a role in zileuton response [19-21]. In particular, SNPs that are significantly associated with both treatment-related clinical outcomes, as well as changes in leukotriene levels may represent important genetic biomarkers.

We hypothesized that genetic variation at the level of unique biological pathways is associated with variability in lung function in response to zileuton. In order to identify sets of genes and pathways that were associated with zileuton response, we first performed a pathway-based analysis using GWAS data obtained from a recent genome-wide investigation of zileuton response [22]. The advantage of using data identified through GWAS to interrogate biological pathways is that this method enriches for potentially functionally important candidate genes and SNPs that may not have been individually strongly associated with clinical response phenotypes. Having identified a set of candidate genes and pathways related to clinical response to zileuton, we then sought to identify specific QTLs associated with changes in LTB4 concentrations post-treatment, as well as changes in global mRNA expression. We next integrated pathway enrichment and QTL information to identify potentially causal candidate genes. Finally, we experimentally validated the expression of selected candidate pathway genes, and determined the activation status of a key pathway of interest.

Materials and methods

Genomic datasets

Clinical cohort

Summary statistics, genome-wide genotype and phenotypic information from a recent GWAS [22] investigating two placebo-controlled clinical trials of zileuton response (309 subjects) were available for interrogation. For this study, 226,890 SNPs with concordant effect size estimates (β coefficient sign) in both the discovery and replication GWAS [22] were pursued for pathway analysis.

Cell lines

Lymphoblastoid cell lines (LCLs) were available from 800 probands from the Childhood Asthma Management Program (CAMP), a 4 year National Heart Lung and Blood Institute-sponsored clinical trial comparing the efficacy of budesonide vs. nedocromil vs. placebo in 1041 patients [23, 24]. Genotypes, microarray expression data, and LTB4 concentrations from 345 LCLs derived from these patients were used to generate the QTL datasets.

Approach to hypothesis testing

Statistical approaches for hypothesis testing for each of the core analyses performed in this study are described in Supplemental Table 1.

Pathway analysis methods

Gene-based pathway analysis

For gene-based pathway analysis, genes were analyzed for their enrichment in pathways from KEGG, Pathway Commons and WikiPathways. Query genes were analyzed relative to all genes annotated in the human genome using the hypergeometric test for evaluation of enrichment, and correcting the resulting P-values for multiple testing by the Benjamini and Hochberg method [25]. The resulting pathways were then ranked according to the lowest P-value (or FDR q-value).

SNP-based pathway analysis

For SNP-based pathway analysis, “Improved Gene Set Enrichment Analysis for Genome-Wide Association Studies” (iGSEA-4-GWAS) [26] (http://gsea4gwas.psych.ac.cn/), was used to generate pathways from SNP and P-value data. The iGSEA-4-GWAS approach applies an improved gene set enrichment analysis (i-GSEA) to GWAS data (SNP P-value, odds ratio or other statistics) using SNP label permutation instead of phenotype to enrich for gene sets with high proportion of significant genes [26]. SNP rsIDs and corresponding –log10 P-values were mapped to +/− 500 KB from genes annotated to canonical pathways in the Molecular Signatures Database (MSigDB). A nominal iGSEA P value cutoff of 0.01 was specified, and FDR was used to correct for multiple testing [26]. Pathways meeting a significance threshold of FDR < 0.25 (to denote the confidence of ‘possible’ or ‘hypothesis’) were prioritized.

Cell studies and functional validation methods

Cell culture and incubation experiments

Cellular analyses were performed using lymphoblastoid cell lines (LCLs) that were derived from CAMP. For this study, 100 paired LCLs were cultured in growth medium (RPMI 1640, 5% FBS, 1% L-Glutamine, and 1% antibiotic–antimycotic) in a T25 flask in a 37 °C, 5% CO2 incubator. The cells were grown to a density of ~200,000 viable cells/mL. For stimulation experiments, 10–13 million cells were resuspended in 32 ml (~350 K/ml) regular growth medium and 2.5 ml were dispensed into each well of two 6-well culture plates (~850 K/well). Three wells of each 6- well plate were treated with 10 mM arachidonic acid or 2.5 uL of 100% ethanol. After incubating for 12 hours, the supernatant was removed and cells were treated with 1 µM zileuton (ZIL) or 2.5 uL of 100% ethanol for 6 hours. After treatment, the cells were harvested and spun down for stimulation. Cell stimulations were performed by resuspending cell pellets in phosphate buffered saline (PBS ( + Ca++, + Mg++)) containing 0.1 mM calcium ionophore and 4 mM arachidonic acid and incubated at 37̊ °C for 5 minutes. Cells were placed on ice and supernatants were harvested for quantification of LTB4 concentration using the LTB4 Express EIA kit (Caymen Chemical), while the cell pellets were lysed for RNA extraction using the RNeasy Mini Kit (Qiagen).

LTB4 enzyme-linked immunosorbent assay (ELISA)

An LTB4-specific ELISA was performed using a commercial Leukotriene B4 Enzyme Immunoassay Kit, according to the manufacturer’s instructions. LTB4 concentrations in each sample were determined using fluorimetry at 405 nm wavelength, and normalized to a standard curve.

Microarray analysis

For microarray, 200 ng RNA specimens were amplified using Illumina’s TotalPrep RNA Amplification kit (Ambion) according to the manufacturer’s protocol, and 750 ng of each amplified specimen were hybridized onto a Sentrix BeadChip Array for Gene Expression HumanRef12_V4, according to the manufacturer’s protocol. The BeadChip Array was scanned using Illumina’s BeadArray Reader. Quality control analysis of array data was performed using Illumina’s GenomeStudio v2009.1 Software, Gene Expression Module version 1.1.1 and iCheck_0.1.4.

Western analysis

For Western blot experiments, cell cultures were harvested after 24 hours of treatment (as described above) and lysed using RIPA buffer with protease and phosphatase inhibitors added. Cell lysates were prepped under reducing conditions per the Life Technologie’s iBind™ Western System protocol. Western blots were either probed with a mouse anti-PI3K primary antibody (Abcam), a rabbit anti-AKT primary antibody (Cell Signaling Technology) or a rabbit anti-phospho-AKT primary antibody (Cell Signalling Technology). AKT and phospho-AKT blots were probed with an HRP conjugated goat anti-rabbit IgG secondary antibody (Abcam), while PI3K blots were probed with an HRP conjugated goat anti-mouse IgG secondary antibody (Life Technologies). β-actin was used as a protein loading control and was probed using a mouse anti- β-actin primary antibody (Life Technologies) and an HRP conjugated goat anti-mouse IgG secondary antibody (Life Technologies). Results were visualized using SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific) and captured using a G:Box (Syngene). Quantitative analysis was performed using GeneTools image analysis software (Syngene).

Quantitative trait loci (QTL) analyses

LTB4 QTL (pQTL) Analysis

LTB4 concentrations from 100 LCLs were obtained as described above. A linear regression model was used to associate individual SNP genotypes (>500,000 SNPs) from each sample with LTB4 concentration values, in corresponding zileuton-treated and ethanol-treated samples. PLINK v.1.0.7 was used to exclude SNPs with >10% missing data, a minor allele frequency <5%, and a paired r 2 > 0.9, resulting in 463,791 SNPs remaining for analysis. The effect of each locus on LTB4 levels was analyzed using PLINK v.1.0.7, assuming an additive genetic model, and adjusting for sex and race. Significance of association with LTB4 concentration values was assessed for each SNP using the Chi-square test, and a Bonferroni corrected significance threshold was set at 1.1 × 10−07 for 463,791 tests. SNPs with nominally significant associations (P < 0.05) in the ethanol-treated group were excluded from further analysis. For comparative analysis in samples from a relevant clinical population taking zileuton, a subset of 22,505 pQTLs that were unique to the zileuton treatment group, and that were at least nominally significant based on association P-values < 0.05, were also investigated for overlap using data from a recent GWAS of zileuton response [22]. SNPs that were at least nominally significant for both pQTL and GWAS analyses (pQTL and GWAS P values < 0.05) were selected for pathway analyses.

Expression QTL (eQTL) analysis

To evaluate the effect of zileuton treatment on gene expression, we evaluated 345 paired LCLs treated with zileuton or ethanol. After stimulation with arachidonic acid, each pair was treated with zileuton or ethanol (sham), as described above. The expression levels of 47,131 gene probes were measured using Illumina’s HumanHT-12_V4 chip, and 36,993 gene probes were retained after quality control procedures. The probe expression levels were log2-transformed and quantile normalized, and the log2 difference in expression level between zileuton treated and sham treated cell lines was obtained for each sample. Paired t-test was used to evaluate differential gene expression in samples treated with zileuton vs. sham, and the false discovery rate (FDR) was applied to correct for multiple testing. Of 36,993 probes, 5123 probe sets achieved nominal significance (P < 0.05) based on paired t-test, of which 918 probes surpassed the FDR threshold of 0.05. We next extracted genotypes (505,066 SNPs) from 345 samples, using the genome-wide genotype data from the subjects from whom the individual LCLs were derived (CAMP dataset). SNPs that failed tests for Hardy-Weinberg, missingness, or with minor allele frequency < 0.05 were excluded from the final data set, resulting in 459,261 remaining SNPs. No population substructure was detected by EIGENSTRAT. To identify eQTLs, we used the MatrixEQTL R package to test the association of SNP genotype with the log2 difference in expression for a gene probe, adjusting for age, gender, and the top four principal components of the expression data. The human genome version 19 annotation for both genes and SNPs was used. For each gene probe, all SNPS within a 50 KB range from flanking genic regions wherein a given probe resided were used to test for cis-effects. All SNPs outside this range, and on other chromosomes, were used to test for trans-effects. To increase the scope of potential inputs that could better inform pathway analyses, with a priori knowledge that we could generate functional models to validate the subsequent results, we chose a raw P value cut-off of 1 × 10−08, and applied a FDR < 0.05 as the multiple test correction threshold.

Results

Pathway analysis of zileuton response GWAS

An overview of the study design is shown in Fig. 1. To identify molecular pathways related to zileuton clinical response, we first evaluated existing data from a recent GWAS of zileuton response in 309 asthmatics (Abbott-Zileuton GWAS) [22]. Zileuton response was defined as a change in lung function, as measured by forced expiratory volume in one second from baseline to 8 weeks of zileuton therapy (∆FEV1) [22]. Data from 226,890 SNPs having the same direction of effect (β) in the discovery and replication GWAS were evaluated using a SNP-based pathway analysis and gene set enrichment analysis. A total of 212,364 SNPs were annotated to 15,117 genes and mapped to 523 gene sets, resulting in four enriched pathways (FDR < 0.25): “PLC PATHWAY”, “CASPASE PATHWAY” and “SIG IL4 RECEPTOR IN B LYMPHOCYTES” and “INNATE IMMUNE RESPONSE” (Table 1)).

Fig. 1

Overview of study design

Table 1

Pathways enriched in GWAS of zileuton response

Pathway rank	Pathway name	Pathway description	Pathway source	Gene symbol	SNP ID	Association logP value
1	PLC pathway	Phospholipase C hydrolyzes the membrane lipid PIP2 to DAG to activates protein kinase C, and IP3, causing calcium influx	BioCarta	PIK3R1	rs1344367	4.03
				PLCB1	rs2327120	3.09
				VAV1	rs347033	2.88
				PRKCA	rs228885	2.47
				PIK3CA	rs1607237	2.34
2	CASPASE pathway	Caspases are cysteine proteases active in apoptosis	BioCarta	LMNB1	rs3763156	5.47
				CYCS	rs10239907	4.05
				APAF1	rs12581119	4.01
				LMNB2	rs4806849	3.60
				GZMB	rs8192917	3.31
				LMNA	rs9919256	2.43
3	SIG IL4 Receptor in B LYMPHOCYTES	Genes related to IL4 receptor signaling in B lymphocytes	UCSD Signaling Gateway	PIK3R1	rs1344367	4.03
				JAK3	rs7250423	3.59
				IRS1	rs7574538	3.38
				IL4R	rs9944340	3.36
				BCL2	rs8096380	3.13
				PDK1	rs9287957	2.96
				SOCS1	rs408665	2.51
				PIK3CA	rs1607237	2.34
				IRS2	rs9521522	2.33
				AKT3	rs3006928	2.32
				JAK1	rs12743599	2.15
				GSK3B	rs3732361	2.05
4	Innate immune response	Innate immune responses	AMIGO	CD1D	rs859005	1.35
				SFTPD	rs2243639	1.82

Association logP value corresponds to the adjusted, log-transformed GWAS P value for the SNP used in the pathway database query

QTL mapping

QTL identification and annotation

In an effort to characterize variability in the leukotriene production pathways that could potentially underlie differential responsiveness to zileuton, we evaluated global gene expression and/or LTB4 production in LCLs treated with zileuton or ethanol (control) in paired LCLs derived from asthmatic patients, and mapped QTLs for both treatment groups. LCLs were first stimulated with arachidonic acid (AA) to stimulate LTB4 production, and were then treated with either zileuton or ethanol. For pQTL analyses, LTB4 concentrations in 100 paired LCLs were determined by ELISA. As expected, following activation with AA, LCLs treated with zileuton showed a significant but variable decrease in LTB4 production ([LTB4]mean_etoh = 270 ± 271 ng/mL vs. [LTB4]mean_zileuton = 74.5 ± 66.9 ng/mL; P < 0.0001) (Supplemental Fig. 1A). We next integrated genome-wide genotypes (463,791 QC-validated SNPs) from each of the LCL samples with LTB4 concentrations from both treatment groups to identify zileuton-related pQTLs. A set of 22,505 SNPs specific to the zileuton-treated group were nominally associated (unadjusted P < 0.05) with differential LTB4 production, of which eight SNPs remained significant following correction for multiple testing; the top 100 pQTLs, ranked by their P-values, are listed in Supplemental Table 2. Of the significant pQTLs, all were associated with persistent LTB4 production despite treatment with zileuton (i.e., poorer responsiveness) (Supplemental Table 2). Two of the top pQTLs (rs10743992 and rs6485597) were annotated to introns within cAMP responsive element binding protein-like 2 (CREBL2) and PR domain-containing protein 11 (PRDM11) genes, while a third SNP (rs2823730) was found within a long intergenic non-protein coding RNA, LINC00478 (Supplemental Table 2).

For eQTL mapping, 345 paired LCLs were treated with AA, followed by treatment with zileuton or ethanol as described above and in Methods, and gene expression levels for both treatment groups were measured using microarray. We identified 262 trans-eQTLs that met an unadjusted significance threshold of P < 1 × 10−08. The top four eQTLs also met FDR < 0.05. The 262 eQTLs regulated the expression of 161 unique genes, of which 132/262 eQTLs (50.3%) were annotated to 81 genes. The top 50 eQTLs are listed in Supplemental Table 4, and the four eQTLs meeting FDR < 0.05 are shown in Figs. 2b–e. Four of the genes in which the 262 eQTLs resided were also significantly differentially expressed in LCLs following zileuton treatment (FDR < 0.05). A quantile-quantile (QQ) plot (Fig. 2a) of the log-difference in gene expression related to zileuton vs. ethanol treatment revealed a significant deviation from the expected for the trans-eQTL analysis (Figs. 2b–e, boxplots). Both rs12762462 and rs12251868 were annotated to G protein-coupled receptor 158 (GPR158), and are predicted to regulate the expression of selectin L (SELL) (Supplemental Table 4). rs10495426 is an intergenic SNP that is predicted to regulate metal response element binding factor 2 (MTF2) (Supplemental Table 4). rs12414020 was annotated to neuregulin 3 (NRG3) and is predicted to regulate the expression of galanin and galanin message-associated peptide (GMAP) prepropeptide (GAL). In LCLs, all of these eQTLs significantly reduced the expression of the genes they are predicted to regulate following treatment with zileuton (Figs. 2b–e). However, these signals appeared to be driven by a small number of subjects who carried homozygous variant genotypes. We identified 10 individuals who carried the homozygous variant genotypes for the four eQTLs (three for rs10495426, three for rs12762462 and rs12251868 (the same individuals for both), and four for rs12414020) whose removal resulted in P-values that remained nominally significant, but were no longer significant after controlling for multiple testing (FDR < 0.05) (data not shown).

Fig. 2

eQTLs for Zileuton Response

. a The Q–Q plot shows the relationship between observed vs. expected (diagonal; gray line) P-values for 735,488 local (red line/dots) and 13,009,811,306 distant (blue line/dots) gene-SNP P values. Jittered boxplots b-e show individual data points (red dots) and medians (±IQR) of log2 mRNA expression values (zileuton:ethanol) (y-axis), stratified by genotype category (“0” = reference allele; “1” = heterozygous variant; “2” = homozygous variant) (x-axis): b Expression of probe 6940358 (gene symbol = SELL) by rs12762462 genotype; c Expression of probe 2760255 (gene symbol = MTF2) by rs10495426 genotype; d Expression of probe 6940358 (gene symbol = SELL) by rs12251868 genotype; e Expression of probe 60452 (gene symbol = GAL) by rs12414020 genotype

QTL integration and pathway analysis

There were 22,505 nominally significant pQTLs that were unique to the zileuton treatment group; these were evaluated to determine whether any were also associated with zileuton-related changes in lung function, using data from a GWAS of zileuton response in 309 asthmatics [22]. A total of 1060 SNPs were nominally associated with differential LTB4 production, as well as zileuton-related 8-week ∆FEV1 based on pQTL and GWAS P values < 0.05. The top 50 pQTLs are shown in Supplemental Table 3. Pathway analysis of the 1060 pQTLs identified two major ‘hubs’ inclusive of the top 10 significantly enriched pathways: phosphatidyl inositol signaling (via phosphatidylinositide 3-kinase  (PI3K)) and catenin signaling. The phosphatidyl inositol signaling (via PI3K) hub included pathways for “FCγR-mediated phagocytosis”, “FcεR1 signaling pathway”, “Focal adhesion”, “Neurotrophin signaling pathway”, “Regulation of actin cytoskeleton”, “Pathways in cancer” and “Type II diabetes mellitus”, while the catenin signaling hub included pathways for “Adherens junction”, “Endometrial cancer”, “Pathways in cancer”, “Arrhythmogenic right ventricular cardiomyopathy (ARVC)” (shown in Table 2). These pathways comprise 31 unique genes, of which PI3K isoforms (PIK3R1 and PIK3R5) were most highly enriched relative to all other genes (Supplemental Fig. 2).

Table 2

Pathways enrichment analysis of pQTLs for LTB4 production in asthmatic LCLs

Pathway rank	Pathway name	Gene symbol	SNP ID	Adjusted pathway P value
1	Fc gamma R-mediated phagocytosis	PRKCE	rs2176347	1.41E-05
		PRKCE	rs1878384
		VAV1	rs347033
		VAV3	rs6680436
		WASF3	rs9512311
		PIK3R5	rs516434
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		PRKCD	rs17052826
		DOCK2	rs2161148
		SYK	rs12550873
		SYK	rs290212
		SYK	rs4744505
		SYK	rs1333630
2	Fc epsilon RI signaling pathway	PRKCE	rs2176347	2.26E-05
		PRKCE	rs1878384
		VAV1	rs347033
		VAV3	rs6680436
		PIK3R5	rs516434
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		PRKCD	rs17052826
		MAPK14	rs2237093
		SYK	rs12550873
		SYK	rs290212
		SYK	rs4744505
		SYK	rs1333630
3	Focal adhesion	VAV1	rs347033	5.52E-05
		VAV3	rs6680436
		PIK3R5	rs516434
		ITGA6	rs7584977
		ITGA6	rs7596259
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		RELN	rs362732
		COL4A1	rs631428
		ACTN1	rs2300875
		ITGAV	rs3738919
		EGFR	rs13244925
		SHC2	rs740871
4	Adherens junction	WASF3	rs9512311	8.58E-05
		ACTN1	rs2300875
		PTPRB	rs2717418
		PTPRM	rs2005039
		EGFR	rs13244925
		TCF7L2	rs3814573
		CTNNA2	rs383452
		CTNNA2	rs1816909
5	Endometrial cancer	PIK3R5	rs516434	0.0001
		EGFR	rs13244925
		TCF7L2	rs3814573
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		CTNNA2	rs383452
		CTNNA2	rs1816909
		AXIN2	rs4541111
6	Neurotrophin signaling pathway	MAP2K5	rs3784716	0.0003
		PIK3R5	rs516434
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		PRKCD	rs17052826
		NTRK3	rs2279409
		NTRK3	rs3784410
		MAPK14	rs2237093
		NTRK2	rs4486281
		SHC2	rs740871
7	Regulation of actin cytoskeleton	VAV1	rs347033	0.0003
		VAV3	rs6680436
		PIK3R5	rs516434
		ITGA6	rs7584977
		ITGA6	rs7596259
		FGF18	rs4620037
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		ACTN1	rs2300875
		ITGAV	rs3738919
		F2R	rs153311
		EGFR	rs13244925
8	Pathways in cancer	GLI3	rs2299142	0.0004
		GLI3	rs846266
		GLI3	rs846335
		PIK3R5	rs516434
		ITGA6	rs7584977
		ITGA6	rs7596259
		FGF18	rs4620037
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		AXIN2	rs4541111
		COL4A1	rs631428
		ITGAV	rs3738919
		ETS1	rs7115613
		ETS1	rs4628683
		EGFR	rs13244925
		TCF7L2	rs3814573
		CTNNA2	rs383452
		CTNNA2	rs1816909
9	Arrhythmogenic right ventricular cardiomyopathy (ARVC)	ACTN1	rs2300875	0.0005
		ITGAV	rs3738919
		CACNB2	rs1571787
		ITGA6	rs7584977
		ITGA6	rs7596259
		TCF7L2	rs3814573
		CTNNA2	rs383452
		CTNNA2	rs1816909
10	Type II Diabetes mellitus	PRKCE	rs2176347	0.0006
		PRKCE	rs1878384
		PIK3R5	rs516434
		CACNA1E	rs1999838
		PIK3R1	rs3730082
		PIK3R1	rs16897561
		PRKCD	rs17052826

We next evaluated whether any of the 262 eQTLs were also present (overlapping) within the set of 22,505 pQTLs related to LTB4 production. We identified 22 overlapping eQTL/pQTL relationships, including three of the top-ranked eQTLs (rs12762462, rs10495426 and rs12251868 (Supplemental Table 4)); 18/22 eQTLs/pQTLs were annotated to 7 genes (Supplemental Table 5). However, none of the 262 eQTLs were also included among the set of 1060 pQTLs that were jointly associated with LTB4 production and zileuton clinical response; as such, they could represent unique determinants of zileuton pharmacological response vs. clinical response.

Pathway analysis of the 262 eQTLs identified seven significantly enriched pathways (Table 3). Notably, these pathways also coincided with those for the pQTLs (Table 3, Supplemental Table 4). Among the pathways, several genes that directly or indirectly interacted with the PIK3 signaling pathway hub were also present, e.g., calcium voltage-gated channel auxiliary subunit beta 4 (CACNB4) and fibroblast growth factor (FGF) isoforms within the “MAPK Signaling Pathway”, and F GF isoforms and integrin subunit beta 6 (ITGB6) within the “Regulation of the actin cytoskeleton” pathway, which act upstream of the phosphatidylinositol signaling pathway (Table 3).

Table 3

Pathway enrichment analysis of zileuton-related eQTLs

Pathway rank	Pathway name	Gene symbol	SNP ID	Adjusted pathway P value
1	Arrhythmogenic right ventricular cardiomyopathy (ARVC)	ITGB6	rs13403826	0.0036
		CDH2	rs8087457
		CDH2	rs643555
		CACNB4	rs16830500
2	Purine metabolism	GDA	rs10781076	0.0162
		FHIT	rs6781028
		PDE4B	rs12032117
3	Endocytosis	SH3GL3	rs2637978	0.017
		FLT1	rs10507385
		RAB31	rs2284137
4	Regulation of actin cytoskeleton	MYH9	rs9610489	0.017
		FGF14	rs8000905
		FGF14	rs1336724
		FGF14	rs7989386
		FGF14	rs1928510
		ITGB6	rs13403826
5	MAPK signaling pathway	FGF14	rs8000905	0.018
		FGF14	rs1336724
		FGF14	rs7989386
		FGF14	rs1928510
		CACNB4	rs16830500
		PLA2G6	rs738322
		PLA2G6	rs132985
6	Hypertrophic cardiomyopathy (HCM)	ITGB6	rs13403826	0.018
		CACNB4	rs16830500
7	Dilated cardiomyopathy	ITGB6	rs13403826	0.018
		CACNB4	rs16830500
8	Focal adhesion	ITGB6	rs13403826	0.05
		FLT1	rs10507385

To determine the functional importance of these pathways in LCLs, the expression of 68 genes that were collectively identified through pathway analyses (listed in Tables 1–3 ) was validated by microarray analysis of zileuton treated vs. ethanol treated LCLs. Of these, 8 genes were differentially expressed (zileuton vs. ethanol) at a nominal significance threshold, with one gene, suppressor of cytokine signaling 1 (SOCS1), remaining significant following multiple test correction (Supplemental Table 7).

Experimental validation of candidate pathway and genes

Pathway prioritization

A pathway analysis of 1060 joint pQTLs had identified phosphatidyl inositol signaling via PI3K as a major hub, where PI3K isoforms were more frequently enriched among pathway genes (Table 2; Supplemental Fig. 2). In addition, by integrating pQTL and eQTL data, we also identified several genes that interact with PIK3 (Table 3 and Supplemental Table 4), namely phosphodiesterase 4B (PDE4B), fragile histidine triad (FHIT), myosin heavy chain 9 (MYH9), cadherin 2 (CDH2), and others. Furthermore, we show that among the 8 differentially expressed pathway genes, PIK3CA expression was reduced, albeit at a nominal significance level (Supplemental Table 7). As the collective analyses pointed toward a rationale for prioritizing the PI3K signaling pathway (which was further emphasized due to its relevance as a potential therapeutic target for asthma [27, 28]), we hypothesized that differential zileuton responsiveness could be related to changes in the expression and activation of genes within this pathway. We therefore next validated the PI3K pathway using a cellular model of zileuton response.

Pathway validation

To validate the PI3K pathway in a clinically relevant cellular model, we first selected LCLs that had demonstrated different magnitudes of treatment responses to zileuton, based on LTB4 production. Four LCLs that had demonstrated significant down-regulation of LTB4 production following zileuton treatment (mean fold change = 6.8 (termed “good responders”)) and four LCLs that demonstrated less than 2-fold change in LTB4 concentrations following zileuton treatment (mean fold change = 1.2 (termed “poor responders”)) were selected for profiling (Supplemental Fig. 1B).

A consequence of PI3K pathway activation is an increase in expression of phosphorylated AKT (phospho-AKT), relative to the un-phosphorylated AKT protein [29, 30]. To evaluate the activation status of this pathway, we compared the expression of AKT, phospho-AKT and PI3K proteins following treatment with fetal bovine serum (FBS), zileuton, or sham (ethanol), in the two groups of LCLs possessing differential LTB4 production following zileuton treatment (i.e., ‘good’ and ‘poor’ responders). For both LCL groups, treatment with fetal bovine serum (FBS), a pathway activator (Fig. 3, red bars), increased expression of all three proteins relative to sham (Fig. 3, green bars), with the greatest increase in AKT and phospho-AKT expression values in the ‘poor’ responders (Fig. 3). Compared with sham treatment, zileuton treatment did not significantly increase expression of phospho-AKT in the ‘good’ responder LCLs (Fig. 3, blue bars). In contrast, relative to ‘good’ responders, expression of phospho-AKT was selectively increased by 64% in the ‘poor’ responder group following zileuton treatment (P = 0.037) (Fig. 3). Otherwise, the expression levels of AKT, phospho-AKT and PI3K did not significantly differ across treatment groups (Fig. 3). These results suggest that, compared with LCLs that responded well to zileuton inhibition, the PI3K pathway is activated in ‘poor’ responder LCLs, at basal levels, as well as in response to zileuton treatment.

Fig. 3

Validation of PIK3CA pathway activation in asthmatic LCLs with differential responses to zileuton

. The barplots show normalized mean expression levels ± SEM for AKT (left plot), phospho-AKT (center plot) and PIK3CA (right plot) in ‘good’ and ‘poor’ responders (labels, x-axis), by treatment category (FBS = red; Sham (ethanol) = green; Zileuton = blue)

Differential gene expression related to zileuton treatment

The top 50 differentially expressed probe sets (genes) are listed in Supplemental Table 6. Among the 918 probe sets, PIK3CA (T-statistic = −2.0; P = 0.044) and AKT1 (T-statistic = −3.2; P = 0.0015) demonstrated an overall reduction of expression in the zileuton treated LCLs relative to sham, across 345 LCLs. However, expression of these genes from the subset of ‘good’ and ‘poor’ responder LCLs was congruent with protein expression values for zileuton vs. ethanol treated LCLS (Supplemental Fig. 3). Notably, the mRNA and protein expression of AKT (zileuton:ethanol) was reduced in both ‘good’ and ‘poor’ responder groups, while the expression of PIK3CA (zileuton:ethanol) was reduced only in the ‘good’ responders (Supplemental Fig. 3).

Discussion

We applied pathway analysis in an integrative approach, combined with experimental validation in a relevant cellular model, to identify candidate genes and pathways underlying clinical and pharmacological (cellular) response to zileuton, a leukotriene inhibitor prescribed for asthma treatment. Using data from a recently published GWAS of zileuton response in asthmatics [22], we first utilized an improved gene set enrichment analysis (i-GSEA) algorithm, which accounts for the combined effects of SNPs/genes with modest effect sizes, and identified four significantly-enriched pathways containing 23 unique genes that collectively represent coordinated pathways related to lipid mediator signaling, leading to cellular growth and survival. Due to the high level of cross talk among these pathways, multiple genes with important regulatory roles in lipid mediator signaling appeared in more than one pathway (e.g., PIKC3A), suggesting that they may be biologically relevant for zileuton treatment responses. The most highly ranked pathway (FDR = 0.022, Table I) comprised phosphoinositide 3-kinases (PI3K) lipid signaling molecules. The PI3K family members phosphorylate the inositol rings of phosphoinositide lipids, influencing vesicle trafficking, cell proliferation, and migration [31]. Akt is also an important co-regulator of PI3K signaling leading to leukotriene production, and the activity of PLA2 is necessary for conversion of AA to 5-LO substrates. Together, these results suggest that key regulatory genes within molecular signaling cascades that regulate the production and signaling of lipid mediators are related to zileuton response in patients.

We subsequently replicated these pathways and identified a set of SNPs that are strongly associated with zileuton-related LTB4 production in LCLs derived from asthmatic patients (pQTLs). Among the top SNP associations that surpassed FDR, we identified rs10743992 and rs6485597 in  cAMP responsive element binding protein-like 2 ( CREBL2) and PR/SET domain 11 (PRDM11) that were associated with a dramatic increase in LTB4 production following zileuton treatment (Table 2). Both rs6485597 and rs10743992 are in complete LD with regulatory and non-synonymous coding SNPs that are also present within these genes, suggesting a potential functional role for CREBL2 and PRDM11 (data not shown). Through pathway analysis of pQTL data, we identified numerous genes within the CREB signaling pathway, including CREB1 (Table 2). These data suggest that CREB-mediated signaling events are involved in zileuton pharmacological response.

Pathway analysis revealed an enrichment of pathways for lipid signaling and cellular migration in response to immune activation. From the 68 genes identified by pathway analysis (shown in Tables 1–3) eight were significantly differentially expressed in zileuton vs. ethanol treated LCLs (Supplemental Table 7). Included with these, as well as a majority of the pathways in Tables 1–3, is PIK3CA, a key mediator of phospholipid signaling that is crucial to the production of leukotrienes. PIK3CA encodes the 110 kDa catalytic subunit of PI3K and relies on ATP to phosphorylate PtdIns, PtdIns4P and PtdIns [4, 5]P2. These phosphorylated lipids activate downstream targets that in turn activate transcriptional regulators such as CREB, leading to the transcription of genes involved in cell proliferation and survival. We also identified PLC activity as a shared pathway in both pQTL analyses. PLC directly affects LTB4 production through conversion of PI to AA, the substrate of 5-LO and precursor of LTB4. Together, these data suggest that PIK3CA pathway activity may represent a functional network that determines the pharmacological and clinical response to leukotriene modifiers. Furthermore, the production of LTB4 in zileuton-treated cells demonstrated high variability (Supplemental Fig. 1). We subsequently validated the PIK3CA pathway in asthmatic LCLs with differential production of LTB4 in response to zileuton. We found that there was significant expression of the activated target of PI3K, phospho-AKT, which was selectively upregulated only in the ‘poor’ responder group (Fig. 3). Relative to ethanol, zileuton-treated LCLs also showed increased mRNA and protein expression of PIK3CA, in ‘poor’ responders only (Supplemental Fig. 3). Collectively, these data suggest that zileuton responses can be regulated by PI3K pathway activity, and that the pathway is differentially activated in individuals with poorer zileuton responses. From a mechanistic perspective, the ‘poor’ responders would be expected to have increased activation of PIK3CA and its targets, in parallel with increased LTB4 production that is resistant to zileuton inhibition. Higher concentrations of LTB4 would be available to bind to cell surface, G-protein coupled receptors, leading to increased PIK3CA signaling and resulting in increased proliferation of pro-inflammatory responses mediated by leukotriene activity. In contrast, our data indicate that ‘good’ responders would appropriately limit LTB4 production following zileuton inhibition, showing correspondingly diminished activity of this pathway due to reduced LTB4 concentrations.

This study is the first to report a role for zileuton in modulating the PI3K pathway in adults with asthma, for which zileuton is clinically used. However, a recent study showed that zileuton could selectively up-regulate phosphorylated Akt expression and down-regulate the over-expression of NF-KB p65 and COX-2, in a rat model of ischemic brain damage [32]. This finding is consistent with our observation that zileuton boosted AKT and phospho-AKT expression in samples from ‘poor’ responders (Fig. 3) and provides additional evidence for interaction of zileuton and the PI3K pathway.

While our results pointed toward a role for the PI3K pathway, we also identified four eQTLs that significantly reduced the expression of SELL, MTF2, and GAL, none of which have been annotated to the canonical PI3K pathway. While these signals were driven by a small number of subjects with homozygous variant genotypes, removal of these subjects did not alter the nominal significance of the four eQTLs. There is evidence that these genes may interact with the PI3K pathway, as well as with the CREB signaling pathway. For example, SELL is part of a co-expression network that includes PI3KCA, PIK3CG and PIK3R1 [33], and its protein product physically interacts with PRKCQ, which in turn interacts with PIK3CA through targeting CD28 and RASGRP3 [34]. MTF2 and CREB1 are also predicted to physically interact [35]; in addition, MTF2 interacts with PIK3R3 [36] which directly targets PIK3CA [36, 37]. Finally, as both GAL and PIK3CA are both activated by binding guanosine triphosphate (GTP), they are co-regulated through the G protein coupled receptor signaling cascade and may also have synergistic molecular roles. For example, GAL receptors can activate the phosphatidylinositol (PtdIns) signaling pathway, and there is evidence that PC12 apoptotic cell death mediated through Gal2R may be dependent on PI3K [38]. Further work is needed to ascertain the roles of these genes in P13K signaling in the context of leukotriene modifier response in asthma.

In summary, we used a multi-tiered, integrative pathway analysis of GWAS, eQTL, pQTL, and clinical data to prioritize pathways of interest, and functionally validated a candidate pathway with implications in asthma biology. This approach provides a means to identify novel biology for complex diseases and drug response phenotypes.

Supplemental Tables

Supplemental Figure Legends

Supplemental Figure 1

Supplemental Figure 2

Supplemental Figure 3