Comprehensive Pharmacogenomic Study Reveals an Important Role of UGT1A3 in Montelukast Pharmacokinetics.

To identify the genetic basis of interindividual variability in montelukast exposure, we determined its pharmacokinetics and sequenced 379 pharmacokinetic genes in 191 healthy volunteers. An intronic single nucleotide variation (SNV), strongly linked with UGT1A3*2, associated with reduced area under the plasma concentration-time curve (AUC0-∞ ) of montelukast (by 18% per copy of the minor allele; P = 1.83 × 10-10 ). UGT1A3*2 was associated with increased AUC0-∞ of montelukast acyl-glucuronide M1 and decreased AUC0-∞ of hydroxymetabolites M5R, M5S, and M6 (P < 10-9 ). Furthermore, SNVs in SLCO1B1 and ABCC9 were associated with the AUC0-∞ of M1 and M5R, respectively. In addition, a candidate gene analysis suggested that CYP2C8 and ABCC9 SNVs also affect the AUC0-∞ of montelukast. The found UGT1A3 and ABCC9 variants associated with increased expression of the respective genes in human liver samples. Montelukast and its hydroxymetabolites were glucuronidated by UGT1A3 in vitro. These results indicate that UGT1A3 plays an important role in montelukast pharmacokinetics, especially in UGT1A3*2 carriers.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ Interindividual variability in the pharmacokinetics of montelukast is high. It is not known how genetic variants in genes encoding drug‐metabolizing enzymes, membrane transporters, and regulatory proteins contribute to this variability.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ This study investigated whether genetic variants in pharmacokinetic genes affect the pharmacokinetics of montelukast and its metabolites.

WHAT THIS STUDY ADDS TO OUR KNOWLEDGE

☑ The results indicate that UGT1A3 plays an important role in montelukast pharmacokinetics, especially in carriers of the UGT1A3*2 allele associated with increased UGT1A3 expression. Furthermore, the results suggest that also CYP2C8 and ABCC9 variants affect the exposure to montelukast.

HOW THIS MIGHT CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE

☑ Genetic variants explain a significant proportion of interindividual variability in montelukast pharmacokinetics. This knowledge may aid in individualizing treatment with leukotriene receptor antagonists.

Montelukast is a leukotriene receptor antagonist, which is widely used in the treatment of asthma.1 After oral administration, montelukast is extensively metabolized and the majority of the metabolites are excreted into the bile.2 Previous studies have shown that the main enzyme involved in the oxidative metabolism of montelukast is cytochrome P450 (CYP) 2C8.3, 4 Also, CYP2C9, CYP3A4, and uridine diphosphate‐glucuronosyltransferase (UGT) 1A3 seem to contribute to the formation of montelukast metabolites.4, 5, 6, 7 Additionally, montelukast has been suggested to be a substrate of organic anion transporting polypeptide (OATP) 1B1, 1B3, and 2B1 transporters.8, 9, 10

High interindividual variability exists in the pharmacokinetics of montelukast. We hypothesized that variation in genes encoding drug‐metabolizing enzymes and membrane transporters, as well as proteins that affect their expression or biochemistry, contributes to this variability. Therefore, the aim of this study was to investigate the possible effects of genetic variability in these pharmacokinetic genes on montelukast pharmacokinetics. To this end, we determined the pharmacokinetics of montelukast after a 10‐mg dose in 191 healthy volunteers and fully sequenced 379 pharmacokinetic genes using massively parallel sequencing.

RESULTS

Montelukast pharmacogenomics

In the present study, substantial interindividual variability was observed in the pharmacokinetic variables of montelukast and its acyl‐glucuronide (M1) and hydroxymetabolites (M5R, M5S, and M6) (Supplementary Table S1). The areas under the plasma concentration–time curve from 0 h to infinity (AUC0‐∞) of montelukast, M1, M5R, M5S, and M6 varied 8.7‐fold, 13‐fold, 30‐fold, 22‐fold, and 23‐fold between individual subjects, respectively.

A total of 105,145 single nucleotide variations (SNVs) were found in the 379 analyzed pharmacokinetic genes (Supplementary Table S2), of which 46,064 had a minor allele frequency (MAF) of at least 0.05. In a stepwise linear regression analysis fixed for demographic covariates, nine common variants (MAF ≥0.05) in three genes were independently associated with montelukast or its metabolite pharmacokinetics at a Bonferroni‐corrected significance level of 1.09 × 10−6 (Table 1). The AUC0‐∞ of montelukast and its metabolites showed the strongest associations with variants in the UGT1A gene (Figure 1, Table 1). For montelukast AUC0‐∞, the strongest association was observed with rs7604115, located in the first intron of UGT1A3. The AUC0‐∞ of montelukast was 18% smaller per copy of the variant allele (P = 1.83 × 10−10). After adjusting for this variant, no other variant remained statistically significantly associated with montelukast AUC0‐∞. The investigated genetic variants had no significant effect on the peak plasma concentration (Cmax) or the elimination half‐life (t½) of montelukast.

Table 1

Results of the stepwise forward linear regression analysis of the effects of 46,064 SNVs in 379 genes on montelukast and its metabolite pharmacokinetics

Pharmacokinetic variable	dbSNP ID	Gene	Location	Nucleotide change	MAF	Effecta		P value
						Average	90%CI
