DNA Sequence Variation in a 3.7-kb Noncoding Sequence 5′ of the CYP1A2 Gene: Implications for Human Population History and Natural Selection

Subjects

DNA sequence was obtained from 113 individuals, together referred to herein as the “Old World sample”: 33 Africans (3 Biaka Pygmies, 7 Mbuti Pygmies, 3 Alur, 3 Nande, 3 Hema, 2 Kenyans, 3 Nigerians, 3 Sotho, 3 Nguni, and 3 San), 25 non-Indian Asians (hereafter termed “Asians”) (3 Cambodians, 5 Han Chinese, 5 Japanese, 3 aboriginal Malayans, 3 Vietnamese, 2 Mongolians, and 4 Daghestanis), 27 Europeans (5 northern Europeans, 5 French, 5 Italians, 5 Poles, and 7 Finns), and 28 South Indians (3 Brahmin, 3 Kshatriya, 3 Yadava, 3 Mala, 3 Madiga, 3 Relli, 3 Irula, 3 Khondadora, 1 Santal, and 3 Chenchu). The Khondadora, Santal, and Chenchu individuals are members of South Indian “tribal” populations, whereas the other Indians are members of South Indian caste populations. Informed consent was obtained from all subjects whose blood was drawn either at the University of Utah or in India.

Sequencing was also performed on the 90-member PDR sample (Collins et al. 1998). This collection of individuals consists of European Americans, African Americans, Hispanic Americans, Native Americans, and Asian Americans. DNA sequencing was also performed for three ape species: five common chimpanzees (Pan troglodytes), one bonobo (P. paniscus), and one gorilla (Gorilla gorilla).

Laboratory Analysis

PCR primer sequences (5′-AACAGGGACTTCTTGGATGC-3′ and 5′-TGTACCAAAGAGTCCCTGCC-3′) were derived from the human CYP1A2 genomic sequence (GenBank accession U02993), spanning 0.5 kb downstream and 3.2 kb upstream from the CYP1A2 promoter. PCR amplification was performed in 50-μl reaction volumes by the Expand Long Template PCR System (Roche). Each reaction contained 50 ng genomic DNA, 350 μM deoxyribonucleoside triphosphates (dNTPs), 0.2 μM each PCR primer, 1 × reaction buffer 2, and 2.6 U Taq/Pwo polymerase mix. Cycling conditions included an initial denaturation at 94°C for 2 min; 10 cycles at 94°C for 10 s, 55°C for 10 s, 68°C for 2 min; and 15 cycles at 94°C for 10 s, 55°C for 10 s, 68°C for 2 min 20 s. Residual primers and dNTPs were removed from PCR products with a Millipore glass fiber filter. The sequence-ready templates were eluted in 70 μl of sterile H₂O. Five microliters of each template was aliquoted to 12 wells of a 384-well sequence dish and was evaporated to dryness in a speed-vac.

Internal primers used for sequencing included the following: 5′-GGATGCTTATGATGTCTCTTG-3′, 5′-GTGCGTGTCAGGTCTCTTCA-3′, 5′-AGAAGGAGCGTAATCCCCC-3′, 5′-AACCTGTGAAGATGCCAAGG-3′, 5′-ACCGAGCCTAACCTCAAACC-3′, 5′-TGTAGAGAGGAGGTCTTG-3′, 5′-TTGTGGTCCCAGCTACTC-3′, 5′-GTTTATCCTTGCTTGAGGG-3′, 5′-CCCCTCAAGCAAGGATAAAC-3′, 5′-GGGACATGTAAGCACGGACT-3′, 5′-GGGGATCATGACACTTCCAT-3′, 5′-ATCAGATTGGCCTGGTTGTC-3′, and 5′-CATGTACCTTCATCCCCAGG-3′. Cycle sequencing was performed in 5-μl reaction volumes by use of ABI BigDye Terminator chemistry. Cycling conditions included an initial denaturation at 96°C for 30 s; followed by 46 cycles of 96°C for 10 s, 50°C for 5 s, and 60°C for 4 min. On completion of cycle sequencing, 20 μl of 62.5% ethyl alcohol and 1 M potassium acetate (with pH 4.5) was added to each reaction, and the sequence plates were centrifuged at 4,000 rpm at 4°C for 45 min in an Eppendorf centrifuge. The samples were resuspended in 15 μl of sterile H₂O and were electrophoresed on an ABI 3700 DNA analyzer prepared with POP-5 capillary gel matrix (ABI).

