Introduction

With approximately 1.6 billion overweight individuals and at least 400 million clinically obese adults world wide, obesity is one of the world’s greatest health concerns. It is a disease which leads to associated health complications including diabetes, cardiovascular disease, musculoskeletal disorders, and cancer (http://www.who.int/topics/obesity). As a result obesity is the focus of much scientific interest. Among several genes that have been associated with obesity, the melanocortin 4 receptor gene (MC4R) has been intensively studied, as mutations in MC4R account for about 6% of clinical obesity (Farooqi et al. 2003), more than any other single gene. Moreover, homozygous Mc4r knockout mice (−/−) become obese as a result of hyperphagia, while heterozygous knockout mice (±) are intermediate in weight and appetite to that of wild-type mice (Huszar et al. 1997). MC4R is a G-protein coupled receptor (GPCR) that is critical for energy homeostasis. It is involved in the melanocortin system, controlling appetite and energy expenditure, is expressed predominantly in the hypothalamic nucleus paraventricularis of the brain and is activated by melanocyte stimulating hormones, which are processed from proopiomelanocortin by prohormone convertases (Butler 2006). Upon activation, MC4R coupled to the Gs/adenylyl cyclase signaling system leads to increased levels of intracellular cyclic AMP (cAMP), which in turn signals the brain to decrease energy intake, i.e. food consumption (Butler 2006).

Many studies have screened the MC4R locus in thousands of obese probands to identify reduced or loss of function mutations associated with dominantly inherited obesity (Biebermann et al. 2003; Branson et al. 2003; Dubern et al. 2001; Farooqi et al. 2003; Gu et al. 1999; Hinney et al. 2003; Hinney et al. 1999; Kobayashi et al. 2002; Larsen et al. 2005; Lubrano-Berthelier et al. 2003; Miraglia Del Giudice et al. 2002; Santini et al. 2004; Srinivasan et al. 2004; Vaisse et al. 2000; Valli-Jaakola et al. 2004; Yeo et al. 2003). To date about 90 such mutations have been identified and functionality determined by an in vitro cAMP accumulation assay. A quantitative effect of MC4R mutations on body mass index (BMI) has also been observed (Dempfle et al. 2004); individuals who are homozygotes or compound heterozygotes for MC4R mutations (non-synonymous, frameshift or premature stop mutations) are more obese than individuals who are heterozygous for the same allele (Hinney et al. 2006), consistent with observations in Mc4r knockout mice. However, it is not uncommon for heterozygous relatives of obese probands to not exhibit early onset obesity (Dempfle et al. 2004). Thus, there are undoubtedly complex interactions between genotype, genetic background and environment which influence the probability and extent of obesity, given a particular MC4R mutation.

Mutations at the MC4R locus appear to be rare in humans. The mean allele frequency of non-synonymous and synonymous mutations in a sample of 4,068 Germans in the KORA study (Hinney et al. 2006) was 0.021%, exemplifying the rarity of MC4R mutations in a random population sample. However, little is known about variation in MC4R in worldwide populations of humans, or about the molecular evolution of MC4R. Here, we analyze MC4R variation in 1,051 individuals from 51 worldwide populations. To determine if MC4R variation in humans is unusual, we sequenced MC4R from eight chimpanzees. In addition to humans and chimpanzees we also sequenced MC4R from 10 other primate species, and analyzed the molecular evolution of MC4R in a total of 41 vertebrate species, to determine the mode of evolution (purifying, neutral or adaptive) that has been acting on MC4R across 420 million years of vertebrate evolution, as well as to determine if the pattern of evolution has changed on any particular lineage during this time. Given previous observations of human MC4R variation (Hinney et al. 2006) and GPCR conservation (Rompler et al. 2007), we hypothesize that MC4R has been largely subject to purifying selection throughout vertebrate evolution, and that MC4R variants in human populations will be rare, deleterious, de novo mutations.

Materials and methods

Human DNA samples and MC4R mutation screening

The MC4R coding sequence from the Human Genome Diversity Panel (HGDP-CEPH), which consists of 1,051 individuals from 51 human populations worldwide (Cann et al. 2002), was analyzed by denaturing high performance liquid chromatography (dHPLC) as described previously (Hinney et al. 1999, 2003). Mutation screens of two overlapping fragments were performed using dHPLC analysis (Yu et al. 2006) on a Transgenomic WAVE® system (Transgenomic, Cheshire, UK). PCRs were performed with primers amplifying MC4R in two overlapping parts with the primer pairs MC4R-1F/MC4R-1R (F1; 615 bp) and MC4R-2F/MC4R-2R (F2, 622 bp); primer sequences and PCR conditions were as described previously (Hinney et al. 1999). The optimal melting temperatures for separation of homo- and hetero-duplexes were deduced from the melting temperature of the PCR-amplicon using WAVEmaker software, version 4.0 (Transgenomic, Cheshire, UK). Five microliters of unpurified denatured–renatured PCR products were loaded onto a preheated specific column (DNASep column, Transgenomic) and eluted from the column by a linear acetonitrile gradient in a 0.1-m triethylammonium acetate (TEAA) buffer, pH 7, and at a constant flow rate of 0.9 ml/min. The gradient was created by mixing eluents A (0.1 m TEAA) and B (0.1 m TEAA, 25% acetonitrile). For analyses of the MC4R F1-fragment, PCR products were eluted from the column (61°C) by increasing buffer B from 54 to 59% within 30 s and from 59 to 68% within 4.5 min. The oven temperature for analyzing the MC4R F2-fragment was set to both 59 and 60°C, to insure detection of variants. For elution of the MC4R F2 products, buffer B was increased from 54 to 59% within 10 s and from 59 to 68% within 4.5 min. All chromatograms were compared with chromatograms of sequenced wild-type samples; PCR amplicons with more than one peak or differences in peak appearance from the wild-type pattern were sequenced (SeqLab, Göttingen, Germany).

MC4R functional assays

Generation of N-terminal HA-tagged mutant receptors

For functional characterization, wild type and mutant MC4 receptors were amplified from genomic DNA using a forward primer containing the Kosac sequence and an N-terminal hemeagglutinine (HA) tag. Wild-type and mutant MC4Rs were cloned into the pcDps expression vector (Prof. Schöneberg, Leipzig) and sequenced (3100 automated sequencer, Applied Biosystems) to ensure correctness.

Cell culture and transfection

COS-7 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM, Sigma, Deisenhausen) supplemented with 10% fetal calf serum (FCS) and 20 mM Glutamin. Cells were incubated at 37°C in humidified air containing 5% CO2. Transfections were carried out using Metafectene (Biontex, Munich) according to the manufacturer’s protocol. Cyclic AMP assays were performed in 12 well plates (2 × 105 cells/well). Cells were transfected with 0.25 μg DNA/well and 1.25 μl Metafectine/well. Cell surface expression assays were performed in 48 well plates (2 × 105 cells/well) and cells were transfected with 0.25 μg DNA/well and 1 μl Metafectene/well.