Montelukast
AUC0‐∞	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	‐17.7%	‐21.6%, ‐13.7%	1.83 × 10−10
Cmax	—
t1/2	—
M1
AUC0‐∞	1. rs73063122	SLCO1B1	intron 11/14	c.1497 + 2246A>C	0.32	26.4%	18.5%, 34.8%	9.20 × 10−9
	2. rs3806592	UGT1A4	upstream	c.‐1531C>T	0.40	24.7%	17.4%, 32.3%	6.02 × 10−9
Cmax	1. rs3806592	UGT1A4	upstream	c.‐1531C>T	0.40	30.4%	20.8%, 40.9%	4.45 × 10−8
	2. rs4149056	SLCO1B1	exon 6/15	c.521T>C	0.22	34.7%	22.2%, 48.4%	9.90 × 10−7
t1/2	—
M1/montelukast AUC0‐∞ ratio	rs3806592	UGT1A4	upstream	c.‐1531C>T	0.40	52.9%	43.5%, 63.0%	4.14 × 10−22
M5R
AUC0‐∞	1. rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−45.7%	−49.5%, −41.7%	1.26 × 10−30
	2. rs704212	ABCC9	intron 12/37	c.1802 + 2622G>A	0.14	−27.6%	−34.6%, −19.8%	4.83 × 10−7
Cmax	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−33.9%	−38.6%, −28.8%	7.13 × 10−17
t1/2	rs7556676	UGT1A3	intron 1/4	c.868‐17430A>G	0.46	−24.8%	−28.2%, −21.2%	2.81 × 10−19
M5R/montelukast AUC0‐∞ ratio	rs4663969	UGT1A3	intron 1/4	c.867 + 16674C>A	0.46	−33.9%	−37.3%, −30.2%	7.06 × 10−27
M5S
AUC0‐∞	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−33.2%	−38.5%, −27.5%	8.14 × 10−14
Cmax	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−26.2%	−31.2%, −20.8%	1.82 × 10−11
t1/2	rs1875263	UGT1A4	upstream	c.‐1845C>T	0.40	−22.1%	−27.4%, −16.4%	2.17 × 10−8
M5S/montelukast AUC0‐∞ ratio	rs2361501	UGT1A3	intron 1/4	c.867 + 51A>T	0.47	−20.1%	−24.7%, −15.3%	1.71 × 10−9
M6
AUC0‐∞	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−39.4%	−43.2%, −35.4%	2.89 × 10−27
Cmax	rs7604115	UGT1A3	intron 1/4	c.868‐17564C>T	0.40	−34.0%	−37.8%, −29.9%	2.22 × 10−23
t1/2	—
M6/montelukast AUC0‐∞ ratio	rs4663969	UGT1A3	intron 1/4	c.867 + 16674C>A	0.46	−25.4%	−28.1%, −22.7%	2.77 × 10−29

AUC0‐∞, area under the plasma concentration‐time curve from 0 h to infinity; CI, confidence interval; Cmax, peak plasma concentration; dbSNP, National Center for Biotechnology Information Short Genetic Variations database; MAF, minor allele frequency; SNV, single nucleotide variation; t½, elimination half‐life.

Per copy of the minor allele.

Figure 1

The associations of SNVs in 379 pharmacokinetic genes (a) and in the UGT1A gene (b) with montelukast AUC0‐∞, and LD plot of the UGT1A gene (c). The Y‐axes in (a,b) describe the negative logarithm of the P value for each SNV and the horizontal lines indicate the Bonferroni‐corrected significance level of 1.09 × 10−6. The X‐axis in (a) shows individual SNVs grouped by protein function. The locations of UGT1A exons and linkage disequilibrium blocks are depicted on the top of (b,c).

In agreement with the effects of UGT1A variants on parent montelukast, the AUC0‐∞ of montelukast acyl‐glucuronide (M1) was 25% larger per copy of the UGT1A rs3806592 variant allele (P = 6.02 × 10−9) (Table 1). Rs3806592 is in a strong linkage disequilibrium with the UGT1A3 rs7604115 SNV associated with the AUC0‐∞ of parent montelukast (r  = 0.95, P = 3.65 × 10−41). UGT1A variants were also significantly associated with the Cmax of M1 and the M1/montelukast AUC0‐∞ ratio. Furthermore, the solute carrier organic anion transporter gene 1B1 (SLCO1B1) SNVs rs73063122 and rs4149056 were significantly associated with the AUC0‐∞ and Cmax of M1, respectively. These two SNVs are in a strong linkage disequilibrium with each other (r  = 0.60, P = 5.73 × 10−27), suggesting that both of these associations are due to the rs4149056 missense SNV known to markedly impair the activity of OATP1B1.11

The effects of UGT1A variants on the AUC0‐∞ of the hydroxylated M5R, M5S, and M6 metabolites of montelukast were larger than what was observed for montelukast (Table 1). The AUC0‐∞ of M5R, M5S, and M6 were 46% (P = 1.26 × 10−30), 33% (P = 8.14 × 10−14), and 39% (P = 2.89 × 10−27) smaller per copy of the rs7604115 variant allele, respectively. UGT1A variants were also significantly associated with the Cmax and the metabolite/montelukast AUC0‐∞ ratios of M5R, M5S, and M6 and the t½ of M5R and M5S. The AUC0‐∞ of M5R was also associated with the rs704212 SNV in the ATP binding cassette subfamily C member 9 (ABCC9) transporter gene.

