A Homogenizing Process of Selection Has Maintained an “Ultra-Slow” Acetylation NAT2 Variant in Humans

Abstract

N-Acetyltransferase 2 (NAT2) is an important enzyme involved in the metabolism of a wide spectrum of naturally occurring xenobiotics, including therapeutic drugs and common environmental carcinogens. Extensive polymorphism in NAT2 gives rise to a wide interindividual variation in acetylation capacity, which influences individual susceptibility to various drug-induced adverse reactions and cancers. Striking patterns of geographic differentiation have been described for the main slow acetylation variants of the NAT2 gene, suggesting the action of natural selection at this locus. In the present study, we took advantage of whole-genome sequence data available from the 1000 Genomes project to investigate the global patterns of population genetic differentiation at NAT2 and determine whether they are atypical compared with the remaining variation of the genome. The nonsynonymous substitution c.590G>A (rs1799930) defining the slow NAT2*6 haplotype cluster exhibited an unusually low F_ST value compared with the genome average (F_ST = 0.006, P = 0.016). It was indicated as the most likely target of a homogenizing process of selection promoting the same allelic variant in globally distributed populations. The rs1799930 A allele has been associated with the slowest acetylation capacity in vivo, and its substantial correlation with the subsistence strategy adopted by past human populations suggests that it may have conferred a selective advantage in populations shifting from foraging to agricultural and pastoral activities in the Neolithic period. Results of neutrality tests further supported an adaptive evolution of the NAT2 gene through either balancing selection or directional selection acting on multiple standing slow acetylation variants.

The human acetylation polymorphism is one of the oldest and best-characterized pharmacogenetic traits that underlie interindividual and interethnic differences in response to xenobiotics. Acetylation catalyzed by N-acetyltransferase 2 (NAT2) is a major biotransformation pathway for aromatic and heterocyclic amines present in the environment and diet, which can be either detoxified or bioactivated into metabolites that have the potential to cause toxicity and cancer (Butcher et al. 2002; Hein 2002). From the clinical point of view, NAT2 acetylation is increasingly recognized as associated with significant health problems. Many clinically useful drugs are excreted by acetylation, some of them crucial in the treatment of diseases representing worldwide concerns, such as tuberculosis, AIDS-related complex diseases, and hypertension. An individual’s acetylation status has proven to be an important determinant of both the effectiveness of prescribed medications and the development of adverse drug reactions and toxicity during drug treatment (Ladero 2008; Meisel 2002). Moreover, epidemiological studies have associated acetylation phenotype with increased susceptibility to various cancers following exposure to aromatic amine carcinogens (Agúndez 2008; Hein 2006; Sanderson et al. 2007; Selinski et al. 2013).

Figure 1.

Distribution of rs1799930 (c.590G>A) allele frequencies in human populations. Allele frequency data for the rs1799930 SNP were collected for 146 worldwide samples by performing an extensive survey of the literature. Description and numbering of samples and retrieved allele frequency data are provided in Supplemental Table S2. Only those samples comprising at least 20 individuals (N = 105) were represented on the map. Six samples could not be localized on the map because of unspecified sampling location (sample 70) or because of divergence between sampling location and region of origin (samples 42, 55, 60, 61, and 102); these samples are displayed in a box beneath the caption.

N-Acetylation activity has been investigated in a wide range of populations, leading to classification of humans into two main phenotypes: fast acetylators, who exhibit the so-called wild-type or normal acetylation activity, and slow acetylators, characterized by decreased enzyme activity. The proportions of rapid and slow acetylators vary remarkably among populations of different ethnic and/or geographic origin (Sabbagh et al. 2011; Walker et al. 2009; Weber and Hein 1985). Depending on the test substrate administered, a trimodal rather than bimodal distribution can be observed, revealing an additional, intermediate phenotype (Cascorbi et al. 1995; Kilbane et al. 1990; Parkin et al. 1997). Moreover, recent results suggest that the slow acetylator phenotype is not homogeneous and that several slow acetylator phenotypes may exist, resulting from allelic heterogeneity and differential functional effects of the slow acetylation alleles (Ruiz et al. 2012; Selinski et al. 2013). A refinement in phenotype inference, notably by the consideration of an “ultra-slow” acetylator category, is advocated to help identify new clinically relevant associations with one or more of these phenotype subcategories.