Sequence-trace files were evaluated using the Phred, Phrap, and Consed programs (Ewing et al. 1998). Potential heterozygotes were identified using PolyPhred, version 3.5 (Nickerson et al. 1997). Polymorphisms were verified by manual evaluation of the individual sequence traces. For most polymorphisms, it was possible to evaluate both the forward and reverse sequences. Detailed information about the SNPs typed in the PDR sample is available at the Web sites dbSNP Home Page and GeneSNPs Public Internet Resource.

Statistical Analysis

Allele frequencies for each SNP were determined by gene counting, and the significance of deviations from Hardy-Weinberg equilibrium was tested using the random-permutation procedure implemented in the Arlequin package (Schneider et al. 2000). Haplotype frequencies were estimated separately for each major population (Africans, Indians, Asians, and Europeans) by the expectation-maximization (EM) algorithm provided in the Arlequin package. This procedure has been shown to yield reliable estimates of haplotype frequencies (Tishkoff et al. 2000). In particular, the EM approach should work well in small regions that have both a high degree of linkage disequilibrium and a high proportion of homozygotes. The data analyzed here meet both of these conditions. Unambiguous haplotypes were inferred directly from genotypes in which only zero or one sites were polymorphic. The inferred haplotype data were used to estimate the D′ linkage-disequilibrium statistic (Lewontin 1964) for all pairs of SNPs.

Haplotypes were used to perform phylogenetic analyses, with the common chimpanzee as an outgroup species. A minimum-spanning tree (MST) relating haplotypes was constructed using the Arlequin package. Trees were also constructed using the neighbor-joining, parsimony, and maximum-likelihood methods implemented in the PHYLIP software package (Felsenstein 1993).

The recombination rate for the genomic region that includes CYP1A2 was estimated by comparison of genetic and physical distances between polymorphic STR markers in a 5-Mb region flanking the gene. The genetic and physical distances between three independent marker pairs were determined using Marshfield and Généthon maps, from NCBI Map Viewier (build 29).

The nucleotide-substitution rate of the CYP1A2 5′ region was calculated in two steps. First, the expected coalescence time of human and chimpanzee lineages was calculated. The expected coalescence time of human and chimpanzee lineages is equal to the sum of the time since human/chimpanzee speciation and the expected coalescence time in the population ancestral to humans and chimpanzees. The time since human/chimpanzee speciation is often estimated at ∼5 million years, and the coalescence time in the population ancestral to humans and chimpanzees has been estimated at 4 million years (Takahata et al. 1995; Sherry et al. 1997). Second, the mean nucleotide difference between humans and chimpanzees was divided by twice the expected coalescence time, to yield a nucleotide-substitution rate measured in substitutions per year. To provide a basis for comparison with previous studies, we estimated the nucleotide-substitution rate by calculating the mean nucleotide difference between human and chimpanzee haplotypes and dividing by twice the estimated time of speciation, 5 million years before the present.

Average nucleotide diversity (π) was assessed in the PDR sample, the Old World sample, and each major population (Africans, Indians, Asians, and Europeans). Because of their distinct geographic location and history, the samples from the Indian subcontinent were treated separately from those from the rest of Asia. Genetic diversity under a mutation-drift equilibrium model (θ) was assessed using the methods of Tajima (1983) and Watterson (1975). These values were in turn used to estimate effective population size, N_e, using the standard equation, θ=4N_eμ. The μ value is estimated as νgL, where g is generation length, L is the length of the nucleotide sequence (3,669 bp), and ν is the mutation rate per nucleotide per year. We assume a generation time of 20 years throughout. This is lower than current human generation times but is probably reasonable for most of human history (Chen and Li 2001).