UGT1A linkage disequilibrium and haplotype analyses

The UGT1A gene region (±20 kb) was found to consist of nine linkage disequilibrium (LD) blocks (Figure 1). The SNVs associated with montelukast or its metabolite pharmacokinetics were located in blocks 3 and 4. Within these blocks, 28 haplotypes were inferred (Supplementary Figure S1). The SNVs showing the strongest associations with montelukast pharmacokinetics were strongly linked to the missense variants rs3821242 and rs6431625 (r ≥ 0.69), which together define the UGT1A3*2 haplotype. Based on the missense variants rs3821242, rs6431625, and rs45449995, the inferred haplotypes were grouped to subtypes of UGT1A3*1 (wildtype), *2 (rs3821242 and rs6431625; MAF 0.39), *3 (rs3821242; MAF 0.060), and *6 (rs3821242, rs6431625, and rs45449995; MAF 0.018). The effects of the UGT1A3*2 haplotype on montelukast and its metabolite pharmacokinetics were similar to the effects of the individual, intronic UGT1A SNVs (Tables 1 and 2, Figure 2).

Table 2

Results of the stepwise forward linear regression analysis of the effects of UGT1A3*2, SLCO1B1 rs4149056, and ABCC9 rs704212 on montelukast and its metabolite pharmacokinetics

Pharmacokinetic variable	Haplotype/dpSNP ID	Effecta
		Average	90% CI	P value
Montelukast
AUC0‐∞	UGT1A3*2	−17.5%	−21.4%, −13.5%	2.75 × 10−10
Cmax	—
t1/2	—
M1
AUC0‐∞	UGT1A3*2	26.4%	19.2%, 34.1%	5.13 × 10−10
	SLCO1B1 rs4149056	34.4%	24.8%, 44.8%	5.20 × 10−10
Cmax	UGT1A3*2	30.0%	20.4%, 40.5%	6.71 × 10−8
	SLCO1B1 rs4149056	35.5%	22.9%, 49.4%	6.54 × 10−7
t1/2	—
M1/montelukast AUC0‐∞ ratio	UGT1A3*2	52.8%	43.3%, 62.8%	5.82 × 10−22
M5R
AUC0‐∞	UGT1A3*2	−45.4%	−49.2%, −41.2%	1.20 × 10−29
	ABCC9 rs704212	−27.6%	−34.8%, −19.7%	6.90 × 10−7
Cmax	UGT1A3*2	−33.8%	−38.5%, −28.7%	7.49 × 10−17
t1/2	UGT1A3*2	−25.3%	−28.9%, −21.6%	1.26 × 10−18
M5R/montelukast AUC0‐∞ ratio	UGT1A3*2	−35.0%	−38.6%, −31.2%	1.01 × 10−26
M5S
AUC0‐∞	UGT1A3*2	−33.1%	−38.4%, −27.4%	9.01 × 10−14
Cmax	UGT1A3*2	−26.0%	−31.0%, −20.6%	2.65 × 10−11
t1/2	UGT1A3*2	−21.8%	−27.1%, −16.1%	2.79 × 10−8
M5S/montelukast AUC0‐∞ ratio	UGT1A3*2	−18.9%	−23.8%, −13.8%	5.87 × 10−8
M6
AUC0‐∞	UGT1A3*2	−39.3%	−43.1%, −35.2%	3.41 × 10−27
Cmax	UGT1A3*2	−33.9%	−37.7%, −29.8%	2.65 × 10−23
t1/2	—
M6/montelukast AUC0‐∞ ratio	UGT1A3*2	−26.1%	−28.9%, −23.1%	3.56 × 10−27

AUC0‐∞, area under the plasma concentration‐time curve from 0 h to infinity; CI, confidence interval; Cmax, peak plasma concentration; dbSNP, National Center for Biotechnology Information Short Genetic Variations database; t½, elimination half‐life.

Per copy of the minor allele.

Figure 2

Geometric mean (90% CI) plasma concentrations of montelukast, its acyl‐glucuronide (M1), and hydroxymetabolites (M5R, M5S, M6) after a single 10‐mg oral dose of montelukast in 191 healthy volunteers with different UGT1A3 genotypes. Open circles indicate noncarriers of UGT1A3*2 (n = 72), solid circles subjects heterozygous for the UGT1A3*2 (n = 90), and solid triangles subjects homozygous for the UGT1A3*2 haplotype (n = 29). The insets depict the same data on a semilogarithmic scale. Plasma concentrations of montelukast were adjusted for BSA and those of M1 and M5S for lean body weight, as according to the linear regression models of their AUC0‐∞ values.

Functional validation

Next, we investigated whether montelukast or its hydroxymetabolites are substrates of UGT1A enzymes in vitro. As the most strongly associated variants are localized around the first exon of UGT1A3, we focused on this enzyme, together with UGT1A1 and UGT1A9, which are also known to catalyze the glucuronidation of carboxylic acids.12 Montelukast, M5R, M5S, and M6 were all metabolized by UGT1A3, but not significantly by UGT1A1 or UGT1A9 (Supplementary Figure S2).

