- Research article
- Open Access
A complex selection signature at the human AVPR1B gene
BMC Evolutionary Biologyvolume 9, Article number: 123 (2009)
The vasopressin receptor type 1b (AVPR1B) is mainly expressed by pituitary corticotropes and it mediates the stimulatory effects of AVP on ACTH release; common AVPR1B haplotypes have been involved in mood and anxiety disorders in humans, while rodents lacking a functional receptor gene display behavioral defects and altered stress responses.
Here we have analyzed the two exons of the gene and the data we present suggest that AVPR1B has been subjected to natural selection in humans. In particular, analysis of exon 2 strongly suggests the action of balancing selection in African populations and Europeans: the region displays high nucleotide diversity, an excess of intermediate-frequency alleles, a higher level of within-species diversity compared to interspecific divergence and a genealogy with common haplotypes separated by deep branches. This relatively unambiguous situation coexists with unusual features across exon 1, raising the possibility that a nonsynonymous variant (Gly191Arg) in this region has been subjected to directional selection.
Although the underlying selective pressure(s) remains to be identified, we consider this to be among the first documented examples of a gene involved in mood disorders and subjected to natural selection in humans; this observation might add support to the long-debated idea that depression/low mood might have played an adaptive role during human evolution.
The neurohypophyseal peptide vasopressin (AVP) is involved in different physiological functions, including stimulation of liver glycogenolysis, contraction of vascular smooth muscle cells, antidiuresis and platelet aggregation (reviewed in ). In addition, AVP plays an important role as a regulator of the hypothalamic-pituitary-adrenal (HPA) axis [2, 3]. AVP receptors are G protein-coupled and can be divided in three subtypes: V1a, V1b, and V2, encoded in humans by AVPR1A, AVPR1B and AVPR2, respectively (reviewed in ). The V2 receptor is primarily expressed in the kidney and it controls renal collecting duct water permeability. AVPR1A has wider expression and it regulates physiological effects such as vascular cell contraction, glycogenolysis and platelet aggregation. The type 1b receptor is mainly expressed by pituitary corticotropes and it mediates the stimulatory effects of AVP on ACTH release. Nonetheless, AVPR1B expression has also been described in many brain areas [4, 5] and in different peripheral tissues , while recent evidences have indicated that AVP can induce glucagone and insulin secretion from isolated rodent pancreatic islets through the V1b receptor [6, 7].
Recently, considerable attention has been placed on the role of AVP and its receptors in complex behavioral tracts. Indeed, variations the AVPR1A promoter region have been shown to influence reproductive and social behavior in voles , as well as complex behavioral traits in humans such as altruism , reproductive attitudes [10, 11] and creative dance performance . Therefore, different studies [8, 13] have analyzed the evolutionary history of the type 1a receptor in different mammalian species. In comparison, AVPR1B has attracted less attention, although data from knock-out mice (V1bR-/-) indicate that it plays central roles in both behavioral and metabolic systems. Its regulatory function on the HPA axis is demonstrated by the reduced levels of circulating ACTH and corticosterone under both stress and resting conditions in V1bR-/- animals . These mice also exhibit limited aggressive behavior  and reduced ultrasonic vocalizations in different social contexts . Interestingly, a selective V1b antagonist produces anxiolytic- and antidepressant-like effects in rodents  and in humans AVPR1B variants have been associated with recurrent major depression , early-onset mood disorders  and panic disorder . In line with these findings, the receptor has been proposed as a possible therapeutic target in stress-related disorders .
DNA samples and sequencing
Human genomic DNA was obtained from the Coriell Institute for Medical Research. The genomic DNA of one gorilla and one gibbon was derived from the European Collection of Cell Cultures (ECACC). All analyzed regions were PCR amplified and directly sequenced; primer sequences are available upon request. PCR products were treated with ExoSAP-IT (USB Corporation Cleveland Ohio, USA), directly sequenced on both strands with a Big Dye Terminator sequencing Kit (v3.1 Applied Biosystem) and run on an Applied Biosystems ABI 3130 XL Genetic Analyzer (Applied Biosystem). Sequences were assembled using AutoAssembler version 1.4.0 (Applied Biosystems), and inspected manually by two distinct operators.
Data retrieval and haplotype construction
Genotype data for Yoruba (YRI) and Europeans (EU) were retrieved from the SeattleSNPs website . Genotype data for 238 resequenced human genes were derived from the NIEHS SNPs Program web site . We selected genes that had been resequenced in populations of defined ethnicity including African American (AA), EU, YRI and East Asians (AS) (NIEHS panel 2). In particular, for each NIEHS gene a 2 kb region was randomly selected; the only requirement was that it did not contain any resequencing gap. Haplotypes were inferred using PHASE version 2.1 [23, 24], a program for reconstructing haplotypes from unrelated genotype data through a Bayesian statistical method. Haplotypes for individuals resequenced in this study are available as supplemental material (Additional File 1).
Phylogenetic relationships among primate AVPR1B genes were reconstructed by obtaining a tree with use of MrBayes . In particular, we run a Markov chain for 1 million cycles under the HKY85 model of DNA substitution with no rate variation across sites.
The ratio of dN over dS was calculated using CODEML in the PAML Package (v.3.15) . We used the so-called "free-ratio" model in which dn/ds is free to vary among branches with no variation among sites.
Linkage disequilibrium analyzes were performed using the Haploview software (v. 4.1)  and blocks were identified through the implemented confidence interval algorithm . In particular, marker pairs are defined to be in "strong LD" if the one-sided upper 95% confidence bound on D' is >0.98 and the lower bound is above 0.7; a block is created when if 95% of informative comparisons are in "strong LD".
Tajima's D , Fu and Li's D* and F*  statistics, as well as diversity parameters θW  and π  were calculated using libsequence , a C++ class library providing an object-oriented framework for the analysis of molecular population genetic data. Coalescent simulations were performed using the cosi package  and its best-fit parameters for YRI, AA, EU and AS populations with 10000 iterations. Additional coalescent simulations were computed with the ms software  specifying the number of chromosomes, the mutation parameter estimated from the data, and the recombination rate with 10000 iterations for each demographic model. The other parameters for each model were set as previously proposed [36, 37]. Estimates of the population recombination rate parameter ρ were obtained with the use of the Web application MAXDIP .
The Maximum-likelihood-ratio HKA test was performed using the MLHKA software  using multi-locus data of 15 NIEHS genes (reference genes) and Pan troglodytes (NCBI panTro2) as an outgroup. The 15 reference genes were randomly selected among NIEHS loci shorter than 20 kb that have been resequenced in the 4 populations (panel 2). The reference set was accounted for by the following genes: ENO1, VNN2, MMP12, GLRX, CHRNA4, SULT1C2, PRDX6, H2AFX, ODC1, MT2A, RETN, CYP4B1, RECQL4, MCL1 and MB. In particular, we evaluated the likelihood of the model under two different assumptions: that all loci evolved neutrally and that only the region under analysis was subjected to natural selection; statistical significance was assessed by a likelihood ratio test. We used a chain length (the number of cycles of the Markov chain) of 2 × 105 and, as suggested by the authors , we ran the program several times with different seeds to ensure stability of results. A second multi-locus HKA test was performed using the "HKA" software distributed by Jody Hey  and the same reference loci reported above; 1000 coalescent simulations were performed with the cosi package . In both cases only neutrally evolving sites were considered.
Median-joining networks to infer haplotype genealogy was constructed using NETWORK 4.5 . Estimate of the time to the most common ancestor (TMRCA) was obtained using a phylogeny based approach implemented in NETWORK 4.5 using a mutation rate based on the number of fixed differences between chimpanzee and humans. A second TMRCA estimate derived from application of a maximum-likelihood coalescent method implemented in GENETREE [42, 43]. Again, the mutation rate μ was obtained on the basis of the divergence between human and chimpanzee and under the assumption both that the species separation occurred 6 MYA and of a generation time of 25 years. Using this μ and θ maximum likelihood (θML), we estimated the effective population size parameter (Ne). With these assumptions, the coalescence time, scaled in 2Ne units, was converted into years. For the coalescence process, 106 simulations were performed. A third TMRCA was calculated as previously proposed  that calculates the average nucleotide diversity between the MRCA and each of the chromosomes.
All calculations were performed in the R environment .
The three-dimensional model of the human V1b receptor (Swiss-Prot entry: P47901) was obtained by comparative modelling using the known crystal structure of the closely related bovine rhodopsin (Protein Data Bank ID 1u19) as a template, as provided by MODBASE . The significance of the alignment between target and template sequences is E = 4e-87. The retrieved structure was rendered with FirstGlance in Jmol .
AVPR1Bevolution in primates
As a first step, we wished to gain insight into the evolutionary history of AVPR1B in primates; to this aim the two exons (including the 5' and 3'UTRs) were sequenced from gorilla and gibbon genomic DNA, while the gene sequences for additional primate species (namely, chimpanzee, orangutan and macaque) were retrieved from public databases. A phylogenetic tree (Fig. 1) was produced for the 6 primates with the use of MrBayes . It is worth noting that the failure to resolve the human-chimpanzee-gorilla trichotomy is likely due to the short span of the analyzed region (a total of ~1.7 kb). Using the free ratio model we calculated the dN/dS ratio (ω) along all lineages; in all cases low ratios are observed indicating that AVPR1B has evolved under purifying selection in primates. Calculation of ω for the human-chimpanzee pairwise alignment resulted in a value of 0.28, comparable to the average value for all human genes (ω = 0.23) .
Nucleotide diversity and neutrality tests
We next aimed at analyzing the evolution of AVPR1B in human populations. We therefore resequenced the two exons of AVPR1B and their flanking sequences (including the putative promoter, part of the intron and the 3'UTR) in two populations of Asian and African American ancestry. Additional data referring to YRI and EU subjects were derived from the SeattleSNPs website. The final data set was accounted for by 96 individuals belonging to 4 ethnically distinct populations. A total of 37 SNPs were identified; among these 3 and 1 nonsynonymous substitutions were located in exon 1 (Lys65Asn, Gly191Arg and Ser267Gly) and exon 2 (Arg364His), respectively. Three synonymous coding variants were also identified (Leu130Leu and His224His in exon 1 and Ser373Ser in exon 2).
Analysis of linkage disequilibrium indicated that AVPR1B lies on a single haplotype block in EU and AS, but not in African populations (see Additional File 2). Nucleotide diversity was assessed using two indexes: θW , an estimate of the expected per site heterozygosity, and π  the average number of pairwise sequence nucleotide differences. In order to compare the values we obtained for the two AVPR1B exons, we calculated θW and π for a set of randomly selected 2 kb windows deriving from 238 genes resequenced by the NIEHS program in the same population samples; the percentile rank corresponding to exon 1 (Ex1) and 2 (Ex2) in the distribution of values for reference 2 kb windows is reported in table 1 and indicates that Ex2 displays high nucleotide diversity in all populations, despite showing a level of human-chimpanzee divergence (0.011) comparable to the genome average ; the same holds true for Ex1 when African subjects are considered, while AS and EU show no unusual nucleotide variation pattern in this gene region. Under neutral evolution, values of θW and π are expected to be roughly equal; for both Ex1 and Ex2 this is not verified in most cases (Tab. 1) and we therefore wished to investigate whether AVPR1B might be subjected to natural selection in humans. Widely used neutrality tests include Tajima's D (DT ) and Fu and Li's D* and F* . DT evaluates the departure from neutrality by comparing θW and π. Positive values of DT indicate an excess of intermediate frequency variants and are an hallmark of balancing selection; negative DT values indicate either purifying selection or a high representation of rare variants as a result of a selective sweep. Fu and Li's F* and D* are also based on SNP frequency spectra and differ from DT in that they also take into account whether mutations occur in external or internal branches of a genealogy. Since, population history, in addition to selective processes, is known to affect frequency spectra and all related statistics, we performed coalescent simulations for all populations using a model that incorporates demographic scenarios . Additional demographic models [36, 37] were used for coalescent simulations and the results, which confirm those reported below, are available as additional file 3. Also, in order to disentangle the effects of selection and population history, we exploited the fact that selection acts on a single locus while demography affects the whole genome: as a control data set we therefore calculated test statistics for the 2 kb reference windows deriving from NIEHS genes. Neutrality tests for Ex1 were consistent with neutrality for African populations, while marginally significant negative values of Fu and Li's F* and D* where obtained for EU. In the case of AS only 3 segregating SNPs are observed in the region.
In contrast, data for Ex2 indicate departure from neutrality in most populations (excluding AS) with significantly positive values for most statistics. In line with these findings, DT, as well as Fu and Li's F* and D* calculated for AVPR1B Ex2 rank above the 95th percentile of the distribution of reference 2 kb windows in non-Asian populations. These latter results suggest that nucleotide diversity in Ex2 has been shaped by balancing selection; conversely, the negative statistics observed for EU at Ex1 can in principle be explained by either purifying selection or directional selection since both processes result in an excess of low frequency variants. Fay and Wu's H  is usually applied to distinguish between the two possibilities. Negative H values indicate an excess of high frequency derived alleles, a finding consistent with the action of directional but not purifying selection. Calculation of Fay and Wu's H for EU resulted in a significantly negative value (H = -9.27, p = 0.0006).
A striking feature of these results is the large difference in DT between Ex1 and Ex2 we observe in the EU sample. Such a marked variation in the allele frequency spectrum is even more impressive in light of the strong linkage disequilibrium between the two exons in Europeans (see Additional File 2). In order to evaluate whether such change in the frequency spectrum might be due to chance alone, we performed 10,000 coalescent simulations by generating gene genealogies for a 8.5 kb region (corresponding to the AVPR1B gene); simulations were performed with the estimated recombination rate for AVPR1B and using the cosi package with its best-fit parameters for EU . The simulated samples were then treated as the gene and DT was calculated for the first 2268 bp and the last 2112 bp (corresponding to Ex1 and Ex2). The results indicated the probability of observing a difference in DT values as large as or larger than that we observe for the two AVPR1B exons amounts to 0.035, therefore rejecting a neutral scenario.
Under neutral evolution, the amount of within-species diversity is predicted to correlate with levels of between-species divergence, since both depend on the neutral mutation rate . The HKA test  is commonly used to verify whether this expectation is verified. Here we applied a Maximum-Likelihood-ratio HKA (MLHKA) test  by comparing polymorphisms and divergence levels at AVPR1B Ex1 and Ex2 with 15 NIEHS genes resequenced in the four populations we analyzed (see methods). The results are shown in table 2 and indicate that for Ex1 a reduction in nucleotide diversity versus divergence is detectable in the AS sample although the result is not statistically significant. The opposite situation is verified at Ex2, a significant excess of polymorphisms being observed in all populations. The MLHKA test is relatively robust to demography given its multi-locus comparison framework ; still, while this method is conservative in cases of population expansion (i.e. for populations of African origin), population size bottlenecks might artificially result in significant p values . In order to evaluate whether this is the case, a second multi-locus HKA test was performed using the "HKA" software  which allows estimation of statistical significance through coalescent simulations. These latter were performed using a previously describe demographic model  as above and significant results were obtained for both EU and AS (Tab. 2); in both cases the test of maximum cell value  indicated AVPR1B Ex2 as an outlier (p = 0.018 and 0.001 for EU and AS, respectively).
The signatures of balancing selection are predicted to extend over short genomic distances [54, 55]; nonetheless we wished to verify that the evolution of Ex2 is not influenced by the presence of a linked balanced polymorphism located elsewhere. To this aim, we exploited the availability in SeattleSNPs of resequencing data covering the whole gene and extending 2 kb downstream in the intergenic region. We calculated the ratio of both θW and π to human-chimpanzee divergence in sliding windows for YRI and EU. As shown in figure 2, a peak in both ratios is observed at Ex2 in both populations. This peak decays in the flanking intergenic region.
As far as coding variants are concerned, analysis of SNPs in Ex1 indicated that the derived Gly191 allele has risen to high frequency in most populations and is fixed among AS. With respect to exon 2, a replacement variant (His364 allele) displays higher frequency in EU and AS compared to African populations. Similarity to the other AVP receptors, AVPR1B belongs to the 7-transmembrane receptor family; residue 191 is located in the second extracellular loop of the receptor, a region which is believed to be important for the binding properties of these molecules , including AVPR1A . Indeed, modeling of the protein (Fig. 3) indicated that the 191 residue is located in close proximity to Asp185 and Cys186, two residues in close contact with AVP . With respect to residue 364, it is located in the intracellular C-terminal tail, a protein region involved in molecular interactions with adaptors and scaffolding proteins .
In order to study the genealogy of Ex1 and Ex2 haplotypes we constructed median-joining networks . The topology of Ex1 genealogy indicates that, while African chromosomes are relatively diverged from one another, most Asian and European haplotypes are clustered and all carry the Gly191 allele (Fig. 4A).
With respect to Ex2 genealogy, two major clades separated by long branch lengths are evident (Fig. 4B), each containing common haplotypes. Some population differentiation is also present along clade 2 since most non-African chromosomes are identical and carry the His364 allele, while YRI and AA also display common haplotypes carrying the Arg364 allele.
In order to estimate the TMRCA (Time to the Most Recent Common Ancestor) of the two Ex2 haplotype clades, we applied a phylogeny-based method  based on the measure ρ, the average pairwise difference between the two haplotype clusters. ρ resulted equal to 6.05, so that using a mutation rate based on 35 fixed differences between chimpanzee and humans and a separation time of 6 million years (MY) , we estimated a TMRCA of 3.19 MY years (SD: 673 KY). Given the low recombination rate in the region, we wished to verify this result using GENETREE, which is based on a maximum-likelihood coalescent analysis [42, 43]. The method assumes an infinite-site model without recombination and, therefore, haplotypes and sites that violate these assumptions need to be removed: in this case, only 2 single segregating sites had to be removed. The resulting gene tree, rooted using the chimpanzee sequence, is partitioned into two deep branches (Fig. 5). A maximum-likelihood estimate of θ (θML) of 2.5 was obtained, resulting in an estimated effective population size (Ne) of 13043, a value comparable to most figures reported in the literature . Using this method, the TMRCA of the Ex2 haplotype lineages amounted to 1.94 MY (SD: 497 KY). A third TMRCA of 3.67 MY was estimated by applying a previously described method  that calculates the average sequence divergence separating the MRCA and each of the chromosomes; coalescence time is then obtained by scaling this average divergence to the mutation rate obtained from human-chimpanzee divergence in the region.
All TMRCA estimates indicate an unusually deep coalescent time, as estimates for neutrally evolving autosomal loci range between 0.8 and 1.5 MY . Deep haplotype genealogies might result from both balancing selection and ancient population structure (reviewed in ). Yet, balancing selection is expected to elongate the entire neutral genealogy, while the effects of ancient population structure are reflected in an increase in the genealogical time occupied by single lineages [63, 64]. A possibility to discriminate between these two scenarios is to calculate the percentage of congruent mutations, meaning those that occur on the basal branches of a genealogy . When we applied this approach to the Ex2 genealogy, a percentage of congruent mutations equal to 25% was obtained; this is much lower than previous estimates under a model of ancient population structure, which ranged from 42 to 45% [65, 66], indicating that balancing selection rather than population subdivision is responsible for the maintenance of the two clades.
The interest in the identification and analysis of genomic regions subjected to non-neutral evolution is at least twofold. First, such studies provide insight into the evolutionary history of our species, with special reference to the adaptive events underlying phenotypic changes between humans and primates/hominids. Second, these analysis are expected to result in the identification of functional variants, which in turn, might influence disease susceptibility or drug response. Analysis of AVPR1B evolutionary history in humans might fit both these aims since the gene has been involved in complex behavioral traits in other mammals [14–16] and it has been associated with psychiatric diseases [17–19]. Also, the recently proposed idea  of the receptor as a potential therapeutic target in stress-related disorders suggests that the identification of functional variants will be valuable in the field of pharmacogenetics.
Here we have analyzed the two exons of the gene and, overall, the data we present suggest that AVPR1B has been subjected to natural selection. Analysis of Ex2 strongly suggests the action of balancing selection in African populations and Europeans: the region displays high nucleotide diversity, an excess of intermediate-frequency alleles, a higher level of within-species diversity compared to interspecific divergence and a genealogy with common haplotypes separated by deep branches. Yet, this relatively unambiguous situation is complicated by the presence of unusual features across Ex1. The possibility that a variant in this region has been subjected to directional selection in EU and has reached fixation in AS is supported by few observations: the large variance of DT in EU and the significantly negative H value in this same population. Yet, population structure and different demographic scenarios can yield statistically significant H values, as well . Also, few theoretical models are available to evaluate the difference in DT at different loci; by performing coalescent simulations we verified that such a shift in the allele frequency spectrum is unlikely to be due to chance; consistently, previous descriptions [68, 69] of marked DT differences at linked regions were ascribed to natural selection. Yet, although the closest upstream gene is located more than 72 kb apart, we cannot exclude that variants in Ex1 are hitchhiking with a distant selected allele. Also, the possibility that non-selective forces are responsible for this finding cannot be ruled out and the features we observe might be due to balancing selection acting on exon 2, with Ex1 being neutrally evolving. Although we cannot provide definitive evidence for directional selection at Ex1, it is worth mentioning that the Gly191Arg variant might be a good candidate to represent a functional SNP: as mentioned above, it is located in the second extracellular loop, in close proximity to residues involved in AVP binding (Fig. 3); support to the possible functional significance of this aminoacid replacement, comes from the observation that mutagenesis of most residues in the second extracellular loop of AVPR1A results in altered AVP binding properties ; also different missense substitutions in the same region of AVPR2 have been described to be pathogenetic and cause nephrogenic diabetes insipidus [70–72]. With respect to the possible selection targets in Ex2, SNPs located in the 3'UTR (figure 4B) represent possible candidates, although the splitting of clade 2 into two haplotype clusters with geographically differentiated distribution suggests that additional variants (either the Arg364His SNP or other UTR polymorphisms) might have been subjected to some selective pressure, as well. Of course, specific experiments will be needed to verify whether this prediction is verified. Similarly, a source of selective pressure must be envisaged to explain our findings. Unfortunately not much is known on the function and regulation of AVPR1B in humans; moreover, the gene is not covered by any HapMap SNP at the moment, so that it is likely to be excluded by most genome-wide association studies. Instead, some information about the gene derives from experiments in rodents where AVPR1B function is abolished by either genetic manipulation or antagonist administration; these experiments imply the inherent difficulty in translating mouse phenotypes to humans, which is even more complicated in the case of genes involved in complex behavioral traits, as AVPR1B is suggested to be. Given this premise, different possibilities can be envisaged for selection to act on AVPR1B, including its role in response to stress, psychological/behavioral manifestations and metabolic control. This latter aspect mainly refers to the expression of AVPR1B in human pancreas; experiments on isolated human islets cells, as well as in rodents, indicate that AVP can induce both glucagon and insulin secretion via binding to AVPR1B [6, 73, 74]. Also, V1bR-/-mice display insulin hypersensitivity  possibly due to altered signaling in adipocytes. The possibility therefore exists that one or more variants in AVPR1B have played an adaptive role as thrifty alleles  by favoring energy storage in times of low food availability. Indeed, long-standing balancing selection at the CAPN10 locus, also involved in insulin secretion and action , has previously been demonstrated . In both cases long-term fluctuations in environmental conditions and food availability might be regarded as an explanation for the maintenance of balanced polymorphisms. Nonetheless, our knowledge on the regulation of AVPR1B signaling in peripheral tissues is too limited to allow extensive speculation on this issue. Conversely, the role of AVP in the stress response is better characterized; in response to stressful stimuli, AVP acts synergistically with corticotrophin-releasing hormone (CRH) to facilitate ACTH release, resulting in a consequent increase in plasma corticosterone/cortisol concentration (reviewed in ). Such AVP effect is mediated by AVPR1B  and ACTH measurements in mice have indicated that AVP plays important roles in response to some acute stressors (e.g. hypoglycemia, forced swim stress, lipopolysaccharyde challenge and ethanol intoxication, change in environment), as well as in times of chronic stress [2, 79–83]. Stress responses have an evident adaptive significance and involve both physical and behavioral strategies, with the signaling pathway of AVPR1B being possibly involved in both. Besides the reported associations between AVPR1B and mood or anxiety disorders [17–19], evidences that the receptor has a role in behavioral traits come from the observation that blocking its activity in mice (either by gene knock-out or pharmacologically) results in decreased aggressive behaviors , reduced social interactions  and anxiolytic- and antidepressant-like effects . In both humans and animals (reviewed in ), behavioral responses can be extremely different among individuals of the same species (reviewed in [84, 85]) and recent observations (reviewed in [84, 85]) indicate the maintenance of different behavioral strategies within populations to be adaptive. Frequency-dependent selection is deemed responsible for the coexistence of behavioral types in populations  and, in general, social interaction behaviors are thought to have frequency-dependent payoffs [86–88]. Frequency-dependent selection is among the possible causes for the maintenance of balanced alleles; it is therefore conceivable that a portion of genes involved in behavioral traits have evolved under a balancing selection regime, AVPR1B possibly representing one such example. An alternative and non mutually-exclusive explanation for the maintenance of variability in behavior-related traits is response to changing environmental conditions (reviewed ), which is also expected to result in balancing selection regimes. Clearly, different behavioral strategies can be advantageous or disadvantageous depending on the environment; for example, in unpropitious situations organisms might benefit from down-regulating effort and risk taking, the opposite being true when the environment is more favorable. These explanations have been regarded as possibly supporting an adaptive origin for depression/low mood in humans , a condition which has been associated to polymorphisms in AVPR1B [17, 18]. The evolutionary hypothesis for depression has been highly debated in recent scientific literature (see  for a review) and is based both on the high incidence of this condition in modern populations  and on different observations. For example, anxiety can be advantageous in perceiving danger , a depressive state correlates positively with harm avoidance  and favors disengagement from unproductive efforts . Obviously, these attitudes are selectively advantageous in some environmental conditions (i.e. when action is futile or dangerous) but not in others. Following this line, balancing selection rather than directional or positive selection might be expected at loci which regulate such behavioral traits.
Up to now the evolutionary theory of depression has received few demonstrations. To our knowledge, the one we report here might be the first documented example of a gene involved in mood disorders and subjected to natural selection during human evolutionary history. Nonetheless, it is worth mentioning that the two reports associating AVPR1B variants to depression/mood disorders identified the same major haplotype as being either protective or predisposing in different population of European descent. Dempster et al.  suggested that the discrepancy might result from the use of slightly different sets of markers; we analyzed the haplotype structure of the entire gene in EU (not shown) and found little support for this possibility. Therefore, further studies will be required before a full knowledge can be obtained of AVPR1B evolutionary role as well as of its association with mental diseases.
time to the most common ancestor
Birnbaumer M: Vasopressin receptors. Trends Endocrinol Metab. 2000, 11: 406-410. 10.1016/S1043-2760(00)00304-0.
Tanoue A, Ito S, Honda K, Oshikawa S, Kitagawa Y, Koshimizu TA, Mori T, Tsujimoto G: The vasopressin V1b receptor critically regulates hypothalamic-pituitary-adrenal axis activity under both stress and resting conditions. J Clin Invest. 2004, 113: 302-309.
Volpi S, Rabadan-Diehl C, Aguilera G: Vasopressinergic regulation of the hypothalamic pituitary adrenal axis and stress adaptation. Stress. 2004, 7: 75-83.
Lolait SJ, O'Carroll AM, Mahan LC, Felder CC, Button DC, Young WS, Mezey E, Brownstein MJ: Extrapituitary expression of the rat V1b vasopressin receptor gene. Proc Natl Acad Sci USA. 1995, 92: 6783-6787. 10.1073/pnas.92.15.6783.
Vaccari C, Lolait SJ, Ostrowski NL: Comparative distribution of vasopressin V1b and oxytocin receptor messenger ribonucleic acids in brain. Endocrinology. 1998, 139: 5015-5033. 10.1210/en.139.12.5015.
Oshikawa S, Tanoue A, Koshimizu TA, Kitagawa Y, Tsujimoto G: Vasopressin stimulates insulin release from islet cells through V1b receptors: a combined pharmacological/knockout approach. Mol Pharmacol. 2004, 65 (3): 623-629. 10.1124/mol.65.3.623.
Yibchok-Anun S, Cheng H, Heine PA, Hsu WH: Characterization of receptors mediating AVP- and OT-induced glucagon release from the rat pancreas. Am J Physiol. 1999, 277 (1 Pt 1): E56-62.
Hammock EA, Young LJ: Microsatellite instability generates diversity in brain and sociobehavioral traits. Science. 2005, 308: 1630-1634. 10.1126/science.1111427.
Knafo A, Israel S, Darvasi A, Bachner-Melman R, Uzefovsky F, Cohen L, Feldman E, Lerer E, Laiba E, Raz Y, Nemanov L, Gritsenko I, Dina C, Agam G, Dean B, Bornstein G, Ebstein RP: Individual differences in allocation of funds in the dictator game associated with length of the arginine vasopressin 1a receptor RS3 promoter region and correlation between RS3 length and hippocampal mRNA. Genes Brain Behav. 2008, 7: 266-275. 10.1111/j.1601-183X.2007.00341.x.
Prichard ZM, Mackinnon AJ, Jorm AF, Easteal S: AVPR1A and OXTR polymorphisms are associated with sexual and reproductive behavioral phenotypes in humans. Mutation in brief no. 981. Online. Hum Mutat. 2007, 28: 1150-10.1002/humu.9510.
Walum H, Westberg L, Henningsson S, Neiderhiser JM, Reiss D, Igl W, Ganiban JM, Spotts EL, Pedersen NL, Eriksson E, Lichtenstein P: Genetic variation in the vasopressin receptor 1a gene (AVPR1A) associates with pair-bonding behavior in humans. Proc Natl Acad Sci USA. 2008, 105: 14153-14156. 10.1073/pnas.0803081105.
Bachner-Melman R, Dina C, Zohar AH, Constantini N, Lerer E, Hoch S, Sella S, Nemanov L, Gritsenko I, Lichtenberg P, Granot R, Ebstein RP: AVPR1a and SLC6A4 gene polymorphisms are associated with creative dance performance. PLoS Genet. 2005, 1: e42-10.1371/journal.pgen.0010042.
Fink S, Excoffier L, Heckel G: High variability and non-neutral evolution of the mammalian avpr1a gene. BMC Evol Biol. 2007, 7: 176-10.1186/1471-2148-7-176.
Wersinger SR, Ginns EI, O'Carroll AM, Lolait SJ, Young WS: Vasopressin V1b receptor knockout reduces aggressive behavior in male mice. Mol Psychiatry. 2002, 7: 975-984. 10.1038/sj.mp.4001195.
Scattoni ML, McFarlane HG, Zhodzishsky V, Caldwell HK, Young WS, Ricceri L, Crawley JN: Reduced ultrasonic vocalizations in vasopressin 1b knockout mice. Behav Brain Res. 2008, 187: 371-378. 10.1016/j.bbr.2007.09.034.
Griebel G, Simiand J, Serradeil-Le Gal C, Wagnon J, Pascal M, Scatton B, Maffrand JP, Soubrie P: Anxiolytic- and antidepressant-like effects of the non-peptide vasopressin V1b receptor antagonist, SSR14 suggest an innovative approach for the treatment of stress-related disorders. Proc Natl Acad Sci USA. 9415, 99: 6370-6375. 10.1073/pnas.092012099.
van West D, Del-Favero J, Aulchenko Y, Oswald P, Souery D, Forsgren T, Sluijs S, Bel-Kacem S, Adolfsson R, Mendlewicz J, Van Duijn C, Deboutte D, Van Broeckhoven C, Claes S: A major SNP haplotype of the arginine vasopressin 1B receptor protects against recurrent major depression. Mol Psychiatry. 2004, 9: 287-292. 10.1038/sj.mp.4001420.
Dempster EL, Burcescu I, Wigg K, Kiss E, Baji I, Gadoros J, Tamas Z, Kennedy JL, Vetro A, Kovacs M, Barr CL: Evidence of an association between the vasopressin V1b receptor gene (AVPR1B) and childhood-onset mood disorders. Arch Gen Psychiatry. 2007, 64: 1189-1195. 10.1001/archpsyc.64.10.1189.
Keck ME, Kern N, Erhardt A, Unschuld PG, Ising M, Salyakina D, Müller MB, Knorr CC, Lieb R, Hohoff C, Krakowitzky P, Maier W, Bandelow B, Fritze J, Deckert J, Holsboer F, Müller-Myhsok B, Binder EB: Combined effects of exonic polymorphisms in CRHR1 and AVPR1B genes in a case/control study for panic disorder. Am J Med Genet B Neuropsychiatr Genet. 2008, 147B: 1196-204. 10.1002/ajmg.b.30750.
Griebel G, Simiand J, Stemmelin J, Gal CS, Steinberg R: The vasopressin V1b receptor as a therapeutic target in stress-related disorders. Curr Drug Targets CNS Neurol Disord. 2003, 2: 191-200. 10.2174/1568007033482850.
SeattleSNPs Programs for Genomic Applications. [http://pga.mbt.washington.edu/]
National Institute of Environmental Health Sciences. [http://egp.gs.washington.edu]
Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001, 68: 978-989. 10.1086/319501.
Stephens M, Scheet P: Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am J Hum Genet. 2005, 76: 449-462. 10.1086/428594.
Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.
Yang Z: PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007, 24: 1586-1591. 10.1093/molbev/msm088.
Barrett JC, Fry B, Maller J, Daly MJ: Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005, 21: 263-265. 10.1093/bioinformatics/bth457.
Gabriel SB, Schaffner SF, Nguyen H, Moore JM, Roy J, Blumenstiel B, Higgins J, DeFelice M, Lochner A, Faggart M, Liu-Cordero SN, Rotimi C, Adeyemo A, Cooper R, Ward R, Lander ES, Daly MJ, Altshuler D: The structure of haplotype blocks in the human genome. Science. 2002, 296: 2225-2229. 10.1126/science.1069424.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.
Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133: 693-709.
Watterson GA: On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975, 7: 256-276. 10.1016/0040-5809(75)90020-9.
Nei M, Li WH: Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc Natl Acad Sci USA. 1979, 76: 5269-5273. 10.1073/pnas.76.10.5269.
Thornton K: Libsequence: a C++ class library for evolutionary genetic analysis. Bioinformatics. 2003, 19: 2325-2327. 10.1093/bioinformatics/btg316.
Schaffner SF, Foo C, Gabriel S, Reich D, Daly MJ, Altshuler D: Calibrating a coalescent simulation of human genome sequence variation. Genome Res. 2005, 15: 1576-1583. 10.1101/gr.3709305.
Hudson RR: Generating samples under a Wright-Fisher neutral model of genetic variation. Bioinformatics. 2002, 18: 337-338. 10.1093/bioinformatics/18.2.337.
Marth GT, Czabarka E, Murvai J, Sherry ST: The allele frequency spectrum in genome-wide human variation data reveals signals of differential demographic history in three large world populations. Genetics. 2004, 166: 351-372. 10.1534/genetics.166.1.351.
Voight BF, Adams AM, Frisse LA, Qian Y, Hudson RR, Di Rienzo A: Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes. Proc Natl Acad Sci USA. 2005, 102: 18508-18513. 10.1073/pnas.0507325102.
Hudson RR: Two-locus sampling distributions and their application. Genetics. 2001, 159: 1805-1817.
Wright SI, Charlesworth B: The HKA test revisited: a maximum-likelihood-ratio test of the standard neutral model. Genetics. 2004, 168: 1071-1076. 10.1534/genetics.104.026500.
Hey Lab Home Page. [http://lifesci.rutgers.edu/~heylab/]
Bandelt HJ, Forster P, Rohl A: Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999, 16 (1): 37-48.
Griffiths RC, Tavare S: Sampling theory for neutral alleles in a varying environment. Philos Trans R Soc Lond B Biol Sci. 1994, 344 (1310): 403-410. 10.1098/rstb.1994.0079.
Griffiths RC, Tavare S: Unrooted genealogical tree probabilities in the infinitely-many-sites model. Math Biosci. 1995, 127: 77-98. 10.1016/0025-5564(94)00044-Z.
Evans PD, Gilbert SL, Mekel-Bobrov N, Vallender EJ, Anderson JR, Vaez-Azizi LM, Tishkoff SA, Hudson RR, Lahn BT: Microcephalin, a gene regulating brain size, continues to evolve adaptively in humans. Science. 2005, 309: 1717-1720. 10.1126/science.1113722.
The R project for statistical computing. [http://www.r-project.org]
Pieper U, Eswar N, Davis FP, Braberg H, Madhusudhan MS, Rossi A, Marti-Renom M, Karchin R, Webb BM, Eramian D, Shen MY, Kelly L, Melo F, Sali A: MODBASE: a database of annotated comparative protein structure models and associated resources. Nucleic Acids Res. 2006, D291-295. 10.1093/nar/gkj059. 34 Database
FirstGlance in Jmol Home. [http://www.bioinformatics.org/firstglance]
Bustamante CD, Fledel-Alon A, Williamson S, Nielsen R, Hubisz MT, Glanowski S, Tanenbaum DM, White TJ, Sninsky JJ, Hernandez RD, Civello D, Adams MD, Cargill M, Clark AG: Natural selection on protein-coding genes in the human genome. Nature. 2005, 437: 1153-1157. 10.1038/nature04240.
Chimpanzee Sequencing and Analysis Consortium: Initial sequence of the chimpanzee genome and comparison with the human genome. Nature. 2005, 437: 69-87. 10.1038/nature04072.
Fay JC, Wu CI: Hitchhiking under positive Darwinian selection. Genetics. 2000, 155: 1405-1413.
Kimura M: The Neutral Theory of Molecular Evolution. 1983, Cambridge: Cambridge University Press
Hudson RR, Kreitman M, Aguade M: A test of neutral molecular evolution based on nucleotide data. Genetics. 1987, 116: 153-159.
Wang RL, Hey J: The speciation history of Drosophila pseudoobscura and close relatives: inferences from DNA sequence variation at the period locus. Genetics. 1996, 144: 1113-1126.
Bubb KL, Bovee D, Buckley D, Haugen E, Kibukawa M, Paddock M, Palmieri A, Subramanian S, Zhou Y, Kaul R, Green P, Olson MV: Scan of human genome reveals no new Loci under ancient balancing selection. Genetics. 2006, 173: 2165-2177. 10.1534/genetics.106.055715.
Wiuf C, Zhao K, Innan H, Nordborg M: The probability and chromosomal extent of trans-specific polymorphism. Genetics. 2004, 168: 2363-2372. 10.1534/genetics.104.029488.
Palczewski K, Kumasaka T, Hori T, Behnke CA, Motoshima H, Fox BA, Le Trong I, Teller DC, Okada T, Stenkamp RE, Yamamoto M, Miyano M: Crystal structure of rhodopsin: A G protein-coupled receptor. Science. 2000, 289: 739-745. 10.1126/science.289.5480.739.
Conner M, Hawtin SR, Simms J, Wootten D, Lawson Z, Conner AC, Parslow RA, Wheatley M: Systematic analysis of the entire second extracellular loop of the V(1a) vasopressin receptor: key residues, conserved throughout a G-protein-coupled receptor family, identified. J Biol Chem. 2007, 282: 17405-17412. 10.1074/jbc.M702151200.
Rodrigo J, Pena A, Murat B, Trueba M, Durroux T, Guillon G, Rognan D: Mapping the binding site of arginine vasopressin to V1a and V1b vasopressin receptors. Mol Endocrinol. 2007, 21: 512-523. 10.1210/me.2006-0202.
Hall RA, Lefkowitz RJ: Regulation of G protein-coupled receptor signaling by scaffold proteins. Circ Res. 2002, 91: 672-680. 10.1161/01.RES.0000037000.74258.03.
Glazko GV, Nei M: Estimation of divergence times for major lineages of primate species. Mol Biol Evol. 2003, 20: 424-434. 10.1093/molbev/msg050.
Tishkoff SA, Verrelli BC: Patterns of human genetic diversity: implications for human evolutionary history and disease. Annu Rev Genomics Hum Genet. 2003, 4: 293-340. 10.1146/annurev.genom.4.070802.110226.
Garrigan D, Hammer MF: Reconstructing human origins in the genomic era. Nat Rev Genet. 2006, 7: 669-680. 10.1038/nrg1941.
Takahata N: A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism. Proc Natl Acad Sci USA. 1990, 87: 2419-2423. 10.1073/pnas.87.7.2419.
Wall JD: Detecting ancient admixture in humans using sequence polymorphism data. Genetics. 2000, 154: 1271-1279.
Barreiro LB, Patin E, Neyrolles O, Cann HM, Gicquel B, Quintana-Murci L: The heritage of pathogen pressures and ancient demography in the human innate-immunity CD209/CD209L region. Am J Hum Genet. 2005, 77: 869-886. 10.1086/497613.
Garrigan D, Mobasher Z, Kingan SB, Wilder JA, Hammer MF: Deep haplotype divergence and long-range linkage disequilibrium at xp21.1 provide evidence that humans descend from a structured ancestral population. Genetics. 2005, 170: 1849-1856. 10.1534/genetics.105.041095.
Przeworski M: The signature of positive selection at randomly chosen loci. Genetics. 2002, 160: 1179-1189.
Hamblin MT, Thompson EE, Di Rienzo A: Complex signatures of natural selection at the Duffy blood group locus. Am J Hum Genet. 2002, 70: 369-383. 10.1086/338628.
Smirnova I, Hamblin MT, McBride C, Beutler B, Di Rienzo A: Excess of rare amino acid polymorphisms in the Toll-like receptor 4 in humans. Genetics. 2001, 158: 1657-1664.
Pan Y, Metzenberg A, Das S, Jing B, Gitschier J: Mutations in the V2 vasopressin receptor gene are associated with X-linked nephrogenic diabetes insipidus. Nat Genet. 1992, 2: 103-106. 10.1038/ng1092-103.
Arthus MF, Lonergan M, Crumley MJ, Naumova AK, Morin D, De Marco LA, Kaplan BS, Robertson GL, Sasaki S, Morgan K, Bichet DG, Fujiwara TM: Report of 33 novel AVPR2 mutations and analysis of 117 families with X-linked nephrogenic diabetes insipidus. J Am Soc Nephrol. 2000, 11: 1044-1054.
Ouweland van den AM, Dreesen JC, Verdijk M, Knoers NV, Monnens LA, Rocchi M, van Oost BA: Mutations in the vasopressin type 2 receptor gene (AVPR2) associated with nephrogenic diabetes insipidus. Nat Genet. 1992, 2: 99-102. 10.1038/ng1092-99.
Fujiwara Y, Hiroyama M, Sanbe A, Aoyagi T, Birumachi J, Yamauchi J, Tsujimoto G, Tanoue A: Insulin hypersensitivity in mice lacking the V1b vasopressin receptor. J Physiol. 2007, 584 (Pt 1): 235-244. 10.1113/jphysiol.2007.136481.
Folny V, Raufaste D, Lukovic L, Pouzet B, Rochard P, Pascal M, Serradeil-Le Gal C: Pancreatic vasopressin V1b receptors: characterization in In-R1-G9 cells and localization in human pancreas. Am J Physiol Endocrinol Metab. 2003, 285: E566-576.
Neel JV: Diabetes mellitus: a "thrifty" genotype rendered detrimental by "progress"?. Am J Hum Genet. 1962, 14: 353-362.
Sreenan SK, Zhou YP, Otani K, Hansen PA, Currie KP, Pan CY, Lee JP, Ostrega DM, Pugh W, Horikawa Y, Cox NJ, Hanis CL, Burant CF, Fox AP, Bell GI, Polonsky KS: Calpains play a role in insulin secretion and action. Diabetes. 2001, 50: 2013-2020. 10.2337/diabetes.50.9.2013.
Molen Vander J, Frisse LM, Fullerton SM, Qian Y, Del Bosque-Plata L, Hudson RR, Di Rienzo A: Population genetics of CAPN10 and GPR35: implications for the evolution of type 2 diabetes variants. Am J Hum Genet. 2005, 76: 548-560. 10.1086/428784.
Aguilera G, Rabadan-Diehl C: Regulation of vasopressin V1b receptors in the anterior pituitary gland of the rat. Exp Physiol. 2000, 85 (Spec No): 19S-26S. 10.1111/j.1469-445X.2000.tb00004.x.
Lolait SJ, Stewart LQ, Jessop DS, Young WS, O'Carroll AM: The hypothalamic-pituitary-adrenal axis response to stress in mice lacking functional vasopressin V1b receptors. Endocrinology. 2007, 148: 849-856. 10.1210/en.2006-1309.
Lolait SJ, Stewart LQ, Roper JA, Harrison G, Jessop DS, Young WS, O'Carroll AM: Attenuated stress response to acute lipopolysaccharide challenge and ethanol administration in vasopressin V1b receptor knockout mice. J Neuroendocrinol. 2007, 19: 543-551. 10.1111/j.1365-2826.2007.01560.x.
Keeney A, Jessop DS, Harbuz MS, Marsden CA, Hogg S, Blackburn-Munro RE: Differential effects of acute and chronic social defeat stress on hypothalamic-pituitary-adrenal axis function and hippocampal serotonin release in mice. J Neuroendocrinol. 2006, 18: 330-338. 10.1111/j.1365-2826.2006.01422.x.
Ma XM, Lightman SL: The arginine vasopressin and corticotrophin-releasing hormone gene transcription responses to varied frequencies of repeated stress in rats. J Physiol. 1998, 510 (Pt 2): 605-614. 10.1111/j.1469-7793.1998.605bk.x.
Stewart LQ, Roper JA, Scott Young W, O'Carroll AM, Lolait SJ: The role of the arginine vasopressin Avp1b receptor in the acute neuroendocrine action of antidepressants. Psychoneuroendocrinology. 2008, 33: 405-415. 10.1016/j.psyneuen.2007.12.009.
Sih A, Bell A, Johnson JC: Behavioral syndromes: an ecological and evolutionary overview. Trends Ecol Evol. 2004, 19: 372-378. 10.1016/j.tree.2004.04.009.
Korte SM, Koolhaas JM, Wingfield JC, McEwen BS: The Darwinian concept of stress: benefits of allostasis and costs of allostatic load and the trade-offs in health and disease. Neurosci Biobehav Rev. 2005, 29: 3-38. 10.1016/j.neubiorev.2004.08.009.
Maynard Smith J: Evolution and the Theory of Games. 1982, Cambridge: Cambridge University Press
Wolf M, van Doorn GS, Weissing FJ: Evolutionary emergence of responsive and unresponsive personalities. Proc Natl Acad Sci USA. 2008, 105 (41): 15825-15830. 10.1073/pnas.0805473105.
Dugatkin LA, Reeve HK: Game Theory and Animal Behavior. 2000, Oxford: Oxford University Press
Nesse RM: Is depression an adaptation?. Arch Gen Psychiatry. 2000, 57: 14-20. 10.1001/archpsyc.57.1.14.
Marks IM, Nesse RM: Fear and fitness: an evolutionary analysis of anxiety disorders. Ethol Sociobiol. 1994, 15: 247-267. 10.1016/0162-3095(94)90002-7.
Akiskal KK, Akiskal HS: The theoretical underpinnings of affective temperaments: implications for evolutionary foundations of bipolar disorder and human nature. J Affect Disord. 2005, 85: 231-239. 10.1016/j.jad.2004.08.002.
We are grateful to Dr Roberto Giorda for useful discussion about the manuscript; we also wish to thank Giorgia Menozzi for technical support. MS is part of the Doctorate School of Molecular Medicine, University of Milan.
RC and SR performed all resequencing experiments and analyzed the data; MF performed population genetics analyzes; MS, MF, RC, GPC and UP analyzed and interpreted the data; NB participated in the study coordination; LP performed protein structure analyzes; MC analyzed AVPR1B evolution in primates. MS conceived the study and wrote the paper.