Measurement of cAMP accumulation

After 48 h of transfection, cells were labeled with 2 μCi/ml of [3H]adenine (Amersham) and incubated overnight. Cells were washed in serum-free DMEM containing 1 mM 3-isobutyl-1-methylxanthine (Sigma, Deisenhausen) and incubated for 1 h at 37°C in the presence of α-MSH, β-MSH and NDP-α-MSH respectively, in a ten-step series with the highest concentration of 1,000 nM (α- and β-MSH) or 100 nM (NDP-α-MSH). Aspiration of the medium and addition of 1 ml ice cold 5% trichloroacetic acid terminated receptor stimulation. The measurement of intracellular cAMP in the supernatant was carried out as described elsewhere (Salomon et al. 1974).

Measurement of cell surface expression

After 72 h of transfection, cells were washed two times with Dulbecco’s phosphate buffered saline (DPBS, Biochrom, Berlin) and fixed for 30 min in 4% formaldehyde in DPBS followed by two times washing in DPBS. After incubation in blocking buffer (10% FCS supplemented DMEM) for 1 h at 37°C followed by a washing step in DPBS, cells were incubated for 2 h in blocking buffer with 1 μg/ml biotin labeled anti-HA monoclonal antibody (Roche, Mannheim) at 37°C followed by three washes in DPBS and incubation in blocking buffer with 1 μg/ml streptavidin labeled peroxidase (Dianova, Hamburg) at 37°C for 1 h followed by three times washing. The color reaction was carried out with 0.1% H2O2 and 10 μg o-phenylendiamine in 0.1 M citric acid and 0.1 M Na2HPO4 at pH 5.2. The reaction was stopped after 10 min with 1 M Na2SO3 in 1 M HCl. Colorimetry was carried out using an anthos reader 2001 (anthos labtech instruments, Salzburg).

MC4R resequencing in primates

The MC4R coding sequence was resequenced in seven Ethiopians, six Germans, eight central chimpanzees (Ptr, Pan troglodytes) from Gabon (Fischer et al. 2004), and one individual from each of ten other primate species, including: bonobo (Ppa, Pan paniscus), gorilla (Ggo, Gorilla gorilla), orangutan (Ppy, Pongo pygmaeus), crab-eating macaque (Mfa, Macaca fascicularis), baboon (Pha, Papio hamadryas), gray langur (Sen, Semnopithecus entellus), marmoset (Cja, Callithrix jacchus), cotton-top tamarin (Soe, Saguinus oedipus), howler monkey (Aca, Alouatta caraya), and spider monkey (Abe, Ateles belzebuth). The Ethiopian and German samples were described previously (Kayser et al. 2003); the bonobo, gorilla, and orangutan DNA samples were provided by Svante Pääbo of the Max Plank Institute for Evolutionary Anthropology; and the macaque, baboon, langur, marmoset, cotton-top tamarin, howler monkey, and spider monkey DNA samples were donated by Christian Roos of the German Primate Center (DPZ).