The selective neutrality of mutations in these sequences was tested using the methods of Tajima (1989), Fu and Li (1993), and Fay and Wu (2000, 2002). Statistical significance was based on 5,000 simulated samples. Hudson/Kreitman/Aguadé (HKA) tests (Hudson et al. 1987) were used to compare diversity patterns in the CYP1A2 5′ region with diversity patterns found at eight other loci: HOXB6 (Deinard and Kidd 1999), ZFY (Dorit et al. 1995), β-globin (Fullerton et al. 1994), DMD (Nachman and Crowell 2000a), ND2 (Wise et al. 1998), CCR5 (Bamshad et al., in press), Xq13.3 (Kaessmann et al. 1999), and three loci on the Y chromosome (Thompson et al. 2000). HKA tests were performed using the HKA computer program (Hey 2001).

Population History

The nucleotide-diversity values found for the CYP1A2 5′ region are similar to those found in most other surveys of human autosomal DNA sequences (Halushka et al. 1999; Wall and Przeworski 2000; Zhao et al. 2000; Jorde et al. 2001). In a survey of SNP densities in 200-kb bins spread across the human genome, the International SNP Map Working Group (2001) found a modal π value of ∼0.0007, with a range from 0.0000 to 0.0020. The value that we find for the CYP1A2 5′ region in the Old World sample (π=0.00043) falls at approximately the 20th percentile in the International SNP Map Working Group’s distribution (see fig. 2a in International SNP Map Working Group 2001). The π value that we observe is only slightly lower than 0.00054, which is the average reported by Stephens et al. (2001) in a survey of >300 genes. The low π values generally observed in humans support the conclusion that humans are genetically less diverse than many other species (Li and Sadler 1991; Nachman et al. 1998; Kaessmann et al. 1999). Under the assumption that evolution in the CYP1A2 5′ region has been neutral, the relatively low level of nucleotide diversity predicts low effective population sizes, estimated here as ∼5,000–15,000. These values are similar to values obtained in previous studies of protein polymorphisms (Nei and Graur 1984), mtDNA (Takahata 1993), Y-chromosome DNA (Hammer 1995; Goldstein et al. 1996), X-chromosome polymorphisms (Zietkiewicz et al. 1998; Kaessmann et al. 1999), and autosomal polymorphisms (Harding et al. 1997; Sherry et al. 1997; Halushka et al. 1999; Zhao et al. 2000).

Nucleotide-diversity values were higher in Africa than elsewhere, and Africans had the largest number of haplotypes. Greater African diversity has been observed in many studies of human genetic diversity—including those of protein polymorphisms (Nei et al. 1993), mtDNA (Vigilant et al. 1991), Y-chromosome DNA (Hammer et al. 1997; Seielstad et al. 1999; Jorde et al. 2000), autosomal microsatellites (Bowcock et al. 1994; Deka et al. 1995; Calafell et al. 1997; Jorde et al. 1997), autosomal Alu polymorphisms (Jorde et al. 2000), and nuclear SNPs and DNA sequences (Nickerson et al. 1998; Halushka et al. 1999; Kaessmann et al. 1999; Rieder et al. 1999; Labuda et al. 2000; Wall and Przeworski 2000; Zhao et al. 2000; Jorde et al. 2001; Yu et al. 2001). This pattern is consistent with an African origin of modern humans, but it can also be explained by a higher effective population size in Africa (Relethford and Jorde 1999).

Additional evidence for an African origin is provided by the geographical distribution of haplotypes: haplotypes found outside Africa are mostly a subset of those found within Africa. This pattern has been observed in a number of other studies (Tishkoff et al. 1996, 2000; Calafell et al. 1997; Hammer et al. 1997; Watson et al. 1997; Alonso and Armour 2001) and agrees with a scenario in which some African populations left the continent to populate the rest of the Old World (Jorde et al. 1998; Pääbo 1999; Yu et al. 2002).