To further elucidate the mechanisms of associations between UGT1A, ABCC9, and SLCO1B1 SNVs and montelukast or its metabolite pharmacokinetics, we then investigated the effects of the SNVs on the respective gene expression in human liver samples (Table 3). UGT1A3 expression showed a strong association with the SNVs associated with montelukast pharmacokinetics. The strongest association was observed with rs4663969, present in both UGT1A3*2 and *3 haplotypes. UGT1A3 gene expression was 24% higher per copy of the variant allele (P = 2.27 × 10−4). Of the UGT1A3*2, *3, and *6 haplotypes, only UGT1A3*2 was significantly associated with increased UGT1A3 gene expression. The UGT1A3 expression was 24% higher per copy of the UGT1A3*2 haplotype (P = 2.08 × 10−4) (Figure 3). ABCC9 expression was 25% higher (P = 3.41 × 10−4) per copy of the rs704212 variant allele. SLCO1B1 expression was 13% lower (P = 0.0146) per copy of the rs4149056 variant allele. Of SLCO1B1 haplotypes, *15 (rs2306283 and rs4149056) was associated with a 15% lower SLCO1B1 expression per copy of the haplotype (P = 0.00733).

Table 3

Results of the linear regression analysis of the effects of UGT1A SNVs and haplotypes on UGT1A3 mRNA expression, ABCC9 rs704212 SNV on ABCC9 mRNA expression, and SLCO1B1 SNVs and haplotypes on SLCO1B1 mRNA expression in liver samples from patients undergoing laparoscopic gastric bypass operation

Gene	n	dpSNP ID/Haplotype	MAF	Effecta		P value
			(n = 185‐197)	Average	90% CI
UGT1A3	136	rs7604115	0.44	22.3%	11.3%, 34.5%	5.65 × 10−4
	136	rs3806592	0.43	22.9%	11.7%, 35.2%	4.78 × 10−4
	136	rs7556676	0.47	23.2%	12.4%, 35.0%	2.57 × 10−4
	136	rs4663969	0.47	23.6%	12.7%, 35.6%	2.27 × 10−4
	136	rs1875263	0.44	22.9%	11.7%, 35.2%	4.78 × 10−4
	136	UGT1A3*2	0.41	23.8%	12.8%, 35.8%	2.08 × 10−4
	136	UGT1A3*3	0.058	−1.2%	−18.7%, 20.0%	0.918
	136	UGT1A3*6	0.010	29.9%	−13.3%, 94.7%	0.286
ABCC9	187	rs704212	0.11	25.2%	13.1%, 38.6%	3.41 × 10−4
SLCO1B1	124	rs4149056	0.24	−13.5%	−21.4%, −4.7%	0.0146
	124	rs2306283	0.47	−7.5%	−14.7%, 0.4%	0.117
	123	SLCO1B1*1B	0.18	10.4%	0.8%, 20.9%	0.0740
	123	SLCO1B1*5	0.019	11.2%	−18.3%, 51.3%	0.570
	123	SLCO1B1*15	0.22	−15.2%	−23.2%, −6.3%	0.00733

MAF, minor allele frequency.

Per copy of the minor allele.

Figure 3

Boxplots of the effects of UGT1A3*2, *3, and *6 haplotypes on UGT1A3 mRNA expression in human liver samples. The horizontal lines inside the boxes represent the median, the box edges show the lower and upper quartiles, and the whiskers show the 10th and 90th percentiles. Individual data points are given as circles for men and as triangles for women.

Candidate gene analysis

We next carried out a candidate gene analysis for montelukast AUC0‐∞, focusing on common (MAF ≥0.05) missense variants in genes suggested to be involved in montelukast pharmacokinetics (CYP2C8, CYP2C9, CYP3A4, SLCO1B1, SLCO1B3, SLCO2B1, UGT1A3), as well as the ABCC9 rs704212 SNV (Supplementary Table S3). In a stepwise linear regression analysis, UGT1A3*2 was associated with a 17% (P = 2.99 × 10−10), ABCC9 rs704212 with a 14% (P = 2.19 × 10−4), and CYP2C8*3 (rs10509681 and rs11572080) with an 11% (P = 0.00659) reduced, and CYP2C8*4 (rs1058930) with a 13% (P = 0.0184) increased AUC0‐∞ of montelukast per copy of each minor allele (Table 4; adjusted R  = 0.41).

Table 4

Results of the candidate gene analysis on montelukast AUC0‐∞

Variable	MAF	Effecta		P value	Adjusted R2 for each step
		Average	90% CI
BSA	—	−11.5%	−14.0%, −9.0%	2.10 × 10−11	0.17
UGT1A3*2	0.39	−16.8%	−20.5%, −12.9%	2.99 × 10−10	0.32
ABCC9 rs704212	0.14	−13.7%	−19.1%, −7.9%	2.19 × 10−4	0.37
CYP2C8*3	0.11	−11.3%	−17.5%, −4.7%	0.00659	0.40
CYP2C8*4	0.07	12.7%	3.7%, 22.5%	0.0184	0.41

BSA, body surface area; AUC0‐∞, area under the plasma concentration‐time curve from 0 h to infinity; CI, confidence interval; MAF, minor allele frequency.

BSA effect per 10% increase; genetic variant effect per copy of the minor allele.

DISCUSSION

