Introduction

ATP-binding cassette (ABC) transporters have diverse and important functions in maintaining cellular homeostasis of endogenous solutes and various xenobiotics, including many clinically important drugs.1, 2 In human beings, 49 ABC proteins have been identified and classified into seven subfamilies (ABCA–ABCG). Most of these are membrane transporters that utilize the energy from ATP hydrolysis to mediate efflux of solutes. The different members exhibit widely varying patterns of tissue and subcellular localization, and have substrate specificities that range from endogenous solutes, such as cholesterol and bile salts (ABCA1, ABCB11, ABCG1, ABCG5 and ABCG8),3, 4, 5, 6 phospholipids (ABCB4)7 and long-chain fatty acids (ABCD1 and ABCD2),8, 9 to drugs and their metabolites (ABCB1, ABCG2 and several members of the ABCC family).2, 10, 11

The important physiological functions of ABC transporters are highlighted by the many genetic disorders resulting from mutations in these proteins, for example, Tangier's disease (ABCA1),4 progressive familial intrahepatic cholestasis (ABCB4 and ABCB11),12, 13 adrenoleukodystrophy (ABCD1),9 sitosterolemia (ABCG5 and ABCG8)3, 6 and cystic fibrosis (CFTR, also referred to as ABCC7).14 Although such severe clinical phenotypes are typically attributable to loss-of-function or reduced-function mutations, natural variation in transporter expression levels may also result in detrimental clinical conditions. For instance, low expression of ABCA1 is associated with reduced plasma levels of high-density lipoprotein, although less extreme than that observed in Tangier's disease patients.15, 16 Similarly, expression level variation in drug-transporting ABCs such as ABCB1, ABCC2 and ABCG2 can affect intestinal absorption, biliary and renal excretion, and tissue distribution of substrate drugs, which in turn can lead to the loss of pharmacological effects, or to drug-induced toxicities. Identification of regulatory elements in ABC transporter genes is therefore essential for understanding the factors that control the expression levels of these important transporters.

Although many cis-regulatory elements have been identified in the genome, discovery of such elements is difficult. The most common approach to the discovery of cis-regulatory regions has involved comparative genomic methods as a first step to identify conserved sequence motifs, followed by further computational analyses to discover transcription factor binding sites.17, 18, 19 Although these methods have resulted in the identification of many regulatory regions, use of evolutionary conservation as a first step may not be reliable. In particular, most regulatory motifs are small, typically consisting of only a few base pairs, and are therefore difficult to identify using sequence comparisons.

Recently, several groups have used publicly available mRNA expression levels in lymphoblastoid cell lines (LBLs) along with genetic variation on a genome-wide scale to identify cis-regulatory variation in the human genome.20, 21, 22, 23 These studies have revealed global patterns of expression quantitative trait loci (eQTLs), including population differences, and support an abundance of polymorphisms in and nearby genes, which associate with gene expression levels. In this study, we leverage these publicly available data sets, with the specific aim of determining population and individual expression level differences in ABC transporters and to discover and functionally annotate cis-regulatory regions that may harbor enhancer or silencer elements that regulate the expression of ABC transporters.

Materials and methods

ABC gene expression data