The linkage-disequilibrium patterns seen in this study are in agreement with a number of other studies that have found less linkage disequilibrium in African than in non-African populations (Tishkoff et al. 1996, 1998, 2000; K. K. Kidd et al. 1998; J. R. Kidd et al. 2000; Reich et al. 2001). Because linkage disequilibrium decreases through time, as a result of recombination, levels of disequilibrium can be correlated with the relative “age” of a population, with older populations having less disequilibrium. Thus, the linkage-disequilibrium results are also consistent with an African origin of modern humans.

Tajima’s D test failed to reject the hypothesis of neutrality and constant population size. However, the D and F tests of Fu and Li, which are statistically more powerful under most circumstances (Fu 1996), yielded results that indicate a significant excess of rare variants. This departure could be the result of either positive natural selection or an expansion from a very small population to a much larger one. The latter possibility is expected to produce an excess of rare alleles, because genetic drift, which tends to remove rare alleles from populations, has little effect during and after an expansion (Slatkin and Hudson 1991; Rogers and Jorde 1995).

The results of Fu’s D and F tests are consistent with several lines of evidence for a population expansion in humans. Among genetic systems, support for a human population expansion has been derived from studies of mtDNA (Di Rienzo and Wilson 1991; Rogers and Harpending 1992), Y-chromosome microsatellites (Pritchard et al. 1999) and SNPs (Shen et al. 2000), and autosomal microsatellites (Shriver et al. 1997; Di Rienzo et al. 1998; Kimmel et al. 1998; Reich and Goldstein 1998; Gonser et al. 2000; Jin et al. 2000; Zhivotovsky et al. 2000). Details vary from study to study—for example, the estimated dates of the expansion differ in some cases, and some studies find evidence for expansions in specific populations, such as Africans, whereas others do not. However, analyses of these different types of genetic systems all detect the signature of a population expansion sometime during the past 100,000 years.

Diversity patterns in the CYP1A2 5′ region combine to form a picture of human genetic diversity that is unremarkable in a number of respects: Africans are more diverse, non-African diversity is largely a subset of African diversity, and evidence for population growth is present. However, some aspects of this variation suggest that the history of CYP1A2 is characterized by more-complex factors as well.

Natural Selection

The failure of Tajima’s D test to reject the hypothesis of neutrality and the presence of significantly negative D and F statistics in Fu and Li’s test are surprising given the large number of high-frequency variants. More than 20% of sites in the Old World sample have frequencies >90%. Under neutral expectations, <5% of sites should attain such frequencies (fig. 1).

The failure of the D and F tests to detect an excess of high-frequency variants is probably due to the insensitivity that these tests have to such variants. Tajima’s D statistic, which compares the number of segregating sites in a sample with the mean pairwise difference (Tajima 1989), is insensitive to the presence of very-high-frequency variants because they contribute little to mean pairwise differences. In a sample of 100 sequences, for example, a derived variant with a frequency of 1/100 will contribute the same to mean pairwise difference as a derived variant with a frequency of 99/100. Fu and Li’s statistics, which compare the number of singletons with the number of nonsingletons (Fu and Li 1993), are insensitive to the presence of high-frequency variants because they are scored simply as nonsingletons. Thus, all three of these statistics are probably unaffected by extreme frequency spectra, such as that observed in the CYP1A2 5′ region.

The application of Fay and Wu’s H test, which is sensitive to the presence of high-frequency derived variants, rejects the hypothesis of neutrality in our data (table 6) (Fay and Wu 2000; Fay 2002). This result is consistent with the observations that two distinct haplotype clusters are present and that the frequency of the more-derived cluster (A) is higher than expected.

The hypothesis of neutrality in the CYP1A2 5′ region is placed further in doubt by the HKA tests, which indicate that within- and between-species diversity patterns in the CYP1A2 5′ region are similar to those of HOXB6, ZFY, and β-globin. This implies that diversity patterns in the CYP1A2 5′ region are similar to those in coding regions with known functions. The results of HKA tests do not conclusively exclude either natural selection or evolutionary neutrality in the CYP1A2 5′ region. Unless HOXB6, ZFY, and β-globin have been evolving neutrally, however, similarities between the CYP1A2 5′ region and these regions make neutrality in the CYP1A2 5′ region seem unlikely.