In this study we used targeted massively parallel sequencing of 379 pharmacokinetic genes to characterize the genetic basis of interindividual variability in montelukast pharmacokinetics. The UGT1A3*2 haplotype and variants located around the first exon of UGT1A3 were strongly associated with the systemic exposure to parent montelukast and its metabolites. We further demonstrated that montelukast and its hydroxymetabolites are substrates of UGT1A3 in vitro and that the variants associated with the pharmacokinetic variables also significantly affect UGT1A3 gene expression in human liver. In addition, our results indicate involvement of ABCC9 and SLCO1B1 in the pharmacokinetics of montelukast metabolites. Moreover, a candidate gene approach suggested that, in addition to UGT1A3 variants, also CYP2C8 and ABCC9 variants affect parent montelukast exposure.

The UGT1A gene encodes the UGT1A family enzymes, which catalyze the formation of hydrophilic glucuronide metabolites.13 Individual UGT1A genes have unique first exons but share exons 2–5. Strong linkage disequilibrium exists throughout the whole UGT1A gene (Figure 1 c). The intronic UGT1A SNVs showing the strongest associations with montelukast or its metabolite pharmacokinetics are strongly linked to the UGT1A3 missense variants rs3821242 (c.31T>C, p.Trp11Arg) and rs6431625 (c.140T>C, p.Val47Ala) that together define the haplotype UGT1A3*2. The UGT1A3*2 haplotype and the individual intronic SNVs similarly reduced the exposure to montelukast and its hydroxymetabolites, and increased the exposure to montelukast acyl‐glucuronide, indicating enhanced glucuronidation. Consistently, UGT1A3*2 has previously been shown to increase the metabolism of the UGT1A3 substrates atorvastatin, telmisartan, and febuxostat in humans.14, 15, 16

The UGT1A3*2 haplotype and the UGT1A SNVs associated with montelukast pharmacokinetics significantly increased UGT1A3 mRNA expression in human liver samples. Similarly, a previous study also showed that UGT1A3 mRNA and protein expression are significantly increased in UGT1A3*2 carriers.14 In addition to UGT1A3*2, the UGT1A3*6 haplotype has been associated with increased UGT1A3 expression.14 In our study, the number of UGT1A3*6 carriers was relatively small, and only a tendency towards increased UGT1A3 expression could be observed (Figure 3).

The causal variant of UGT1A3*2 affecting UGT1A3 expression has remained unknown.14 We identified several intronic SNVs that are strongly linked to UGT1A3*2 (Supplementary Figure S1). Of these, rs3806597 (c.‐204A>G) is located on a proposed farnesoid X receptor (FXR) binding site upstream of UGT1A3, but the variant allele has not affected FXR‐mediated induction of UGT1A3 by the bile acid chenodeoxycholic acid.17 None of the other strongly linked variants appear to be located in a transcription factor binding site upstream of UGT1A3.18 In addition to being associated with expression, the UGT1A3*2 missense variants might also alter the enzymatic activity of UGT1A3. However, in vitro studies with these variants have shown conflicting results, with reduced, increased, and unchanged activity.19, 20, 21 Altogether, although the causal variant cannot be identified, the increased UGT1A3 expression in association with UGT1A3*2 provides a mechanistic explanation for our pharmacokinetic results.

In accordance with a recently published in vitro study,7 our results demonstrate that montelukast is efficiently glucuronidated by UGT1A3. We also showed that montelukast hydroxymetabolites M5R, M5S, and M6 are efficiently glucuronidated by UGT1A3 in vitro, which could explain why the concentrations of the hydroxymetabolites were affected more by the UGT1A3 variants than those of parent montelukast. Taken together, our results indicate that glucuronidation plays an important role in the metabolism of both parent montelukast and its hydroxymetabolites. However, due to low plasma concentrations, the hydroxymetabolites are unlikely to contribute to the pharmacological effects of montelukast.

In addition to UGT1A3 variants, SNVs in the SLCO1B1 gene were associated with significantly higher AUC and Cmax of montelukast acyl‐glucuronide, M1. SLCO1B1 encodes OATP1B1, an influx transporter mediating the hepatic uptake of its substrates from sinusoidal blood.11 Of the associated SNVs, the rs4149056 missense variant (c.521T>C, p.Val174Ala, SLCO1B1*5 or *15) markedly impairs the function of OATP1B1, as demonstrated both in vitro and in vivo in humans.11, 22, 23, 24 Many glucuronide conjugates are OATP1B1 substrates and our results suggest that also M1 is a substrate of OATP1B1. In a previous study, the rs2306283 (c.388A>G, p.Asn130Asp) missense variant was associated with increased SLCO1B1 expression in human liver samples (n = 143), whereas rs4149056 had no effect.25 In our study, however, rs4149056 was associated with decreased SLCO1B1 expression, but rs2306283 appeared to have no effect. Nevertheless, the SLCO1B1*1B haplotype, which contains rs2306283 without rs4149056, showed a tendency for a 10% increased expression of SLCO1B1 per copy of the haplotype.

A recent study suggested that montelukast is a substrate of OATP1B1, 1B3, and 2B1, and that OATP‐CYP2C8 interplay is an important determinant of parent montelukast pharmacokinetics.10 Our finding that the OATP1B1 function‐impairing rs4149056 SNV did not significantly affect the pharmacokinetics of parent montelukast, even in the candidate gene analysis with no correction for multiple testing, indicates that OATP1B1 is not of major importance for the hepatic elimination of montelukast in vivo in humans. Moreover, we fully sequenced the SLCO1B3 and SLCO2B1 genes and found no association between montelukast AUC and variants in these genes, indicating that these genes are also not important determinants of interindividual variability in montelukast pharmacokinetics.