PCR reactions were carried out in a RoboCycler 96 Gradient Cycler (Stratagene-Agilent Technologies, La Jolla, CA, USA), consisting of an initial denaturing and Taq activation step at 95°C for 15 min, followed by 35 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 2 min, and a final extension step at 72°C for 10 min. Reactions (25 µl) consisted of 7.5 ng of template DNA, 1× PCR buffer (Applied Biosystems, Foster City, CA, USA), 200 μM dNTPs (Amersham Biosciences, Uppsala, Sweden), 400 nM of each primer (MC4R_Fa: GGATTGGTCAGAAGGAAGC and MC4R_Ra: GATATTCTCAACCAGTACCCT ACAC) (Biotez, Berlin, Germany) and 1.25 U of AmpliTaq Gold DNA polymerase (Applied Biosystems, Foster City, CA, USA). PCR products were subsequently sequenced in both directions in a total volume of 10 μl. PCR primers and four other primers (MC4R_Fb: ACACTTCTCTGCACCTCTGGAAC, MC4R_Rb: CAGCAATCCTCTTAATGTGAAGC, MC4R_Fc: AACATTATGACAGTTAAGCGG GT and MC4R_Rc: GGAGAGACAAAAAGTTGCTCGTA) were used as sequencing primers in reactions consisting of 1 μl of Applied Biosystems Big Dye Terminator v 1.1, 0.75× sequencing buffer, 320 nM primer (Biotez, Berlin, Germany) and ~12 ng of amplified template DNA. Sequencing was carried out in an ABI 3730 DNA Analyzer. Sequence data were analyzed manually with Applied Biosystems SeqScape Software V2.5 (Applied Biosystems, Foster City, CA, USA) and with comparison to a reference sequence (NCBI accession number NT_025028, bases: 5829430..5830428).

Additional coding sequence sources