One factor that can cause the presence of distinct haplotype clusters is population subdivision, which can allow divergent clusters to evolve in different demes. Under these conditions, divergent clusters should be restricted to different geographical areas. This is not the case in the CYP1A2 5′ region. Haplotype clusters A and B both are widespread geographically (fig. 2). Both are found in several Asian populations (Cambodians, Chinese, Japanese, and Indians) and several African populations (Sotho, Biaka Pygmy, Mbuti Pygmy, Nande, Hema, and Nigerians). The presence of the two clusters does not therefore seem to be due to population subdivision.

Another factor that can cause excesses of high-frequency variants is balancing natural selection. Balancing natural selection can maintain old evolutionary lineages at low frequencies by protecting them from genetic drift (Lewontin and Hubby 1966). Positive selection with low levels of recombination can have similar effects (see fig. 2 in Fay and Wu 2000). Positive selection can have the additional effect of reducing the overall levels of genetic variation, which could explain the relatively low levels we observe in the CYP1A2 5′ region. These low levels of variation are also consistent with the reduced recombination observed in the region (Aquadro et al. 2001; Nachman 2001). Finally, the presence of four high-frequency derived SNPs at frequencies >0.9 matches patterns observed by Fay and Wu (2000) in simulations of positive natural selection. Such selection would not need to be acting directly on the CYP1A2 5′ region. Ding et al. (2002), for example, found positive selection in markers linked to DRD4. The variation patterns that we see could be caused by the action of natural selection on CYP1A2 exons or even on other loci (Braverman et al. 1995; Fay and Wu 2000).

The hypotheses of population growth and natural selection are not mutually exclusive. It is possible for both to occur at the same time in the same populations. Whether both occur in the CYP1A2 5′ region and relative importance of each are still unclear. Further analyses of diversity in CYP1A2, especially in coding regions, may provide more information about the relative importance of demographic and selective effects by allowing comparisons between synonymous and nonsynonymous nucleotide substitutions (Fay et al. 2001).

Implications

Patterns of population genetic variation in the CYP1A2 5′ region shed some light on questions, about the relationship between CYP1A2 diversity, ethnicity, and geography, raised by earlier association studies. There are appreciable regional differences in allele frequency. Haplotype 9, for example, was found in ∼20% of Asian and African samples but was rare in the European and Indian samples. Furthermore, although haplotype cluster B was found with a frequency of ∼10% in the Africans, Indians, and Asians, Europeans lacked haplotypes from cluster B altogether.

Our sequences do not extend far enough upstream or downstream from the CYP1A2 promoter to allow comparisons with most polymorphisms identified in intron 1 (Sachse et al. 1999) and in 5′ regions (Aitchison et al. 2000). However, one of the polymorphic sites identified in our sample, position −2963, was studied by Nakajima et al. (1999), who found the A nucleotide at a frequency of 0.23 in a Japanese study population. We find the same variant at a similar frequency (0.22) in the Asian population that we studied. On the basis of both analyses of rates of caffeine metabolism in Japanese subjects and gel-retardation assays that detected differential protein binding to the A and G variants at position −2963, Nakajima et al. (1999) concluded that variation at that position accounts for some of the phenotypic variation in CYP1A2 expression. Our results do not contradict those of Nakajima et al. (1999), but they do suggest that other sites may be important as well. The sites that distinguish haplotype clusters A and B may be of particular interest.

Diversity patterns in the CYP1A2 5′ region also suggest a number of hypotheses about phenotypic variation in CYP1A2 expression and activity. To date, most association studies on CYP1A2 have compared the mean values of phenotypes, rather than their variances. However, the geographic distribution of genetic diversity that we observe suggests that phenotypic variances should be analyzed as well. One hypothesis is that, if haplotypes in the CYP1A2 5′ region are associated with phenotypes, then phenotypic variance should be greatest in Africans, followed by Asians and then Europeans.