Interestingly, the ABCC9 rs704212 SNV was associated with reduced plasma concentrations of montelukast hydroxymetabolite M5R, and with increased mRNA expression of ABCC9 in human liver samples. Based on these results, the SNV was selected for the candidate gene analysis, where it was found to be associated also with reduced AUC of montelukast. The ABCC9 gene encodes the sulfonylurea receptor 2 (SUR2), which is a subunit of an adenosine triphosphate‐sensitive K+ channel and has no identified transporter function.26 Notably, another ABCC9 SNV, rs1283807, which is in a strong linkage disequilibrium with rs704212 (r 2 = 0.82, P = 1.97 × 10−37), has previously been associated with impaired efficacy of angiotensin II receptor antagonists in a genome‐wide association study.27 The same SNV has also been associated with the risk of virologic failure during efavirenz‐containing antiretroviral therapy via an epistatic interaction with an SNV in transporter 1, ATP‐binding cassette subfamily B member (TAP1).28 Our study is the first to show an association between ABCC9 and drug pharmacokinetics, suggesting a pharmacokinetic explanation also for the findings concerning angiotensin II receptor antagonists and efavirenz.

It has been estimated that 80% of montelukast metabolism is mediated by CYP2C8.6 In the analysis of all the 379 genes, CYP2C8 variants were not significantly associated with montelukast pharmacokinetics. However, in the candidate gene analysis, the AUC of montelukast was reduced by the CYP2C8*3 allele and increased by the CYP2C8*4 allele. These findings are consistent with previous studies indicating increased CYP2C8 activity with CYP2C8*3 and reduced activity with CYP2C8*4 in vivo in humans.29, 30, 31, 32, 33, 34, 35 In the present study, UGT1A3*2 reduced the AUC of montelukast by 18% per allele copy, equivalent to a 22% increase in oral clearance. The simultaneous increase in the M1/montelukast AUC ratio suggests an about 50% increase in UGT1A3‐mediated clearance per UGT1A3*2 allele. Thus, it seems plausible that the contribution of UGT1A3 to the metabolism of montelukast varies between 35–60%, depending on the UGT1A3 genotype. In previous studies, the CYP2C8 inhibitor gemfibrozil has markedly increased the AUC of montelukast.3, 6 Interestingly, gemfibrozil inhibits UGT1A3 in vitro.7 Therefore, it is possible that UGT1A3 also contributes to the gemfibrozil‐montelukast interaction. In any case, the UGT1A3 genotype might affect the magnitude of this interaction.

Considerable interindividual variability exists in the efficacy of montelukast, but no firm clinical or genetic predictors of montelukast response have yet been identified.36, 37 In clinical practice, montelukast is administered on a regular basis. The steady‐state plasma concentrations for montelukast are predictable from its pharmacokinetic parameters measured after a single dose.38 Therefore, the effects of genetic variants on the steady‐state plasma concentrations of montelukast should be similar to the effects on the AUC0‐∞ of montelukast observed in our study after a single dose. Because montelukast dose‐dependently improves chronic asthma,39 the reduced plasma concentrations of montelukast due to the UGT1A3*2 haplotype might impair its efficacy. Even though the effect of the UGT1A3*2 haplotype on montelukast pharmacokinetics is modest, the haplotype is common (Supplementary Figure S3). Therefore, it may be an important factor explaining the variability in montelukast response at the population level. Body surface area and the UGT1A3*2 allele together explained 32% of interindividual variability in montelukast exposure. When the ABCC9 and CYP2C8 alleles were added to the model, this percentage increased to 41%. Together with pharmacodynamic markers,37 this knowledge might aid in individualizing treatment with leukotriene receptor antagonists.

In conclusion, genetic variability in UGT1A3 significantly affects montelukast pharmacokinetics. This indicates that glucuronidation via UGT1A3 has a larger role in the metabolism of montelukast than previously thought, especially in subjects carrying the UGT1A3*2 haplotype. These results also further confirm that the UGT1A3*2 haplotype enhances the glucuronidation of UGT1A3 substrates. Moreover, the candidate gene analysis suggested that also CYP2C8*3, CYP2C8*4, and ABCC9 rs704212 affect the pharmacokinetics of montelukast.

METHODS

Study participants

In all, 201 healthy unrelated Finnish Caucasian volunteers participated in the pharmacokinetic study after giving written informed consent. Their health was confirmed by medical history, clinical examination, and laboratory tests. Participants were not on any continuous medication nor were tobacco smokers. The study was approved by the Coordinating Ethics Committee of the Hospital District of Helsinki and Uusimaa and the Finnish Medicines Agency Fimea. Ten participants discontinued the study before montelukast administration and thus the pharmacokinetic data were obtained from 191 participants. Of these, 92 were women and 99 men. Their mean ± SD age was 24 ± 4 years, height 174 ± 9 cm, weight 70 ± 12 kg, and body mass index (BMI) 22.8 ± 2.5 kg/m2.