Acetylation polymorphism arises from allelic variations in the NAT2 gene, which result in the production of arylamine N-acetyltransferase 2 (NAT2) proteins with variable enzyme activity or stability. The NAT2 gene contains two exons with a relatively long intronic region of about 8.6 kb. Exon 1 is very short (100 bp), and the entire protein-coding region is contained within the 870-bp exon 2. Extensive polymorphism has been described in exon 2, with 38 nucleotide variations registered to date (see the NAT database at http://nat.mbg.duth.gr/). Of these, four common nonsynonymous substitutions at positions 191, 341, 590, and 857 are the most studied and characterize the major NAT2 slow haplotype clusters (NAT2*14, NAT2*5, NAT2*6, and NAT2*7, respectively). Individuals who are homozygous or compound heterozygous for two of these low-activity haplotypes are classified as slow acetylators.

NAT2 acetylation has attracted much research interest in evolutionary biology, and several population genetic studies have attempted to clarify the role that slow acetylation could have played in the adaptation of our species (Fuselli et al. 2007; Luca et al. 2008; Magalon et al. 2008; Mortensen et al. 2011; Patin et al. 2006; Sabbagh and Darlu 2006; Sabbagh et al. 2011). The high prevalence of slow acetylators in humans (well above 50% worldwide) is thought to be a consequence of the shift in modes of subsistence and lifestyle in the last 10,000 years, which triggered significant changes in diet and human exposure to xenobiotic compounds. Several surveys of NAT2 sequence variation have indeed provided compelling evidence that at least some of the slow acetylation variants of NAT2 have been driven to present-day frequencies through the action of natural selection (Luca et al. 2008; Magalon et al. 2008; Mortensen et al. 2011; Patin et al. 2006; Sabbagh et al. 2011). The slow acetylation phenotype may thus have been a key adaptation to increase our species fitness in response to the transition from foraging to farming.

Striking patterns of geographic differentiation have been described for the major NAT2 slow acetylation variants (García-Martín 2008; Sabbagh et al. 2011). The function of NAT2 in mediating the interactions between humans and their chemical environment, which varies depending on diet and lifestyle, makes it an excellent candidate for population-specific selection pressures. Notably, an unusually high level of population differentiation between East Asians and other populations (F_ST values ~ 0.40) has been described for the c.341T>C slow acetylation variant (rs1801280), as well as the two linked c.481C>T (rs1799929) and c.803A>G (rs1208) nonfunctional SNPs, compared with an empirical distribution of F_ST computed across a 400-kb region encompassing the whole human NAT gene family (Sabbagh et al. 2008). In contrast, the slow 590A variant (rs1799930) was found to occur at roughly similar frequencies among widely dispersed populations (Luca et al. 2008; Sabbagh et al. 2011). Figure 1 describes the frequency distribution of this variant in 105 worldwide samples. Such a low level of geographic differentiation may suggest a homogenizing process of natural selection, promoting the same allelic variant in otherwise disparate populations (through either directional or balancing selection). Although many polymorphisms have been described in other regions of the NAT2 gene (Mortensen et al. 2011; Patin et al. 2006), limited data exist on the geographic distribution of these variants in worldwide populations.

In this study, we took advantage of wholegenome sequence data available from the 1000 Genomes (1KG) project to explore global patterns of population genetic differentiation for the whole set of variants occurring in the entire NAT2 gene sequence (~10 kb). An outlier approach was used to determine whether the patterns of geographic differentiation at this locus were atypical compared with those observed for the remaining variation of the genome. Selection tests based on the site frequency spectrum and extended haplotype homo-zygosity (EHH) were further applied to determine the possible role of natural selection in shaping the atypical patterns observed.

Materials and Methods

Data Retrieval

Whole-genome variation data generated by the 1KG project in 1,089 unrelated individuals were directly downloaded from the 1KG website (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/), using the phase 1 integrated release, version 3, dated April 2012 (1000 Genomes Project Consortium et al. 2012). The 1,089 individuals are drawn from 14 different populations in sub-Saharan Africa, Europe, East Asia, and the Americas (see Table 1). From the obtained .vcf (variant call format) files, we extracted exclusively the low-coverage VQSR (variant quality score recalibrator method) single nucleotide variant (SNV) calls to avoid any bias that might result from differences between low-coverage whole-genome calls and high-coverage exome SNV calls. Indels were not used. Functional annotation of the 36,382,866 SNVs retrieved was performed using classification from the dbSNP database, build 137. SNVs were assigned to two main classes: genic and nongenic. Genic SNVs were further classified as intronic, 5’-UTR, 3’-UTR, coding synonymous, coding non-synonymous, or splice site.

Population Genetic Differentiation

Global levels of population genetic differentiation at NAT2 (chr8: 18248755–18258723 in the human GRCh37/hg19 assembly) were evaluated by using the fixation index F_ST (Wright 1951), which quantifies the proportion of genetic variance explained by allele frequency differences among populations. F_ST ranges from 0 (for genetically identical populations) to 1 (for completely differentiated populations). F_ST scores were computed for all NAT2 SNVs occurring with a minor allele frequency (MAF) ≥ 0.05 in at least one of the 14 1KG populations, using the BioPerl module PopGen (Stajich et al. 2002). Extreme values of F_ST can result from natural selection but also from nonselective events linked to the demography of populations, such as genetic drift. Because such nonselective processes randomly act on the genome, they are expected to have the same average effect across the genome, in contrast to natural selection, which impacts population differentiation in a locus-specific manner. The genome-wide variation data provided by the 1KG project can thus be used to infer the action of natural selection by adopting an outlier approach (Kelley et al. 2006). For that purpose, we built eight empirical distributions of the F_ST statistic by considering different subsets of SNVs defined according to their physical location and/or functional impact. To obtain distributions of likely independent observations, a linkage disequilibrium (LD)-based pruning procedure was applied to each of these eight subsets using Plink (Purcell et al. 2007) with default parameters (pruning based on a variance inflation factor of at least 2 within each sliding window of 50 SNVs with a step of five SNVs). This resulted in a total of 25,532,386 independent autosomal SNVs included in the genome-wide empirical distribution. These numbers were 15,141,160 nongenic, 11,282,100 genic, 10,477,050 intronic, 24,395 5’-UTR, 198,718 3’-UTR, 107,644 coding synonymous, 146,572 coding nonsynonymous, and 1,912 splice-site distributions. These eight distributions were then used as reference to assess whether the patterns of genetic differentiation observed at NAT2 are atypical. Empirical P-values were estimated as the proportion of F_ST scores in the empirical distribution that are either higher (diversifying selection) or lower (homogenizing selection) than the value observed at the locus of interest. Since F_ST strongly correlates with heterozygosity (Barreiro et al. 2008; Beaumont and Nichols 1996; Elhaik 2012), empirical P-values were calculated within bins of SNVs grouped according to their global MAF. A total of 27 bins were considered for the whole MAF range: ten bins of size 0.001 for MAF = 0–0.01, nine bins of size 0.01 for MAF = 0.01–0.10, and eight bins of size 0.05 for MAF = 0.10–0.50.

Selection Tests

To determine whether natural selection has played a role in the unusual patterns of geographic differentiation disclosed, we used two complementary approaches based on the allele frequency spectrum of segregating sites and on the local haplotype structure. Tajima’s D (Tajima 1989) is a classical neutrality test that compares estimates of the number of segregating sites and the average number of pairwise differences between nucleotide sequences (π). A zero value of the test statistic D is expected under the null hypothesis of selective neutrality, a positive D indicates balancing selection, and a negative value indicates directional selection. Tajima’s D scores were computed across the whole NAT2 coding region by using a sliding window approach with a window size of 1 kb and a step size of 100 bp. Statistical significance of the test statistic was assessed using an empirical approach. From the genome-wide data available from the 1KG project, we selected a set of unlinked noncoding regions expected to be mostly neutrally evolving. A total of 100 autosomal regions of 1 kb were selected that met the following criteria: (a) at least 100 kb away from any known or predicted genes or expressed sequence tag or region transcribed into mRNA; (b) outside any segmental duplication or region transcribed into a long noncoding RNA or conserved noncoding element, as defined in Woolfe et al. (2007); (c) distant from each other by at least 100 kb and not in LD with each other; and (d) containing a number of SNVs equal to the mean number of SNVs included in the 1-kb sliding windows spanning the NAT2 coding region. Tajima’s D scores were computed for these 100 regions so as to obtain the null (neutral) distribution of the test statistic in each population sample. An empirical P-value was estimated at each sliding window position within NAT2 by considering the proportion of regions showing a test statistic higher (excess of intermediate-frequency variants) or lower (excess of low-frequency variants) than the value observed at that specific position.

We next used methods based on the extended haplotype homozygosity (EHH) measure, that is, the sharing of identical alleles across relatively long distances by most haplotypes in a population sample (Sabeti et al. 2002). We calculated the integrated haplotype score (iHS) (Voight et al. 2006), which compares the rate of EHH decay observed for both the derived and ancestral allele at each core SNV. An extremely positive or negative value at the core SNV provides evidence of positive selection with unusually long haplo-types carrying the ancestral or the derived allele, respectively. The raw iHS scores were computed for all NAT2 SNVs using the iHS option implemented in WHAMM! software (Voight et al. 2006), which we slightly modified to speed up computation times: thresholds for EHH decay were modified from 0.25 to 0.15, and the size of the analyzed region was set to 0.2 Mb instead of 2.5 Mb. Information on ancestral allele state was obtained from a four-way alignment of human, chimpanzee, orangutan, and rhesus macaque species, provided by the 1KG consortium (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/supporting/ancestral_alignments/).

We also applied a cross-population test by computing the cross-population EHH (XP-EHH) statistic (Sabeti et al. 2007), which compares the integrated EHH computed in a test population with that of a reference population. XP-EHH scores were computed using the same EHH decay parameters and window size as for iHS. The Yoruba (YRI) sample was used as a reference for samples outside Africa, and the Utah residents of European ancestry (CEU) as a reference for African samples.

iHS and XP-EHH scores were also computed for all available SNVs in the 100 neutral regions described above, thus providing a reference distribution for each test statistic to estimate empirical P-values. Raw scores of iHS and XP-EHH in NAT2 were standardized in bins of derived allele frequency (step size of 0.05) using the corresponding distribution of each statistic. For both iHS and XP-EHH, we used the 1KG phased data from the integrated phase 1 release, version 3 (April 2012), and genetic distances were obtained from the high-density genetic combined map based on 1KG pilot 1 data.

Estimation of the Age of the Derived Allele at rs1799930

The age of the derived A allele at SNP rs1799930 was estimated using the maximum likelihood method implemented in the Estiage software (Génin et al. 2004). This method was originally developed to estimate the age of the most recent common ancestor of a rare mutation involved in rare diseases using microsatellite data. It is based on the length of ancestral haplotype segments around the mutation shared by mutation carriers. Haplotypes in a 20-kb region centered around rs1799930 were obtained from the phased 1KG data (see Supplemental Table S1). Only SNVs were considered, and their genetic positions were obtained from the genetic maps provided on the Beagle website (http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/). Only one SNV every 0.01 cM was kept in the final analysis with Estiage. Age estimates were derived in each population separately. The mutation rate is a key parameter in Estiage that accounts for the fact that the most proximal marker where a different allele from the ancestral haplotype is observed in an individual might in fact indicate not a recombination but a mutation or a genotyping error. For SNVs, mutation rates are assumed to be very low, on the order of 10^–6, but genotyping errors can be relatively high. Therefore, analyses were performed using mutation rate values between 10^–6 and 10^–4 per marker per individual per generation. Similar age estimations were obtained with the different mutation rates, and results are presented for the lower and upper bound values (Table 1). Population allele frequencies at the different markers were estimated on the haplotypes not carrying the A alleles at rs1799930.

In Silico Prediction of SNV’s Functional Effects

The F-SNP method (Lee and Shatkay 2008; see http://compbio.cs.queensu.ca/F-SNP/) was applied to assess the potential functional effect of SNVs. This integrative scoring method combines assessments from 16 independent computational tools and databases, using a probabilistic framework that takes into account both the certainty of each prediction and the reliability of the different tools, depending on the physical and functional annotation of the specific variant tested. It provides a functional significance (FS) score that quantitatively measures the possible deleterious effect of the tested SNV at the splicing, transcriptional, translational, and posttranslational levels. An FS score of 0.5 is considered the cutoff point for predicting a deleterious effect (Lee and Shatkay 2009).

Results and Discussion

Global patterns of population genetic differentiation were examined in the genomic region spanning the entire NAT2 gene (~10 kb), using the F_ST statistic (Wright 1951) and the sequence variation data provided by the 1KG project (1000 Genomes Project Consortium et al. 2012). P-values were estimated from empirical distributions built from the background genomic variation (see Materials and Methods). We assigned to each genetic variant of the NAT2 gene a main P-value derived from the genome-wide empirical distribution and a subset P-value derived from the distribution including the subset of SNVs having a similar location and/or functional impact as the SNV of interest (i.e., nongenic, intronic, 5’-UTR, 3’-UTR, coding synonymous, coding nonsynonymous, and splice site).

No SNVs in NAT2 exhibited significantly high F_ST values compared with the genomic background, when considering the global differentiation among the 14 worldwide populations from 1KG (all P > 0.05). No significant P-values were observed when contrasting East Asia to the rest of the world (data not shown). Although high F_ST scores were observed for the three SNVs c.341T>C (rs1801280), c.481C>T (rs1799929), and c.803A>G (rs1208) in this specific pairwise comparison (F_ST ≈ 0.30), they could not be considered atypical compared with the rest of the genome (P-values = 0.06–0.09). This contrasts with previous findings indicating an atypical pattern of differentiation for these three variants when considering HapMap data and only the set of variants located within a 400-kb region surrounding the NAT2 gene as a reference distribution (Sabbagh et al. 2008). This difference may be due to the different set of populations surveyed or to the more accurate empirical distribution used to represent the background genomic variation in the present study.

Table 1.

Estimation of the Age of the Derived Allele at rs1799930

In contrast, five NAT2 SNVs exhibited unusually low F_ST values compared with the genome average, with both main and subset P-values below 0.05 (Figure 2). Four of them (rs6984200, rs2087852, rs11996129, and rs1112005) are located in the intronic region of NAT2, whereas the fifth one (rs1799930) is a nonsynonymous substitution defining the NAT2*6 slow haplotype cluster (c.590G>A resulting in R197Q). Note that the four intronic SNVs are in high LD with the rs1799930 variant (r² = 0.80–0.88; Table 2). They all occur at high frequencies in the global human population (within the 0.23–0.26 MAF range). Such low levels of population genetic differentiation suggest that at least one of these polymorphisms may be subject to balancing or species-wide directional selection, the rs1799930 being the most likely target given its gene location and functional impact.

To determine whether another putative candidate in the genomic region surrounding NAT2 might explain the patterns observed, through a significant LD with these variants, we extended the analysis to a 600-kb region centered on the human NAT multigene family on chromosome 8 (Figure 3). All the variants exhibiting an r² value above the 0.10 threshold with the rs1799930 SNV are located within a 56-kb region (chr8: 18229877–18285763 in hg19) in the direct vicinity of the NAT2 gene (Figure 3A), making the involvement of another gene in the region unlikely. Table 2 provides the list of variants that show a significantly low interpopulation F_ST value for either the main or subset P-value and are in moderate to strong LD with the rs1799930 variant (r² > 0.50). They are all either intergenic (located up to ~3 kb upstream or 22 kb downstream of the NAT2 gene) or intronic to NAT2. The prediction

Table 2.

List of Variants in the 600-kb Region Surrounding the NAT2 Gene with Significantly Low F_ST Values and in High Linkage Disequilibrium with rs1799930

Figure 2.

Distribution of –log₁₀ (P-values) of interpopulation F_ST scores across the human NAT2 gene. F_ST scores were computed among the 14 worldwide populations from the 1000 Genomes project. P-values were estimated from the lower tail of the empirical distributions. (A) Main P-value derived from the genome-wide empirical distribution (including 25,532,386 SNVs). (B) Subset P-value derived from the empirical distribution, including the subset of SNVs having a similar location and/or functional impact as the SNV of interest (i.e., intronic, 3’-UTR, coding synonymous, and coding nonsynonymous). The red dotted line indicates the 0.05 significance threshold.

Figure 3.

Distribution of –log₁₀ (P-values) of interpopulation F_ST scores across a 600-kb region centered on the human NAT gene family for those variants in linkage disequilibrium (r² > 0.10) with the rs1799930 polymorphism. (A) Level of linkage disequilibrium, as measured by the r² statistic (Hill and Robertson 1968), with the rs1799930 genetic variant. (B) Main P-value estimated from the lower tail of the genome-wide empirical distribution (including 25,532,386 SNVs). (C) Subset P-value estimated from the lower tail of the empirical distribution including the subset of SNVs having a similar location and/or functional impact as the SNV of interest (i.e., nongenic, intronic, 5’-UTR, 3’-UTR, coding synonymous, coding nonsynonymous, and splice site). The red dotted line indicates the 0.05 significance threshold. SNVs are displayed in different colors according to their r² value with rs1799930, ranging from dark blue (r² = 0.10) to dark red (r² = 1.0). The rs1799930 polymorphism is represented as a black triangle. Coding genes and pseudogenes in the 600-kb region are represented below as dark gray and light gray boxes, respectively. The genomic position (in megabases) on chromosome 8 is indicated on the horizontal axis (human GRCh37/hg19 assembly).

of their functional impact with the F-SNP method revealed high FS scores for some of them (FS = 0.50), denoting a potentially deleterious effect by affecting either the transcriptional or splicing regulation. However, the highest FS score was observed for the rs1799930 polymorphism (0.87), which also displayed the lowest subset P-value (P = 0.016; Table 2). Altogether, these results point to the rs1799930 polymorphism as the most likely target of homogenizing selection in the genomic region surveyed.

Interestingly, recent evidence suggests that the NAT2*6 haplotype cluster (characterized by the rs1799930 A slow acetylation allele) is related with the slowest acetylation capacity in vivo and that the homozygous genotype NAT2*6/*6 thus defines a new category of “ultra-slow” acetylators (Ruiz et al. 2012; Selinski et al. 2013). Ultra-slow acetylators have about 30% lower activities for caffeine metabolism compared with other slow acetylators. This is of the same order of magnitude as the reduction in enzyme activity between rapid and intermediate acetylators (Ruiz et al. 2012). These findings are consistent with a previous study by Cascorbi et al. (1995) that demonstrated a markedly decreased NAT2 activity in vivo in NAT2*6/*6 compared with NAT2*5/*5 genotypes. Indirect evidence is also provided by clinical association studies related to both drug toxicity and cancer risk. Anti-tuberculosis drug-induced hepatotoxicity risk has been shown to be particularly high in carriers of the NAT2*6/*6 genotype (An et al. 2012; Huang et al. 2002; Lee et al. 2010; Leiro-Fernandez et al. 2011; Selinski et al. 2014; Teixeira et al. 2011). Similarly, the ultra-slow genotype, and not the common slow NAT2 genotype, has been significantly associated with an increased risk of urinary bladder cancer (Selinski et al. 2013). Altogether, both in vivo phenotyping studies and clinical reports suggest that the NAT2*6 variant is likely to be associated with a specific acetylation phenotype. This could explain why this particular NAT2 slow acetylation variant, and not another one, may have been a specific target of natural selection.

Although a highly homogeneous distribution of the rs1799930 A allele is observed across worldwide populations (with a global frequency of 0.246), resulting in an unusually low F_ST value (interpopulation F_ST = 0.006), a threefold lower frequency of this allele has been reported in hunter-gatherers (~0.08) compared with agriculturalists and pastoralists (~0.25) in a comprehensive survey of human NAT2 variation comprising 128 population samples classified according to their major subsistence strategy (Sabbagh et al. 2011). This significant difference in the frequency of NAT2*6 alleles (P < 0.0001) was identified as the main genetic cause of the higher prevalence of the slow acetylation phenotype in populations practicing farming and herding compared with those relying mostly on hunting and gathering (46% vs. 22%, respectively) (Sabbagh et al. 2011). Given this marked correlation between the rs1799930 A allele and the subsistence strategy adopted by past populations in the last 10,000 years, it has been suggested that this slow acetylation allele may have conferred a selective advantage in populations shifting from foraging to agricultural and pastoral activities in the Neolithic period. New or more concentrated NAT2 substrates introduced in the chemical environment of food-producing communities have likely promoted a slower acetylation rate in these populations. This hypothesis is further supported by the age estimation of the rs1799930 A allele in the 1KG populations provided by a maximum-likelihood method implemented in the Estiage software (Génin et al. 2004; Table 1). Our estimations showed that this allele started to increase in frequency at similar and recent times in all populations, less than 10,000 years ago (except in Japanese, where estimations are just over 10,000 years ago), thereby supporting global expansion since the emergence of agriculture in the Neolithic. Consequently, the markedly low level of population differentiation observed at the rs1799930 locus may result from the convergent selection of the rs1799930 A allele in agriculturalist and pastoralist populations that are now present in most parts of the world.

Several lines of evidence support the hypothesis that the rs1799930 G>A nonsynonymous substitution (R197Q) has specifically occurred in the human lineage. First, the NAT2 197R residue appears to be highly conserved throughout primate evolution, with 100% of the orthologous NAT2 sequences generated in 19 distinct simian species harboring an arginine (R) at this position (Sabbagh et al. 2013). Second, the 197R position was found to be monomorphic in 103 individuals from six great ape species (Pan troglodytes, Pan paniscus, Gorilla beringei, Gorilla gorilla, Pongo abelii, Pongo pygmaeus) (Prado-Martinez et al. 2013; E. S. Poloni, unpublished observations), as well as in 28 rhesus monkeys (Macaca mulatta) fully sequenced for the NAT2 gene (A. Sabbagh, unpublished observations), making the R197Q polymorphism a specific feature of the human lineage. The hypothesis of a transspecies polymorphism maintained for several million years, through shared balancing selection pressures, therefore seems unlikely.

Assuming that the rs1799930 A allele has conferred a selective advantage to populations shifting from food collection to farming and animal breeding in the Holocene, this could have happened either through directional selection or balancing selection. A gene-dose effect has indeed been described for this variant, with a significant trend toward a slower acetylation capacity in individuals carrying an increasing number of NAT2*6 haplotypes (0, 1, or 2) (Ruiz et al. 2012; Selinski et al. 2013). Therefore, heterozygous individuals for this allele display an intermediate metabolic phenotype that may have been advantageous if one considers the competing needs of both maintaining an efficient detoxification of harmful xenobiotics and avoiding the damaging effects of the putative carcinogens that can be activated through NAT2 acetylation. In an attempt to provide further insights into the evolutionary mechanisms that might have driven and maintained the rs1799930 A allele at high frequencies in most human populations worldwide, we carried out several tests of selective neutrality based on the allele frequency spectrum, Tajima’s D (Tajima 1989), and haplotype structure, iHS (Voight et al. 2006) and XP-EHH (Sabeti et al. 2007). An empirical approach using sequence variation data from 100 unlinked noncoding regions was adopted to assess statistical significance.

Table 3.

Results of Selection Tests for the NAT2 rs1799930 Polymorphism

All iHS and XP-EHH scores computed for the rs1799930 SNV in all individual populations from 1KG were not significant at the 0.05 threshold (Table 3). This precludes a clear signal of positive selection for this variant as the one expected under a “hard sweep model,” which assumes the rapid fixation of a single newly arisen advantageous mutation (Pritchard et al. 2010). In contrast, we found significant Tajima’s D scores in the 1-kb regions encompassing the rs1799930 variant in five population samples: British, Finnish, Tuscans, Utah residents of European ancestry, and Puerto Ricans (P < 0.05; Tables 3, 4). We acknowledge that these results become nonsignificant when a correction for multiple testing is applied, but we also note that the ratio of five significant tests out of 14 is higher than the expected 5% proportion of false positives. Nonsignificant scores, but with P-values getting closer to the 5% threshold, were observed in two additional samples (Colombians, P = 0.07; Luhya, P = 0.08). Furthermore, although nonsignificant scores prevent rejection of the null hypothesis of selective neutrality, it is noteworthy that all populations tested but one (Japanese) gave positive Tajima’s D values, suggesting a trend toward an excess of intermediate-frequency variants compatible with the action of balancing selection. Such consistent results for populations with different demographic pasts make it unlikely that they are due to demography rather than to balancing selection. This is also in agreement with previous findings demonstrating globally positive and significant Tajima’s D values in different continental populations and the absence of any signature of positive selection, as detected by EHH-based tests (Fuselli et al. 2007; Luca et al. 2008; Magalon et al. 2008; Mortensen et al. 2011). A notable exception concerns the c.341C>T slow acetylation variant for which a selective sweep was detected in western and central Eurasian populations (Patin et al. 2006) with the long-range haplotype test (Sabeti et al. 2002). We did not confirm such signature of positive selection at this locus in any of the 14 populations from 1KG (both iHS and XP-EHH scores not significant at the 0.05 level). No significant scores were observed for any of the other slow acetylation variants (c.191G>A, c.341T>C, and c.857G>A; data not shown). Therefore, patterns of diversity at NAT2 seem compatible with either balancing selection or a more complex model of “multiallelic” directional selection where different slow variants of NAT2 may have simultaneously become targets of directional selection, thereby generating an excess of intermediate-frequency alleles. This would explain why our conventional tests of selection based on EHH, more suited to detect classical selective sweeps, failed to detect a signature of positive selection at the rs1799930 locus. The signature of selection at this individual position could have been weakened by the global increase in frequency of other NAT2-altering mutations. Note, however, that contrary to c.191G>A (NAT2*14), c.341T>C (NAT2*5), and c.857G>A (NAT2*7), which mainly cluster in specific continental regions (sub-Saharan Africa, Europe, and Asia, respectively; Sabbagh et al. 2011), the cosmopolitan distribution of the c.590G>A variant (NAT2*6) suggests that it may have been positively selected in globally distributed food-producing communities. Finally, the hypotheses of balancing and directional selection are not mutually exclusive, and multiple modes of selection may have operated at the NAT2 locus on a population-specific basis, as previously suggested (Mortensen et al. 2011).

Table 4.

Tajima’s D Scores in the 1-kb Sliding Windows Spanning the Whole NAT2 Coding Region

Conclusion

We have described an atypical pattern of geographic differentiation for five genetic variants of the NAT2 gene: the functional rs1799930 SNP defining the slow NAT2*6 haplotype series and four intronic SNPs in high LD with it. An extended analysis of a 600-kb region surrounding NAT2 pointed to the rs1799930 polymorphism as the most likely target of a homogenizing process of natural selection promoting the same allelic variant in most human populations, resulting in an unusually low F_ST value (F_ST = 0.006). The rs1799930 A allele has been associated with the slowest acetylation capacity in vivo and is much more frequent in agriculturalists and pastoralists compared with hunter-gatherers, suggesting it may have been positively selected in food-producing communities that are now present in most parts of the world. Neutrality tests based on the allele frequency spectrum revealed a trend toward an excess of intermediate-frequency variants at NAT2, compatible with either balancing selection or a more complex model of multiallelic directional selection. Our findings provide further insights into the functional importance of the rs1799930 polymorphism and the role it may have played in human adaptation to fluctuating xenobiotic environments.