A simple test of the hypothesis that CYP1A2-related phenotypes are more diverse in Africans than in Europeans can be performed with the data of Basile et al. (2000), who used the abnormal-involuntary-movement scale (AIMS) to measure phenotypes in African Americans and Europeans who were treated with antipsychotic drugs thought to interact with CYP1A2 to cause tardive dyskinesia. Basile et al. found significant differences between individuals with genotypes defined by an A/C SNP (Basile et al. 2000, but see Schulze et al. 2002). Within genotypes, Basile et al. (2000) found no significant difference between the mean AIMS scores of African Americans and Europeans. We found, however, that the phenotypic variance measured using the AIMS is significantly greater in African Americans than in Europeans, in two of the three reported genotypes. An analysis of variance yields P<.01 for Basile et al.’s AA and AC types (see fig. 2 in Basile et al. 2000). The greater phenotypic variance found in Africans could be associated with the greater diversity in African CYP1A2 haplotypes.

Whether clusters A and B contain functionally divergent haplotypes cannot be determined directly from our analyses. More important is the fact that diversity patterns in the CYP1A2 5′ region reflect patterns of population genetic variation, including broad-scale similarities and differences among geographical regions, that can be used to improve the design and implementation of future studies. The presence of these clusters is likely to be of some significance in further studies of the association between CYP1A2 genotypes and phenotypes.

DNA sequences from chimpanzee and gorilla were a key source of information in our analyses. We anticipated that inclusion of these species would resolve the ambiguity of ancestral and derived character states in only a small number of variable positions. However, these outgroup comparisons were crucial in the identification of high-frequency variants. The usual assumption in the absence of outgroup data is that the minor allele is in the derived state. With such an assumption, we would have described all four high-frequency variants as low-frequency variants, and the outcomes of statistical tests would have been flawed. The importance of outgroup comparisons is well recognized, and chimpanzees are becoming widely adopted as a standard in human genetic analyses (Hacia et al. 1999). Our results reemphasize the importance of the inclusion of outgroups in analyses of human genetic variation.

The PDR is being extensively used to detect new human SNPs, so the comparison between results based on this sample and those based on the Old World sample is of some significance. It is expected that high-frequency alleles should be ancient and thus should be widely shared among populations. This expectation was confirmed in the comparison undertaken here: the more common polymorphisms were indeed found in both samples, indicating that the PDR sample should be effective for the detection of most such polymorphisms. One polymorphism, with a frequency of 11% in the African sample, was not, however, detected in the PDR sample. This likely reflects a broader sampling of African diversity in the Old World sample than in the PDR sample.

The wider representation of diversity in the Old World sample also resulted in higher average nucleotide diversity than in the PDR sample and in a larger number of detected low-frequency polymorphisms (although the latter is also partially due to a larger number of chromosomes in the Old World sample [226] than in the PDR sample [180]). Although low-frequency polymorphisms are less informative for linkage and linkage-disequilibrium analyses (Ott and Rabinowitz 1997), they can sometimes be useful for association-based mapping. For example, if a disease-causing mutation occurs on a chromosome that contains a rare marker, the disease-marker association will be stronger than it would be if the mutation were to occur on a chromosome that contains a polymorphism that is common in the population (Jorde 2000).

Patterns of genetic diversity observed in the CYP1A2 5′ region bring attention to some basic limitations of the PDR as a tool for the study of patterns of human genetic variation. The PDR is intended primarily to “facilitate the discovery of DNA sequence variants” (Collins et al. 1998, p. 1229), not for population genetic analysis. We found that, although patterns of genetic diversity in the Old World and PDR samples are similar in some respects, a clear advantage of the Old World sample is that genetic diversity can be measured and compared in different subpopulations. Between-population comparisons can be an important step in the characterization of genetic diversity, as was demonstrated by our comparisons of haplotype-cluster frequencies in different geographical regions. In the case of the CYP1A2 5′ region, the effects of population substructuring and of natural selection would probably be indistinguishable without information about geographic provenance. These results underscore the continued importance of the sampling of diverse populations across broad geographical regions.