A whole‐blood DNA sample and a liver biopsy were obtained from 201 patients undergoing laparoscopic gastric bypass operation at the Kuopio University Hospital, as part of the Kuopio Obesity Surgery Study.40, 41, 42 Good quality RNA expression data and genotypes were obtained from 188 patients. These patients consisted of 129 women and 59 men (mean ± SD: age 48 ± 9 years and BMI 44 ± 6 kg/m2). Seventy‐six patients had type 2 diabetes, 34 had nonalcoholic fatty liver, 32 had nonalcoholic steatohepatitis, and 58 used lipid‐lowering medication. The degree of liver steatosis was graded from 0 to 3 and that of lobular inflammation from 0 to 2. Written informed consent was obtained from all patients and the study protocol was approved by the Ethics Committee of the Northern Savo Hospital District.

Montelukast pharmacokinetics

After fasting overnight, the healthy volunteers ingested a single 10‐mg dose of montelukast (Singulair tablet, Merck Sharp & Dohme, Haarlem, The Netherlands) with 150 ml of water at 8 am. Standardized meals were served at 4, 7, and 10 h after montelukast ingestion. Timed blood samples (4–9 ml each) were collected into light‐protected ethylenediaminetetraacetic acid (EDTA) tubes prior to and up to 24 h after montelukast administration. Tubes were immediately placed on ice and plasma was separated within 30 min. Samples were stored at –70°C until analysis.

The concentrations of plasma montelukast, montelukast acyl‐glucuronide (M1), montelukast 1,2‐diol (M6), 21R‐hydroxy montelukast (M5R), and 21S‐hydroxy montelukast (M5S) were measured using a Nexera X2 liquid chromatography instrument (Shimadzu, Kyoto, Japan) interfaced with a 5500 Qtrap tandem mass spectrometer (AB Sciex, Toronto, ON). Prior to quantification, the plasma sample was purified from proteins and phospholipids using a Phree phospholipid removal plate (Phenomenex, Torrance, CA) according to the manufacturer's instructions. In short, the plasma sample was mixed with acetonitrile containing 1% formic acid and internal standards (1:4 v/v) and drawn through the Phree cartridges. The chromatographic separation was achieved on a reversed‐phase Kinetex C18 analytical column (100 × 2.1 mm internal diameter, 2.6 μm particle size; Phenomenex) using 2 mM ammonium acetate (A) (pH 4.0) adjusted with 98% formic acid and acetonitrile (B) as a mobile phase. The injection volume was 3 μl and the column temperature was held at 30°C. The gradient profile was set as follows: a linear increase from 35% B to 62% B over 2.2 min, held at 62% B for 2 min, a gradient from 62% B to 95% B over 0.6 min, and maintained at 95% B for 1.4 min followed by a reversion to the initial conditions. The mass spectrometer was operated in a positive electrospray ionization mode (ESI+) employing scheduled multiple reaction monitoring (MRM) for the optimal dwell time for each analyte. A corresponding deuterated reference compound served as internal standard for each analyte, except for M6, which utilized montelukast‐d6. The target mass‐to‐charge ratios (m/z) for montelukast, M1, M5R/M5S, and M6 were 586→422, 762→422, 602→147, and 602→438, and the limits of quantification (ng/ml) were 1.0, 0.2, 0.1, and 0.1, respectively. The day‐to‐day coefficient of variation (CV) was below 10% at relevant concentrations for all analytes. The AUC0‐∞, Cmax, and t½ values were calculated for montelukast, M1, M5R, M5S, and M6 with standard noncompartmental methods using Phoenix WinNonlin, v. 6.3 (Certara, Princeton, NJ).

DNA sequencing and genotyping

Genomic DNA was extracted using the Maxwell 16 LEV Blood DNA Kit on a Maxwell 16 Research automated nucleic acid extraction system (Promega, Madison, WI; pharmacokinetic study) or DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany; liver samples). DNA concentration and absorbance 260/280 ratio (A260/A280) were determined with the NanoDrop 2000 UV‐Vis Spectrophotometer (Thermo Fisher Scientific, Waltham, MA).

Targeted massively parallel sequencing of the 379 pharmacokinetic genes ±20 kb (Supplementary Table S2; genome build GRCh37) was performed in all pharmacokinetic study participants (n = 201). For library preparation, 3 μg of genomic DNA was processed according to the NEBNext DNA Sample Prep protocol (New England BioLabs, Ipswich, MA). Target enrichment capture was performed using the NimbleGen SeqCap EZ Choice capture protocol (Roche Sequencing, Pleasanton, CA). Sequencing was done on the Illumina HiSeq2000 platform with 100 bp paired‐end reads (Illumina, San Diego, CA). Quality control, short read alignment, and variant calling and annotation were carried out using an in‐house‐developed pipeline, as described previously.43 The sequencing and bioinformatics pipelines were carried out at the Technology Centre at the Institute for Molecular Medicine Finland (Helsinki, Finland). Mean coverage depth was 37.2×. Coverage depth ≥10×, Hardy–Weinberg equilibrium P < 3.15 × 10−7 (Bonferroni‐correction), and proportion missing ≤0.05 were employed as quality thresholds for including genotype data in statistical analysis.

The pharmacokinetic study participants and liver samples were genotyped for the UGT1A rs7604115, rs3806592, rs7556676, rs4663969, rs1875263, rs3821242, rs6431625, and rs45449995, ABCC9 rs704212, and SLCO1B1 rs4149056 and rs2306283 SNVs with TaqMan genotyping assays on a QuantStudio 12K Flex Real‐Time PCR System according to the manufacturer's protocol (Thermo Fisher Scientific). Call identity with sequencing data was 99–100%. In case of discordant results, genotypes obtained by sequencing were used in the statistical analysis.