MC4R Coding sequences (cds) from 29 species [bushbaby (0ga, Otolemur garnettii), tree shrew (Tbe, Tupaia belangeri), mouse (Mmu, Mus musculus), mouse (Msp, Mus spretus), rat (Rno, Rattus norvegicus), guinea pig (Cpo, Cavia porcellus), squirrel (Str, Spermophilus tridecemlineatus), red fox(Vvu, Vulpes vulpes), artic fox (Ala, Alopex lagopus), raccoon dog (Npr, Nyctereutes procyonoides), dog (Cfa, Canis familiaris), cat (Fca, Felis catus), horse (Eca, Equus caballus), boar (Ssc, Sus scrofa), cow (Bta, Bos Taurus), microbat (Mlu, Myotis lucifugus), hedgehog (Eeu, Erinaceus europaeus), elephant (Ema, Elephas maximus), opossum (Mdo, Monodelphis domestica), platypus (Oan, Ornithorhynchus anatinus), chicken (Gga, Gallus gallus), goose (Aan, Anser anser), clawed frog (Xtr, Xenopus tropicalis), fugu (Tru, Takifugu rubripes), stickleback (Gac, Gasterosteus aculeatus), medaka (Ola, Oryzias latipes), zebrafish (Dre, Danio rerio), goldfish (Cau, Carassius auratus)] were downloaded from ENSEMBL (http://www.ensembl.org) and the NCBI core nucleotide (http://www.ncbi.nlm.nih.gov) databases. Species and accession numbers can also be found in Supplemental Table 1.

Data analysis

Coding sequences were converted into amino acid sequences with MEGA 3.1 (Kumar et al. 2004), aligned with MUSCLE 3.6 (Edgar 2004), exported in CLUSTAL W format, and converted into a cds alignment with a custom PERL script. The alignment is provided in Supplemental Table 3. A species tree (Benton and Donoghue 2007; Bininda-Emonds et al. 2007) was generated and the transition/transversion bias (κ) was estimated using the CODEML program in the PAML 4a package (Yang 2007). The PHYLIP package (Felsenstein 2004) was used to construct nucleotide and amino acid neighbor-joining (NJ) trees. The program dnadist was used to estimate nucleotide distances, using the estimated κ (2.58) and both the F84 and Kimura substitution models. Amino acid distances were estimated with the program protdist using both the JTT and PAM substitution matrixes, and assuming a gamma distribution for the rate variation. Finally, the program neighbor was used to generate NJ trees. Sequences were randomized on input. The species tree and constructed nucleotide and amino acid trees were analyzed by the K–H and S–H test (Kishino and Hasegawa 1989; Shimodaira and Hasegawa 1999), as implemented in CODEML of the PAML 4a package, to determine if any one tree significantly fit the data better than alternative trees.

The inferred coding sequence for each individual from the HGDP-CEPH panel, as well as for each resequenced Ethiopian, German, chimpanzee and other primate, were analyzed in DNAsp 4.0 (Rozas et al. 2003). Summary statistics for 319 genes, including sequence length, number of segregating sites, nucleotide diversity, and Tajima’s D, were downloaded from the SeattleSNPs database (http://pga.mbt.washington.edu/) in May 2008, and used as a comparative empirical dataset. The estimated number of non-synonymous and synonymous sites for the genes in the SeattleSNPs database was taken from previously published estimates (Bustamante et al. 2005); such estimates were available for 268 genes. Neutral simulations were performed with Hudson’s ms (Hudson 2002). For the human data, simulations were carried out using demographic parameters found to provide the best fit, based on empirical data from three human populations: Africans, Europeans, and East Asians (Schaffner et al. 2005). The model used allows for migration between each population, includes a bottleneck for non-African populations, growth in Africans, and a post-agricultural expansion in each population. For chimpanzee data, simulations included a current effective population size of 118,000 (N2), an ancestral population size of 16,000 (N1), and growth from 500,000 years ago to today from N1 to N2 (Caswell et al. 2008). Simulations were based on the observed number of segregating sites (S), and for humans on a recombination rate of 0.00103/base (Ptak et al. 2004) across the entire sequence length of 999 bp. The significance of observed summary statistics was then determined by comparing empirical observations to 10,000 simulated loci.

Amino acid sequence evolution

The role of selection across the MC4R phylogeny was investigated by analyzing patterns of evolution at the codon level. The non-synonymous to synonymous substitution rate ratio (ω) was estimated by maximum likelihood methods (Goldman and Yang 1994) using the CODEML programs of the PAML4a package (Yang 2007). Three different codon substitution models were used in each analysis to insure the robustness of the results: the F3X4 codon substitution model (Goldman and Yang 1994) derives equilibrium codon frequencies from each of the three codon positions using three nucleotides, and thus has 9 free parameters; the F61 substitution model (Goldman and Yang 1994) uses each codon as a free parameter, estimating codon frequencies from the empirical distribution; and the recently developed FMutSel codon substitution model (Yang and Nielsen 2008) specifically models mutational biases and selection at silent sites by introducing a parameter of codon fitness. This latter model was developed to address possible selection at silent sites resulting in codon usage bias. Furthermore, several models were fit to the data which allow ω to vary across sites, across branches (Yang 1998; Yang and Nielsen 1998), and across clades (Bielawski and Yang 2004). These models and analyses allow inferences concerning the type of selection (neutral, negative, or positive) which has acted across a phylogeny, as well on specific branches and/or sites. A likelihood ratio test (LRT = 2(l 1 − l 0), where l 1 and l 0 are the log likelihoods of the two models) was used to determine if one model fit the data significantly better than another model.

Results

HGDP-CEPH screening

The HGDP-CEPH panel was screened for polymorphisms at the MC4R locus by dHPLC (WAVE) melting curve technology, and aberrant PCR products were re-sequenced to identify polymorphisms. Fourteen segregating sites (S) were found across 1,074 bps in 2,030 chromosomes, with an average variant allele frequency of only 0.17% (Table 1). The most frequent polymorphism (V103I) had a variant allele frequency of just 1.33%. Eleven of the segregating sites are non-synonymous mutations, two are synonymous mutations and one is found in the 5′ region. The mean nucleotide diversity estimate (π) for MC4R is 0.00005, which is lower than the π value for any of the 319 genes in the SeattleSNPs database (mean π = 0.0009, minimum π = 0.00016). The π value for MC4R is also significantly lower than expected (Pvalue = 0.001) based on simulations that incorporate the best-fitting demographic parameters (Schaffner et al. 2005) for three human populations (Africans, Europeans, and Asians). Since the HGDP-CEPH includes more population structure than would be observed in just these three human populations (Li et al. 2008; Rosenberg et al. 2002, 2005), the comparison of π values from the HGDP-CEPH to the empirical SeattleSNPs data and to the simulations based on three populations is conservative, making the extremely low level of diversity at MC4R all the more unusual. In addition to the dHPLC screening of the HGDP-CEPH panel, the MC4R locus was amplified and directly sequenced in six Germans and seven Ethiopians. No mutations were identified in these 26 chromosomes, consistent with the low levels of diversity found in the HGDP-CEPH screening. Thus the low level of diversity found in this worldwide sample is not an artifact of the dHPLC technology.

Table 1 CEPH panel mutations
Full size table

This extremely low level of diversity could reflect either very recent directional selection in all human populations, or continuous purifying selection. Under neutral expectations, the rate of mutations at non-synonymous and synonymous sites are expected to be equal, i.e. the intra-species ratio of the non-synonymous to synonymous polymorphism rate per codon site should equal 1 (pN/pS = ~1). The observed pN/pS value for MC4R is 1.98, but this is not significantly greater than 1 (χ 2 = 0.826, P = 0.362). However, the pN/pS ratio for MC4R is significantly larger (P = 0.008) than pN/pS estimates for 268 genes from the SeattleSNPs database (SeattleSNPs), and moreover it is in the top 1% of estimated values from the Celera exon re-sequencing dataset (Bustamante et al. 2005), when considering genes with informative polymorphic sites at both non-synonymous and synonymous positions. The excess of mutations at non-synonymous sites would seem to suggest that directional selection (or relaxed constraints) may be acting on MC4R, resulting in the accumulation of non-synonymous mutations that may alter protein function. However, a caveat to this observation is that no polymorphism in the CEPH data is above 1.4% in frequency. It is common for low-frequency alleles to be removed from analyses of polymorphic synonymous and non-synonymous mutations, as slightly deleterious mutations segregate at such frequencies and bias estimates of intra-population variation and tests of selection, including the MK test (Charlesworth 1994; Eyre-Walker 2006; Fay et al. 2001). However, removing low-frequency alleles in the MC4R dataset would remove all of the variants. Eight of the observed polymorphisms are singletons, three are doubletons, two are observed in four chromosomes, and V103I, which has been dubbed a protective allele (Heid et al. 2008), is the most common non-synonymous allele observed in 27 out of 2,030 chromosomes (frequency = 1.33%). Thus all non-synonymous mutations are rare and are likely to be transient, deleterious mutations which will ultimately be removed from the population (Boyko et al. 2008). In regard to synonymous polymorphisms, a possible explanation for their rarity in MC4R, beyond the persistent problem of background selection, is codon usage bias resulting in selection on silent sites; we explore this possibility in more detail below.

Functional characterization of non-synonymous mutations

Eight of the eleven observed non-synonymous mutations in the HGDP-CEPH screening were tested in cAMP accumulation and cell surface accumulation assays to measure the effect of each mutation on MC4R function as compared to the wild-type (Table 2); the remaining three sites (R7C, I251L and V103I) have been previously analyzed (Srinivasan et al. 2004; Stutzmann et al. 2007; Xiang et al. 2006). Only one of the 11 non-synonymous mutations (N62S) exhibited a reduced function, as indicated by a strongly reduced cell surface expression, and by signal transduction properties that did not result in an enhancement of intracellular cAMP levels (Table 2). Some mutants (C40Y, I170V and S191T) displayed a reduced cell surface expression; however, maximal stimulation and other parameters with all three tested ligands were comparable to the wild-type receptor, and therefore we define these mutants as having wild-type function.

Table 2 MC4R functional assays
Full size table

Pan troglodytes re-sequencing

To determine if the unusually low π and unusually high pN/pS values are restricted to humans, we sequenced the MC4R coding region in eight central chimpanzees. Four segregating sites were observed (2 synonymous and 2 non-synonymous mutations), from which π was estimated to be 0.00111. This observed π value is not unusual (P = 0.23) when compared to 26 non-coding regions which were sequenced in the same chimpanzees (Fischer et al. 2006), nor is it unusual (P = 0.185) when compared to simulated data based on a demographic model appropriate for chimpanzees (see “Materials and methods”). Thus, unlike humans, the amount of diversity at the MC4R locus is not unusual in central chimpanzees. Further, the estimated pN/pS value for MC4R in chimpanzees is 0.36, which is smaller than but not significantly different from that of humans (Table 3).

Table 3 Polymorphic non-synonymous versus synonymous mutations in homo and pan
Full size table

Primate resequencing

The coding sequence of MC4R was determined from a total of 13 primate species, including humans (Fig. 1, Supplemental Table 1). Published data were added for another primate species, the galago (Otolemur garnettii, Oga) as well as for the tree shrew (Tupaia belangeri, Tbe) a species basal to the primate phylogeny. The conservation of the coding sequence is high across primates. The largest p-distance (0.088; proportion of nucleotide differences) is found between Papio hamadryas and Otolemur garnettii. Only a few amino acid sequence differences (23 in 334 amino acid sites) are found among primates (Fig. 1), across approximately 60 million years of evolution.

Fig. 1
figure 1

Amino acid sequence alignment of all primates used in this study as well as the primate outgroup (Tbe, tree shrew). Variable positions are shaded in gray. Species names corresponding to each three letter code can be found in Supplemental Table 1

Full size image

Phylogenetic analysis by maximum likelihood (PAML)

The MC4R coding sequences for 29 species, covering all of vertebrate evolution, were retrieved from available databases (accession numbers can be found in Supplemental Table 1). Nucleotide and protein sequence distance matrixes were computed using the F84 and Kimura nucleotide substitution models and the JTT and PAM protein substitution models, from which neighbor-joining trees were generated. A “species tree” was created according to the accepted taxonomy (Benton and Donoghue 2007; Bininda-Emonds et al. 2007). To determine if any one tree statistically fits the data better than any other, we used the likelihood-based S–H and K–H tests (Goldman et al. 2000; Kishino and Hasegawa 1989; Shimodaira and Hasegawa 1999), as implemented in the PAML package. The species tree had the largest log likelihood, but no tree was statistically better than any other. Given this result we chose to conduct all further analyses with the species tree (Fig. 2).

Fig. 2
figure 2

MC4R species tree. Omega estimates for each lineage, from the free-ratio (M1) model, across the unrooted MC4R phylogeny. Estimates equal to infinity (Inf) are designated as such when the estimation of dS was equal to 0.0, i.e. no silent substitutions were observed. Foreground and background labels are for the MC clade model analysis of divergent selection pressures; all branches descended from the branch highlighted in red are foreground branches

Full size image

The goal of these analyses is to investigate various ways of estimating ω (the ratio of the rates of non-synonymous to synonymous substitutions per codon site) across the phylogeny to determine if the MC4R locus has been evolving neutrally or rather has been subject to functional constraints, and if the latter is the case, to then identify essential residues which are important to protein function. Each analysis was performed using one of three different codon substitution models: F3X4, F61 and FMutSel (see “Materials and methods” for details). The results are largely concordant between substitution models and thus results will only be presented for the F3X4 model, unless otherwise specified. Further, because alignment gaps are ignored in the analyses, each analysis was conducted with both the full data set of 41 species and with a reduced data set which includes the primate species and the basal tree shrew (Tbe) lineage. The full data set provides more power to detect changes in evolutionary rates, whereas the reduced data set permits the evaluation of evolutionary rates at sites that are present in all primate species but missing in one or more of the other vertebrate species. As the results from both the full dataset and the reduced dataset of only primate sequences were largely concordant, we present and discuss the results based on the full dataset, and present those for the reduced dataset only when they differ.

First, a global ω value was estimated for the entire tree, using the M0 model, which constrains each lineage in the tree to evolve with the same ω value. The global ω value for the MC4R phylogeny is 0.043, and is a significantly better fit to the data (P < 0.001) than when the ω value is set to 1. This ω estimate suggests strong negative selection acting throughout the MC4R phylogeny.

We next allowed the ω value to vary from lineage to lineage by fitting the free-ratio model, which estimates the ω value for each lineage. A likelihood ratio test (LRT, two times the difference in log likelihoods of the two models) indicates that this model provides a significantly better fit to the data, compared to the global M0 model (LRT = 162, P < 0.001, df = 78), although this is not true for the primate phylogeny alone (LRT = 28.5, P = 0.39, df = 27). This difference between the full and reduced dataset may reflect a change in the evolutionary rate of MC4R between primates and other vertebrates, a possibility we investigate below. Overall, even as negative selection predominates across the phylogeny there is significant variation in ω values (i.e. levels of constraint) across lineages; the lineage-specific ω values are presented in Fig. 2.

We next investigated variation in ω values across sites, by fitting the M3(discrete) sites model to the MC4R phylogeny. This model defines three categories of sites (0, 1 and 2) each with a free and independent estimate of ω. Site class 0 has an ω value (ω0) of 0.009 with 71% (p 0) of the sites designated to this category; site class 1 has p 1 = 0.24 and ω1 = 0.094; and site class 2 has p 2 = 0.05 and ω2 = 0.483. This model is a better fit to the data than the M0(one-ratio) model (LRT = 536, P < 0.001, df = 4). Thus, there is variation in the rate of evolution at different sites across the MC4R phylogeny, with most sites (71%) experiencing very low rates of non-synonymous substitutions (ω = 0.009). There is also significant evidence of site variation in evolutionary constraint when using the primate phylogeny alone (LRT = 15, P < 0.001, df = 2). However, here the M3(discrete) model with three categories of sites is not a better fit than the same model with two categories of sites. As such, when analyzing the primate phylogeny alone 93% of sites (p 0) have a ω0 estimate of 0.026, while the remaining 7% of sites (p 1) have a ω1 estimate of 0.63.

Moreover, the three site classes (0, 1, and 2, corresponding to sites under strong, moderate, and weak purifying selection) are not evenly distributed over the topology of the MC4R protein (Fig. 3). The occurrence of each site class in the three topological categories (extracellular, transmembrane, and intracellular; Table 4) shows a significant deviation from random expectations (χ 2 = 71.22, P < 0.001). The extracellular region of MC4R has a deficiency of class 0 sites (very strong purifying selection), and an excess of class 2 sites (those with the largest ω value). This observation, along with the fact that the extracellular N-terminus has many alignment gaps (Fig. 3) and displays the lowest degree of evolutionary conservation, suggests that this section of the protein can tolerate more mutations, relative to other regions of MC4R, while still maintaining protein function. The analogous analysis with the primate phylogeny alone shows a similar, although non-significant (χ 2 = 4.14, P = 0.11) pattern (Supplemental Fig. 1).

Fig. 3
figure 3

Bar plot. Site class probability bar plot for each codon across the MC4R sequence. White spaces are alignment gaps

Full size image
Table 4 Site and protein topology matrices
Full size table

To determine if these site classifications are informative in identifying functionally critical residues, we compared the association of MC4R non-synonymous mutations found in obese cohorts as well as in population based samples, including the CEPH panel, to the site classifications of the M3 analysis. If the level of evolutionary constraint is informative in predicting the deleterious effect of a mutation then we would predict that mutations occurring at class 0 sites (i.e. those under very strong purifying selection) would be more likely to reduce MC4R function in cAMP assays. Correspondingly, we would predict that mutations that function like the wild type (LWT) MC4R in cAMP assays should occur more frequently at class 1 and 2 sites (i.e. those under less evolutionary constraint). Using our own data (Table 2) and published data we identified 42 codons (Supplemental Table 2) for which a non-synonymous mutation was observed, functional tests were performed, and the site classification could be determined (Farooqi et al. 2000; Gu et al. 1999; Hinney et al. 1999, 2003, 2006; Kobayashi et al. 2002; Larsen et al. 2005; Nijenhuis et al. 2003; Santini et al. 2004; Srinivasan et al. 2004; Stutzmann et al. 2007; Vaisse et al. 2000; Valli-Jaakola et al. 2004; Xiang et al. 2006; Yeo et al. 2003). Of these 42 non-synonymous substitutions, 23 caused a reduction in function (RF), and 22 (96%) of them occur at a class 0 site, under strong evolutionary constraint (Supplemental Table 2). Of the 19 non-synonymous substitutions which retain LWT function, nine (47%) occurred at a class 0 site. Overall, there is a significant association of class 0 sites with mutations which cause a reduction in function (Fisher’s exact test, P < 0.001). Hence, non-synonymous mutations that have reduced function are significantly more likely to occur in codons under “very strong” purifying selection, while non-synonymous mutations that retain LWT functionality are over-represented in classes experiencing less evolutionary constraint (site class 1 and 2).

Finally, we investigated the possibility that certain codons in the MC4R gene have experienced recurrent positive selection across the MC4R phylogeny, possibly providing an adaptive function for particular lineages. We did this by fitting two neutral models of codon evolution (M1a and M7) which allow ω values to range from 0 to 1 across codon sites, and two positive selection models (M2a and M8) which allow a proportion of the sites to have ω values larger than 1 (Nielsen and Yang 1998; Wong et al. 2004; Yang et al. 2005). The M1a and its corresponding M2a model are discrete models, where ω estimates take a discrete value, while the M7 and M8 models estimate ω based on a continuous beta distribution. No selection model was significantly better than the corresponding neutral model, indicating that there is no evidence for recurrent positive selection at particular amino acid residues across the MC4R phylogeny.

Divergent selection pressures in primates

The above analyses indicate that overall, the MC4R protein has been under strong negative selection throughout vertebrate evolution. Yet, is there any evidence for a change in selection pressure on MC4R during primate evolution? To investigate this, a clade model of codon evolution (MC) was fit to the data to determine if there are divergent selection pressures between the primate clade and all other lineages in the MC4R phylogeny (Bielawski and Yang 2004). This model has 3 site categories: ω at class 0 sites can range from 0 to 1 (purifying selection); ω at class 1 sites are set to 1 (neutral evolution); and class 2 sites are under divergent selection pressure where ω is freely estimated. Further, class 2 sites have two estimations of ω: one estimate (ω2) is the ω estimate for the background branches, i.e. all branches not included in the primate clade; the other estimate (ω3) is the ω estimate for the primate clade. Thus, it is the class 2 sites, which freely estimate ω for the foreground branch (primate lineage) and the background branch (all other lineages) independently, that allow the determination of possible divergent selection pressures acting on these two opposing branches of the MC4R phylogeny (Fig. 2). To construct an LR-test this clade model is compared to the neutral sites model (M1a) with 3 degrees of freedom. The results of this model indicate that 22% of the sites are designated to site class 0 (purifying selection) where ω0 is equal to 0.134, 3% of the sites belong to site class 1 (neutrally evolving), and 75% of the sites are under divergent selection pressure where ω2 (background) equals 0.011 and ω3 (primate clade) equals 0.003. This model is significantly better than the M1a model (LRT = 350, P < 0.001), and thus estimates that 75% of the sites across the MC4R locus are evolving in the primate clade under a level of purifying selection that is three times stronger than in the rest of the phylogeny. While MC4R has been subject to strong purifying selection throughout vertebrate evolution, it appears that the level of evolutionary constraint has increased significantly during primate evolution.

Selection on codon usage bias

Recently, a method was developed to address selection on codon usage bias, which would influence variation at silent sites (Yang and Nielsen 2008). This model, FMutSel, explicitly models mutational biases and selection at silent sites in comparison to a null model (FMutSel0) which only models mutational biases. The model not only determines if there is selection at silent sites creating codon usage biases, but also estimates the proportion of all possible mutations that would be advantageous and the mean selection coefficient for advantageous and disadvantageous mutations. With respect to the MC4R phylogeny, the FMutSel model is a significantly better fit to the data than the null model (LRT = 111, P < 0.001), suggesting codon usage biases across the MC4R phylogeny. Furthermore, the model estimates that 80% of all observed substitutions at silent sites are advantageous. The mean selection coefficient estimates across the MC4R locus are 1.98 for positively selected sites, and −2.33 for negatively selected sites; these values are approximately five and four times larger, respectively, than the average values estimated for 5,600 genes (Yang and Nielsen 2008). Thus, mutations which occur at silent sites at the MC4R locus are either largely advantageous or largely disadvantageous, and hence would be expected to be quickly fixed or lost via selection.

Discussion

The MC4R locus is highly conserved and exhibits evidence of strong purifying selection across the entire vertebrate phylogeny, as evident by the very low global ω estimate of 0.043. Furthermore, there is no evidence of any one lineage experiencing an acceleration of protein evolution that may be adaptive (Fig. 2). The level of functional constraint appears high across the entire MC4R phylogeny, which is consistent with the MC4R locus, and GPCRs in general, being vital to the organisms carrying them (Rompler et al. 2007; Schoneberg et al. 2007).

Despite the high level of purifying selection across the entire MC4R phylogeny, the amount of functional constraint has apparently increased even more during primate evolution. Approximately 75% of all codons exhibit three times more evolutionary constraint in the primate lineage than all other lineages, which would suggest that the physiological role of MC4R is more constrained in primates than in other vertebrates. This may be a consequence of changes in MC4R regulation and its effect on animal behavior. For example, fish deprived of food increase MC4R expression in the liver (Kobayashi et al. 2008); only mammals express MC4R exclusively in the brain where there is an apparent gene–environment interaction influencing eating behavior (Butler and Cone 2003). It would, however, be necessary to determine if MC4R regulation and function vary between primates and other mammals to fully understand the observation of increased functional constraint in primates.

Against this background of even higher purifying selection on MC4R in primates, observations from human polymorphism data are consistent with a continuing increase in evolutionary constraints on the human lineage. The MC4R locus has a significantly low level of diversity in humans, when compared to other empirical data (Bustamante et al. 2005; SeattleSNPs) and to simulations. This observation appears to be restricted to humans, as MC4R diversity in chimpanzees is not unusually low as compared either to simulations or to empirical data from non-coding regions (Fischer et al. 2006).

Humans do exhibit an unusually large pN/pS ratio, as compared to empirical data, which would suggest relaxation of functional constraints. This observation, at first, seemed at odds with the overall low level of diversity at MC4R. However, the observed level of variation at synonymous and non-synonymous sites is not statistically different from that observed in chimpanzees, although this conclusion may very well be the result of so few variable sites, and hence a corresponding lack of power in the statistical test. More importantly, all observed human mutations are rare (mean = 0.17%) suggesting that most of these mutations are transient deleterious mutations on their way to being removed from the population, consistent with strong purifying selection on MC4R and the functional data of the observed non-synonymous mutations. This, in turn, predicts that the polymorphisms observed in the human data are all new mutations, and hence should reflect the presumably random mutational spectrum at non-synonymous and synonymous positions, driving the unusually large pN/pS. To test this prediction, we performed a very simple simulation: we took the full length of the MC4R coding sequence (996 bp, with an estimated 732 non-synonymous positions and 264 synonymous positions), randomly mutated 13 positions, computed the pN/pS value, and then repeated this process 2,000 times. The probability, based on this empirical distribution, of observing a pN/pS value equal to or greater than the value of 1.98 is 0.271. Therefore, the observed pN/pS value for MC4R is not significantly different from expectations based on a purely random distribution of mutations.

Although the observed pN/pS value for MC4R is not different from purely random expectations, the occurrence of only two synonymous substitutions, and at such low frequency, remains puzzling. Certainly, continuous negative selection at linked non-synonymous sites will maintain a low level of diversity at synonymous sites. However, negative selection alone may not be sufficient to explain the observed data; another factor may be codon usage bias. As demonstrated by the FMutSel analyses, the majority of mutations at silent sites are expected to be either highly advantageous or highly disadvantageous, with average estimated selection coefficients for MC4R much larger than those for typical genes. Thus, selection will either quickly remove or fix a silent mutation, which may help explain the lack of silent polymorphisms.

Finally, non-synonymous mutations which have been identified as loss of function mutations tend to occur preferentially at class 0 sites (under very strong purifying selection) of the M3(discrete) model (Supplemental Table 2). Moreover, non-synonymous mutations that do not have an effect on MC4R function in the in vitro assay tend to occur preferentially at sites that are under less stringent constraints (i.e., weaker purifying selection). This correlation may be useful in assessing the potential functional consequences of newly identified mutations, and should assist further studies of the epidemiology of the MC4R locus.

In conclusion, MC4R has experienced strong negative selection throughout vertebrate evolution and an apparent increase in selective constraint along the primate phylogeny. Moreover, humans have a significant paucity of diversity at MC4R, which is not observed in chimpanzees. This suggests that the magnitude of negative selection has continued to increase along the human lineage. Further, many of the observed non-synonymous mutations reduce MC4R protein function, suggesting a detriment to an individual’s fitness. Given these combined observations we conclude that mutations at MC4R, across vertebrates and especially in primates and humans, are overwhelmingly deleterious to protein function, resulting in prominent sequence conservation and significantly low levels of diversity in humans.