Gene expression data were downloaded from the Wellcome Trust Sanger Institute web site (http://www.sanger.ac.uk/humgen/genevar). The published expression data had been log-normalized using a quantile normalization method24 across experimental replicates for each individual, and then by a median normalization method across all 270 HapMap individuals.23 Only data for unrelated individuals were used in the analyses, resulting in a total of 60 US Caucasians of Northern and Western European origin (CEU), 60 Yoruba from Ibadan, Nigeria (YRI), 45 Japanese from the Tokyo area (JPT) and 45 Chinese from the Beijing area (CHB) in the final data set. For genes where multiple probes were available, analyses were performed using both individual probes and the mean expression across all probes. The differences were generally negligible and, for clarity, only the results based on mean probe intensities are presented in the paper. Probes were considered detectably expressed in LBLs if their average intensity was above the 25th percentile of all probes in that population.

Population-dependent expression analyses

Analysis of variance was used to assess population- and sex-dependent differences in gene expression, using either a single-variable (population) or multi-variable (population+sex) general linear model fit in the R statistics package (http://www.r-project.org). A Bonferroni-corrected P<0.001 was considered significant.

Genotype data and tag single-nucleotide polymorphism determination

Genotype data from the r27 release of combined data from HapMap Phase II and III were downloaded from the International HapMap Consortium web site (http://www.hapmap.org). For each population, local single-nucleotide polymorphisms (SNPs), that is, SNPs located within a region spanning from 50 kb upstream of the transcription start site to 50 kb downstream of the end of the longest transcript (in the UCSC Genome Browser human genome release hg18; http://genome.ucsc.edu) were extracted and used in subsequent analyses. Tag SNPs were determined for the CEU, YRI and combined JPT+CHB (ASI) populations separately, using the Tagger algorithm25 in HaploView 4.1 (http://www.broadinstitute.org/haploview/haploview) with r2⩾0.8 and minor allele frequency ⩾0.05.

Association analyses and multiple test corrections

We used a linear model, as implemented in PLINK v.1.0526 (http://pngu.mgh.harvard.edu/purcell/plink), to fit normalized log scale expression data for each gene to each of its local tag SNPs. Genotypes were encoded as 0, 1 or 2, depending on the number of minor alleles present. Two types of multiple test corrections were used: permutation correction, in which the order of gene expression levels to genotypes was randomly permutated 25 000 times, and a per-gene Bonferroni correction. A corrected P<0.05 was considered significant. The association analyses were performed in each of the three populations separately. To account for the Asian population being a combination of the Japanese and Han Chinese populations, conditional permutations were performed such that genotypes and phenotypes were only combined within each subpopulation.

Multiple SNP associations using partial least-squares projection

Partial least-squares (PLS) projection, as implemented in the PLS R package,27 was used to approximate the total expression level variance for each gene explained by its local SNPs (r2). To achieve comparable results for all genes and to minimize the risk of overfitting the regression models, we used only the first PLS component to estimate the total r2. The relative contribution of each local SNP to the total r2 was calculated as the fractional contribution of that SNP's variable importance to the total variable importance summed over all SNPs.

Prioritizing SNPs for experimental validation

All tag SNPs significantly associated with ABC gene expression levels were considered for experimental follow-up, as were all HapMap SNPs within their linkage disequilibrium (LD) blocks (r2⩾0.8). Multiple data sources were used to prioritize SNPs for experiments: an SNP was given a higher probability of being functional, that is, driving the observed expression variation, if it was (i) better correlated with the phenotype than were other SNPs in the LD block, (ii) located in regions with higher a priori probabilities of being functional, such as proximal to the transcription start site, in the 3′-untranslated region, or in the first intron, (iii) predicted to alter binding of transcription factors by the P-Match (http://www.gene-regulation.com), TF-Search (http://www.cbrc.jp/research/db/TFSEARCH.html) or the DeltaMatch software,28 or (iv) located in a region conserved between human beings and other placental mammals (based on the phyloP44wayPlacMammal track in the UCSC Genome Browser; http://genome.ucsc.edu).

Cloning of putative regulatory elements

The regions surrounding SNPs selected for experimental validation were amplified by a PCR of pooled human genomic DNA (SOPHIE cohort; http://pharmacogenetics.ucsf.edu), using previously reported PCR protocols.29, 30 PCR primers were designed using the Primer 3 software (http://frodo.wi.mit.edu/primer3) by inputting 350 bp upstream and downstream sequence relative to each SNP. The primers used to construct the putative regulatory regions are listed in Supplementary Table 3. The PCR-amplified regions (approximately 500 bp centered on the SNP) were purified by gel extraction (QIAquick Gel Extraction Kit’ Qiagen, Germantown, MD, USA) and cloned into pENTR-dTOPO vectors (Invitrogen, Carlsbad, CA, USA). Following plasmid DNA purification (QIAprep Spin Miniprep Kit; Qiagen), proper orientation was confirmed by sequencing. Clones in the correct orientation were subcloned into the multiple cloning site of the pGL4.23[luc2/minP] vector (Promega, Madison, WI, USA), upstream of a minimal TATA promoter and a luc2 firefly luciferase reporter gene. The orientation and sequence of the inserts were then re-verified by sequencing. The use of pooled genomic DNA from individuals of different ethnicities enabled us to obtain both the reference and variant alleles from the same PCR amplification for most regions. Site-directed mutagenesis (QuikChange Site-Directed Mutagenesis Kit, Stratagene, Santa Clara, CA, USA) was used to obtain the variant allele for the remaining regions. For follow-up studies of negative regulatory elements, putative regulatory regions were re-amplified from the pGL4.23[luc2/minP] plasmids and subcloned into the pGL3 Promoter vector (Promega), which contains a strong SV40 promoter and a luc+ firefly luciferase reporter gene downstream from the multiple cloning site.

Reporter assays for determining transcription regulatory effects

The human renal cell adenocarcinoma cell line ACHN was purchased from American Type Culture Collection (ATCC, Manassas, VA, USA). The human hepatoblastoma cell lines HepG2 and Huh-7, human colorectal cell line HCT-116 and human renal adenocarcinoma cell line HEK293T were supplied by the UCSF Cell Culture Facility. All cell lines were maintained in a culture medium consisting of Dulbecco's modified Eagle's medium with 4500 mg l−1 glucose, supplemented with 10% fetal bovine serum, 100 U ml−1 penicillin and 100 μg ml−1 streptomycin. In the reporter assay, the five cell lines were seeded into 48-well culture plates and were transfected 16–24 h later using the Lipofectamine LTX transfection reagent (Invitrogen), according to the manufacturer's protocols. The pGL4.23[luc2/minP] vectors (200 ng), with or without the region of interest, were co-transfected with pGL4.74[hRluc/TK] vector (5 ng; to normalize for transfection efficiency) into each well using 1.0–1.5 μl Lipofectamine LTX transfection reagent. Cells were lysed with passive lysis buffer 24 h after transfection, and were assayed for firefly and Renilla luciferase activity in a GloMax 96 Microplate Dual Injector Luminometer, using the Dual-Luciferase Reporter Assay system kit (Promega). To account for differences in cell type, transfection efficiency and cell density, reporter activities were normalized as log2 ((firefly luciferase/Renilla luciferase)/average of the negative control), where negative control refers to cells transfected with the pGL4.23[luc2/minP] vector or pGL3-Promoter vector not containing the PCR-amplified regions.

Accession numbers

The expression data used in this study have been deposited previously at the Wellcome Trust Sanger Institute web site (http://www.sanger.ac.uk/humgen/genevar) and in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo; accession no. GSE6536).

Results

Population-dependent expression of ABC transporters

Our analyses followed the steps outlined in Figure 1 and began with expression level analyses followed by association analyses of gene expression with genetic variants in and nearby each ABC transporter gene. For expression level analyses, first, we assessed the extent of inter-individual, population- and gender-dependent variation in the expression levels of each ABC transporter (Figure 1a). Genome-wide transcript levels, measured in LBLs from the 270 participants in the original HapMap European American (CEU), Han Chinese (CHB), Japanese (JPT) and Yoruba African (YRI) populations, were obtained from Stranger et al.23 Probes were available for all 49 human ABC transporters (Table 1), and 45 of these were detectably expressed in at least one of the HapMap populations.

Figure 1
figure 1

Schematic of the workflow applied to identify regulatory regions and single-nucleotide polymorphisms (SNPs) in ATP-binding cassette (ABC) genes. The number of SNPs assessed in each analysis step is shown in the shaded panel to the right. (a) Gene expression data for ABC genes were downloaded and extracted from the Wellcome Trust Sanger Institute GeneVar database. Expression levels for all 210 unrelated individuals were used to determine population differences in ABC gene expression. (b) Genotype data was downloaded from the HapMap database, and all SNPs within 50 kb of each ABC gene was used in the further analyses. (c) Tag SNPs were determined for each of the US Caucasians of Northern and Western European origin (CEU), Yoruba from Ibadan, Nigeria (YRI) and combined Japanese from the Tokyo area+Chinese from the Beijing area (ASI) populations separately. (d) Associations between ABC expression levels and SNPs were determined in each population. (e) Tag SNPs significantly associated with expression levels, together with all HapMap SNPs in linkage disequilibrium (r2⩾0.8) with them, were prioritized for experimental follow-up based on individual association significance, SNP location, species conservation and predicted altered transcription factor binding. (f) The ∼500 bp regions surrounding each selected SNP were cloned into a luciferase reporter vector (pGL4.23[luc2/minP] vector). (g) Luciferase activity was assessed in five different cell lines (ACHN, HEK293T, HCT-116, HepG2 and Huh-7) for each reporter construct. Constructs resulting in decreased luciferase activity compared with the negative control were subsequently cloned into an alternative vector carrying a strong SV40 promoter (pGL3-Promoter vector). (h) Allele-specific effects on luciferase activity were determined for all significant regulatory regions.

PowerPoint slide

Table 1 Number of genes expressed in lymphoblastoid cell lines from HapMap populations in the seven ABC transporter families

The inter-individual variance of ABC transporter transcript levels was slightly higher than the average for all detectable transcripts, although not statistically significant (1.2–1.5 × genome average; P<0.15–0.5, depending on the population). Compared to within-population variation, population differences contributed significantly less to the overall variability (Supplementary Table 1). On average, 10% (range: 1–40%) of the total variance for the different transporters was related to population. Expression levels generally correlated well (r2=0.985–0.997) between pairs of populations (Figure 2). To increase power, expression level data from the Japanese and Chinese populations were combined in subsequent analyses (ASI population).

Figure 2
figure 2

Comparisons of expression of ATP-binding cassette (ABC) transporter genes in lymphoblastoid cell lines among four populations. The correlations of normalized probe intensities are shown across all pairs of populations. Probe intensities are shown as filled circles, and the standard error of the mean as gray ellipses. For most probes, the standard errors were smaller than the size of the plot symbols.

PowerPoint slide

Of the 45 ABC transporters expressed in any population, 19 showed statistically significant population differences (Bonferroni-corrected P<0.001; Figure 3). However, these differences were generally modest, with an average difference in log-normalized expression of 0.27±0.18 between the populations with the lowest and highest expression (corresponding to an average relative difference of 21%; range 5–81%) (Figure 3). No major gender differences were observed (Supplementary Table 1); the only ABC transporter with a statistically detectable gender effect was ABCG1, which was expressed on average at 22% higher levels in female subjects than in male subjects (Bonferroni-corrected P<0.05).

Figure 3
figure 3

ATP-binding cassette (ABC) transporter genes exhibiting significant population differences in expression levels in lymphoblastoid cell lines. Analysis of variance was used to determine differential expression between the US Caucasians of Northern and Western European origin (CEU), Yoruba from Ibadan, Nigeria (YRI) and combined Japanese from the Tokyo area+Chinese from the Beijing area (ASI) populations. The plot shows ABC genes with significant differences after correcting for multiple tests (Bonferroni-corrected P<0.001). Bottom panels show log2 normalized expression values for each population. Top panels show relative differences compared with the overall mean across all populations, in original untransformed units.

PowerPoint slide

Associations between genotypes and gene expression

Next, we performed association analyses to identify local SNPs that associated with expression levels in each of the ABC transporters (see Figures 1b–d). Genome-wide analyses have suggested that the majority of regulatory polymorphisms are located within or proximal to the regulated genes,20, 22, 23, 31 and based on this, we decided to include SNPs in a ±50 000 bp region surrounding each transporter to balance coverage of putative regulatory regions against the statistical power to detect associations.

We used genotype data from the combined first, second and third phases of the HapMap project, which resulted in an increased coverage of genetic variation compared with previous genome-wide association analyses for this expression data set.23 In total, 13 242, 13 120 and 13 335 SNPs with a population-specific minor allele frequency of at least 5% were examined in the CEU, YRI and ASI populations, respectively. The number of SNPs per gene varied between 79 and 963 depending on the transporter and the population, but when accounting for the large variation in gene size, the coverage was more similar among genes and populations: between 0.54 and 4.71 SNPs per kb, with an average of 1.46, 1.45 and 1.43 in the Asian, Caucasian and African populations, respectively.

To increase our power to detect association signals, we exploited LD patterns to limit the number of statistical tests. As generally observed, genetic diversity was largest in the Yoruba population, for which a total of 3882 tag SNPs were needed to cover all common HapMap SNPs at an r2⩾0.8. In the Caucasian and Asian populations, the corresponding numbers were 2069 and 1744. The genome-wide association signals obtained for these tag SNPs are shown in Figure 4. Each population harbored some SNPs that were associated with expression levels of individual ABC transporters at extraordinarily low P-values. After correcting for multiple tests using a stringent Bonferroni correction or the permutation of phenotype data to determine empirical P-values, 24 tag SNPs were significantly associated with the expression levels of 10 different ABC transporters (Table 2). On average, the tag SNPs explained 21±8% (range: 11–38%) of the variation in the corresponding transcript levels. Figure 5a shows the most significant association for each of these 10 transporters.

Figure 4
figure 4

Genome-wide distribution of loci associated with ATP-binding cassette (ABC) gene expression. The uncorrected significance of each SNPs’ association with expression levels in lymphoblastoid cell lines in the US Caucasians of Northern and Western European origin (CEU), Yoruba from Ibadan, Nigeria (YRI) and combined Japanese from the Tokyo area+Chinese from the Beijing area (ASI) populations is shown as a function of relative genomic position, divided per chromosome. The relative locations of human ABC genes are shown in the top panel.

PowerPoint slide

Table 2 Summary of SNPs significantly associated with ABC gene expression in lymphoblastoid cell lines from three populations
Figure 5
figure 5

Single-nucleotide polymorphisms (SNPs) that significantly associate with expression levels of ATP-binding cassette (ABC) transporters in lymphoblastoid cell lines from three populations. (a) The most significant association across all populations is shown for each ABC gene with a significant local SNP association. (b) Population differences in ABCA12 expression (see Figure 3) can be attributed to a private SNP in the US Caucasians of Northern and Western European origin (CEU) population. Genotypes are shown for the coding strand.

PowerPoint slide

Notably, for some transporters, multiple significant tag SNPs were detected in a single population. For example, five significant associations were detected in ABCA1 in the CEU population, with r2 ranging from 0.21 to 0.38. As a cutoff r2 of 0.8 was used to define LD blocks, intermediate intercorrelations often still exist among selected tags (for example, the five significant tags in ABCA1 exhibit pairwise r2=0.03–0.51 and D′=0.31–1.0). Consequently, the tag SNPs are not completely independent, and partial correlations can result in cases in which a true causal effect of one tag is detected as a weaker effect in another. We therefore adopted a multivariate approach, using PLS projection, to estimate the overall contribution of all local SNPs to the expression variation.

For ABC transporters having at least one significant association, the full set of local SNPs explained 25±10% (range: 13–49%) of the total variance in expression levels (Supplementary Table 2). Very similar results were obtained when only tag SNPs were included in the models, indicating that the PLS projection methodology accounts for SNP intercorrelation caused by LD. The multi-SNP analyses typically recaptured the same relative importance of local SNPs as that determined through single SNP associations: the average correlation between PLS-derived SNP weights and correlation coefficients derived from the single SNP analyses was r2=0.92±0.06. In general, the total expression variance explained by all local SNPs was only moderately larger than that explained by the most significant single SNP (Supplementary Figure 1), suggesting that for most of the studied genes, a single regulatory locus explains all of the expression variation related to the interrogated local SNPs.

Population-dependent associations

We examined whether SNPs that associated with expression levels in any one population also explained the observed population differences. For seven of the SNPs associated with within-population expression variation, the allele frequencies were highly correlated with population average transcript levels (r2⩾0.95; P<0.1). In particular, the significantly higher expression of ABCA12 in the CEU compared with the YRI and ASI individuals could be attributed to an SNP (rs2888327) private to the Caucasian population (Figure 5b). Similarly, population average expression levels of ABCC10 were well explained by the allele frequencies of rs12195350 (r2=1.0; P<0.002), with expression levels CEU>ASI>YRI and allele frequencies CEU>ASI>YRI (data not shown).

Detection of transcriptional regulatory regions

SNPs that associated with the expression level of the transporters were used to identify putative regulatory regions. These regions were then experimentally validated in medium-throughput functional studies using a transient transfection reporter system (see Figures 1e–h). We prioritized all SNPs in LD with each significant tag based on multiple criteria, including their association with additional phenotypes, their location in genomic regions with higher a priori potential for regulatory function (for example, near the transcription start site, in untranslated regions, or in the first few introns of the gene of interest), conservation to other placental mammals or predicted alteration of transcription factor motifs (Figure 1e; see Materials and methods section).

A total of 16 SNPs in seven ABC transporters were selected for experimental verification. For each of these, a region of approximately 500 bp centered on the SNP was cloned from human genomic DNA into a luciferase reporter vector, upstream of a minimal TATA promoter. The constructs were transiently transfected on three to five separate occasions into five different cell lines, derived from human colorectal epithelium (HCT-116), renal epithelium (ACHN, HEK293T) or hepatocytes (HepG2, Huh-7) (Figure 6a). For most constructs, similar results were obtained across the panel of cell lines. Five of the selected 16 variants resulted in a more than twofold increased luciferase activity in multiple cell lines, three of which showed more than fourfold increased activity in specific cell types. Notably, several of the cloned regions resulted in a markedly reduced luciferase activity (Figure 6a). As the weak promoter is not specifically designed to study negative regulatory effects, we subcloned these six regions into an alternative vector upstream of the stronger SV40 promoter. Two of the regions gave a more than fourfold reduced luciferase activity also for the strong promoter (Figure 6b), suggesting a binding of repressive regulatory factors to these regions. In summary, guided by statistical associations between genetic polymorphisms and expression levels, we identified and validated a total of five regions that were designated enhancer or silencer regions, and two that were designated moderate enhancers. Through affecting the expression levels of proximally located ABC transporters, these regulatory regions may influence diverse physiological and pharmacological functions, including cholesterol efflux to high-density lipoprotein,3, 4 peroxisomal transport of very-long-chain fatty acids,9, 32 presentation of cytosolic peptide antigens to major histocompatibility complex class I molecules in the endoplasmatic reticulum33, 34 and resistance to anticancer drugs such as docetaxel and vinorelbine35, 36 (Table 3).

Figure 6
figure 6

Regulation of luciferase reporter activity by genomic regions surrounding single-nucleotide polymorphisms (SNPs) in or near ATP-binding cassette (ABC) transporter genes. (a) Putative regulatory elements were cloned into the weak promoter luciferase vector pGL4.23[luc2/minP]. (b) Regions resulting in negative regulation in the weak promoter assay were assayed in an alternative vector containing the strong SV40 promoter (pGL3[Promoter]). Luciferase activity is shown relative to that of a vector only containing the promoter and the luciferase gene. The activity of each region was assessed in multiple cell lines (from left to right, for each construct: ACHN, HCT-116, HEK293T, HepG2 and Huh-7; the latter was only used in the weak promoter assay). Red, yellow, light blue and dark blue bars: luciferase activity >4 × , >2 × , <0.5 × and <0.25 × background promoter activity, respectively. Values are shown as mean±standard error of the mean from three to five separate experiments, each performed in triplicate.

PowerPoint slide

Table 3 Regulatory regions in ABC transporter genes identified experimentally

For the seven regulatory regions (five of which resulted in >2-fold enhanced luciferase activity and two of which resulted in a >4-fold reduced luciferase activity), we compared luciferase activities of the reference and variant alleles to determine the effect of the variants on transcription. Significant differences were observed for five of the allelic variants (Table 4). For most of these, statistical differences were specific to one or two of the tested cell lines. Four polymorphisms showed the same directionality as that observed in the association with mRNA expression in LBLs, although for two of these the difference did not reach statistical significance. Notably, the SNP rs17091297 in ABCD1 showed large differences across all tested cell lines, reproducing the directionality of the LBL association: the major allele resulted in an almost eightfold lower reporter activity than the minor allele, mirroring the lower average LBL gene expression for the major allele (Figure 7). In contrast, the SNPs rs12227668 in ABCD2 and rs9349256 in ABCC10 showed the opposite directionality to that expected from LBL expression levels (Table 4); these results were consistent across multiple experiments and cell lines. Computationally predicted transcription factor binding affinities were consistently greater for the genetic variants that showed the higher reporter activity (Table 4), suggesting that the observed allelic differences are real, and that the different directionality could be the result of combinatorial regulation by additional regions in vivo, which were not included in the cloned regions.

Table 4 Direction of association of selected M versus m alleles of ABC transporter genes with expression levels in lymphoblastoid cell lines and with luciferase activities
Figure 7
figure 7

Association of a single-nucleotide polymorphism (SNP) in ABCD1 with expression levels in lymphoblastoid cell lines from the Yoruba from Ibadan, Nigeria (YRI) population, and with luciferase reporter activity. (a) The SNP rs17091297 is correlated with ABCD1 mRNA levels in the YRI population (Bonferroni-corrected P<0.05). (b) The major allele results in decreased luciferase activity compared with the control, when inserted upstream of a strong SV40 promoter. Luciferase activity is shown relative to the background promoter activity. Values are shown as mean±standard error of the mean from four separate experiments, each performed in triplicate.

PowerPoint slide

Thus, by using statistical associations between SNPs and expression levels to guide the experiments, five regions were identified that resulted in a markedly altered reporter activity, and an additional six having intermediate effects. To assess if this was more than expected by chance, we estimated the total number of regulatory regions in the examined ABC transporters based on a combination of species conservation data and experimental data for DNAse I hypersensitivity or histone 3-lysine 4 monomethylation (H3K4me1). Each of these measures have been shown to correlate with regulatory function,37, 38, 39, 40 but used separately they are likely to overestimate the number of regulatory active sites. We used the fraction of the interrogated genomic DNA that had two or more of these characteristics as an estimate of true regulatory function, resulting in a background chance of hitting a regulatory region of 6.4% (Supplementary Table 4). In contrast, our rate of 29–69%, depending on the cutoff used to call a significantly altered reporter activity, is significantly higher (P<0.0002–1 × 10−14).

Discussion

Predisposition to disease and the response to drugs are highly variable in human populations,41, 42 with expression level differences in various genes contributing to interindividual differences in disease risk and drug response.43, 44 We assessed the influence of ethnicity, gender and local genetic variation on the expression levels of ABC transporters in LBLs, and used the detected eQTL to guide experimental assays of putative regulatory genomic regions.

For 19 of the 49 ABC genes in the human genome, ethnicity was significantly associated with expression levels. As ABC transporters are known to affect the disposition of many drugs, metabolites and endogenous ligands,1, 2 such population-dependent expression can ultimately lead to differential susceptibility to disease or outcome of drug therapy among populations. For example, we observed that the expression level of the hepatic bile salt efflux pump (ABCB11/BSEP), which when mutated is associated with familial, drug- and pregnancy-induced intrahepatic cholestasis,13, 45, 46 was significantly lower in the Asian than in the European and African samples. Although our data cannot be directly extrapolated to expression level differences in the liver, if such population differences in hepatic expression levels of ABCB11 occurred, our data would suggest that Asians may be more susceptible to intrahepatic cholestasis mediated by BSEP. Intriguingly, the incidence of pregnancy-induced cholestasis is higher in Asians than in European women, and the possible connection with population-dependent BSEP expression thus warrants further studies.45

The largest difference among populations was observed for ABCG1, with expression levels on average 1.8-fold higher in the Asian compared with the Caucasian samples. Notably, in contrast to the several cases of statistically significant population-dependent expression, ABCG1 was the only ABC transporter for which gender was significantly associated with expression. ABCG1 mediates the efflux of cellular cholesterol to high-density lipoprotein,5, 47 and the higher average expression levels in women, and in the Asian population, suggest that ABCG1 expression might contribute to the relatively higher high-density lipoprotein plasma levels observed in these populations.48, 49 Again, some caution in extrapolation of our data in LBLs to liver must be exercised.

Previous studies have indicated that genetic differences account for a substantial part of observed gene expression differences among populations.50, 51 This was also observed here for some of our genes, with population differences in allele frequencies generally well correlated with population average expression. Particularly, the SNP rs2888327 was only observed in the CEU population, which had a higher average expression of the proximally located ABCA12 gene than the ASI and YRI populations (Figure 5b). On the whole, however, the contribution of population and gender was small compared with the total ABC gene expression variation, consistent with other studies of population effects on gene expression,50, 51, 52 and also with the observation that 85–95% of human genetic variation is due to within-population variation.53, 54

With current advances in sequencing techniques, access to near-complete personal genomes is not far from reality. However, for these large-scale data sets to be truly useful, much work is needed in functionally annotating the human genome. So far, our knowledge of non-coding regions is greatly lagging behind that of protein coding regions. Especially little is known about the extent and function of distantly located regulatory elements such as enhancers and silencers. Evolutionary constraints on non-coding sequences have been successfully used to identify regions with regulatory function.17, 18, 19 However, not all evolutionarily conserved regions have detectable functions in in vitro and in vivo enhancer assays.18, 19, 55 Furthermore, many regulatory elements are likely located in regions not detectable using comparative genomics.56, 57 Alternatively, DNAse I hypersensitivity,37, 38 binding of enhancer-associated proteins such as the acetyltransferase and transcriptional coactivator p30040, 55 and epigenetic modifications such as histone methylation and acetylation39, 40 have been used to predict the genome-wide locations of regulatory regions. Although highly promising, additional experimental validation of predicted regulatory regions is needed to determine the predictive value of these methods.

Here, we used a complementary approach to identify putative regulatory regions, through the association of genetic polymorphism with expression levels of ABC transporters. We identified 24 tag SNPs that were significantly associated with ABC transporter expression levels in LBLs. These were in LD with a total of 104 interrogated SNPs with a minor allele frequency >5% in at least one population. After prioritizing the initial list of possible eQTLs based on association strength, genomic location, predicted alteration of transcription factor binding sites and species conservation, 16 regions within a ±50 000 bp range surrounding seven ABC transporters were selected for experimental validation. This eQTL-centered approach led to the identification of five genomic regions that resulted in at least fourfold increased or decreased reporter activity, and six additional regions with at least twofold altered activities. Allele-specific regulatory function was shown for five of the regions. Genetic variation in these regions can thus contribute to differential expression among individuals, and can ultimately affect drug response or susceptibility to disease.

Previous analyses of expression quantitative traits in LBLs20, 22 have focused on global properties of regulatory polymorphisms, and have not included experimental verification of identified associations. In our gene family focused effort, 5–11 of 16 tested associations were experimentally replicated depending on the reporter assay cutoff used. However, five regions did not alter reporter activity. Although significant signals in the reporter assay suggest causative regulatory effects of the examined regions, interpreting the lack of signal is less straightforward. It could be a consequence of false positives in the association analysis, but could also result from LD between the examined region and the actual causative site. In addition, combinatorial regulation of gene expression in vivo is difficult to replicate in an artificial system that does not include the larger genomic context.

We recently used a combination of comparative genomics and transcription factor binding site predictions to select 50 putative regulatory regions in important hepatic drug transporters from the ABC and SLC families.58 The enhancer activities of these regions were assayed in mice using hydrodynamic tail vein injection, which results in a highly localized expression in the liver. Using this technique, 12/50 (24%) of the prioritized regions showed experimental activity in vivo. Notably, however, none of the regions with experimental regulatory impact in this study were detectably conserved among placental mammals, suggesting that eQTL mapping provides additional information compared with comparative genomics approaches.

This work thus shows the general utility of expression level QTL in defining genomic regions with regulatory potential, and, specifically, shows transcriptional regulation of human ABC transporters by several hitherto uncharacterized genomic regions. Genetic and epigenetic variation in the identified regions can thus lead to altered transport of drugs and their metabolites, and of endogenous solutes such as cholesterol, fatty acids and peptide fragments. Ultimately, this can cause inter-individual variation in the pharmacological response to drugs and in important physiological processes, including cellular lipid metabolism, lysosomal protein degradation and in antigen presentation on the cell surface.