Reverse transcription quantitative real‐time PCR

RNA from liver samples was extracted using the miRNeasy Mini Kit (Qiagen, Chatsworth, CA) and stored at –80°C. RNA was reverse‐transcribed using the SuperScript VILO cDNA Synthesis Kit, according to the manufacturer's instructions (Thermo Fisher Scientific). The cDNA samples were preamplified (14 cycles) with a custom TaqMan pre amp pool containing the assays for UGT1A, ABCC9, SLCO1B1 and reference genes before quantitative real‐time PCR (qPCR), according to the manufacturer's instructions (Thermo Fisher Scientific).

The qPCR was carried out using OpenArray technology on the QuantStudio 12K Flex Real‐Time PCR System. The custom OpenArray plate contained UGT1A3 (Hs04194492_g1), ABCC9 (Hs00245832_m1), SLCO1B1 (Hs00272374_m1), and reference gene assays, allowing the within‐sample normalization with multiple reference genes. The reference genes were comprised of actin beta (ACTB; Hs01060665_g1), ribosomal protein lateral stalk subunit P0 (RPLP0; Hs99999902_m1), glyceraldehyde‐3‐phosphate dehydrogenase (GAPDH; Hs02758991_g1), and beta‐2‐microglobulin (B2M; Hs00984230_m1). Prior to thermal cycling, a volume of 2.5 μL of TaqMan OpenArray Real‐Time PCR Master Mix (Thermo Fisher Scientific) was combined with 1.3 μL nuclease‐free water followed by 1.2 μL of cDNA, and the sample was loaded to OpenArray plate. The samples were processed on the Freedom EVO 150 automated liquid handling system (Tecan Group, Männedorf, Switzerland) and the OpenArray AccuFill instrument (Thermo Fisher Scientific). Data were analyzed with the ExpressionSuite software (v. 1.0.3, Thermo Fisher Scientific) using the comparative ΔCq method.

In vitro studies with UGT1A1, 1A3, and 1A9 recombinant enzymes

The incubation mixtures (triplicate samples) contained 100 mM phosphate buffer, pH 7.4, 5 mM MgCl2, 0.3 mg/ml of supersomal protein (Corning, Woburn, MA), 2 mM UDP‐glucuronic acid (UDPGA, triammonium salt), and 0.1 μM montelukast, M5R, M5S, or M6, in a total volume of 500 μl. Control incubations (duplicate samples) were performed with control supersomes, which do not contain active UGT enzyme. The reactions were initiated by the addition of the substrate and incubated in a shaking water bath at 37°C. Reactions were stopped by moving 50 μl samples to 50 μl acetonitrile containing 1% formic acid and internal standards at timepoints 0, 6, 12, 20, 40, and 60 min. After centrifugation, the supernatants were analyzed using SCIEX API 2000 tandem mass spectrometer (AB Sciex), as described previously.4 Depletion (% remaining) of the substrates at each timepoint was calculated for UGT1A1, 1A3, 1A9, and control samples.

Statistical analysis

The data were analyzed with the statistical programs JMP Genomics 7.0 (SAS Institute, Cary, NC) and IBM SPSS 22.0 for Windows (Armonk, NY). The pharmacokinetic variables were logarithmically transformed before analysis.44 Possible effects of demographic covariates on pharmacokinetic variables and gene expression were investigated using stepwise linear regression analysis, with P‐value thresholds of 0.05 for entry and 0.10 for removal. Sex, body weight, lean body weight,45 and body surface area, (BSA)46 were tested as demographic covariates for pharmacokinetic data, and age, sex, BMI, type 2 diabetes, degree of lobular inflammation and hepatic steatosis, and use of lipid‐lowering medication for gene expression data (Supplementary Table S4). Possible effects of genetic variants on pharmacokinetic variables were investigated using stepwise linear regression analysis fixed for significant demographic covariates. A Bonferroni‐corrected P‐value threshold of 1.09 × 10−6 was employed for the 379 gene and UGT1A3 haplotype‐based analysis and thresholds of 0.05 for entry and 0.10 for removal for the candidate gene analysis. The effects of genetic variants on gene expression were investigated with linear regression analysis fixed for significant demographic covariates, with P < 0.05 considered statistically significant. Additive coding was employed for genetic variants and multiallelic variants were expanded. The LD blocks for UGT1A (±20 kb) SNVs (MAF >0.05) were created using a previously described algorithm,47, 48 as implemented in JMP Genomics. Haplotype computations were performed with PHASE v. 2.1.1.49, 50

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

AUTHOR CONTRIBUTIONS

P.H. and M.Ni. wrote the article; P.H., A.T., J.T.B., and M.Ni. designed the research; P.H., A.T., M.Ne., T.T., M.P‐H., V.K., V.T.M., J.P., J.T.B., and M.Ni. performed the research; P.H. and M.Ni. analyzed the data.

Supporting Information 1

Click here for additional data file.

Supporting Information 2

Click here for additional data file.

Supporting Information 3

Click here for additional data file.

Supporting Information 4

Click here for additional data file.

Supporting Information 5

Click here for additional data file.

Supporting Information 6

Click here for additional data file.

Supporting Information 7

Click here for additional data file.

Supporting Information 8

Click here for additional data file.